Skip to content

GitLab

Commit 1939a038 authored by Matthijs's avatar Matthijs
Browse files

partial merge from orbital_tomography-3d branch

parent 56ec4dcb
Hide whitespace changes
Inline Side-by-side
proxtoolbox/Problems/OrbitalTomog/Graphics/Phase_graphics.py
View file @ 1939a038
......@@ -35,47 +35,34 @@ from .stack_viewer import XYZStackViewer, SingleStackViewer
def Phase_graphics(config, output):
# algortihm=config['algorithm']
# beta0 = config['beta_0']
# beta_max = config['beta_max']
# u_0 = config['u_0']
if output['u1'].ndim == 2:
u = output['u1']
u2 = output['u2']
else:
u = output['u1'][:, :, 0]
u2 = output['u2'][:, :, 0]
# iter = output['iter']
change = output['change']
f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2, figsize=(12,8))
f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2, figsize=(9, 6))
im = ax1.imshow(np.abs(u), cmap='gray')
f.colorbar(im, ax=ax1)
ax1.set_title('best approximation amplitude - physical constraint satisfied')
ax1.set_title('best approximation amplitude: \nphysical constraint satisfied')
im = ax2.imshow(complex_to_rgb(u))
# im = ax2.imshow(np.real(u), cmap='gray')
# f.colorbar(im, ax=ax2)
ax2.set_title('best approximation phase - physical constraint satisfied')
ax2.set_title('best approximation phase: \nphysical constraint satisfied')
im = ax3.imshow(np.abs(u2), cmap='gray')
f.colorbar(im, ax=ax3)
ax3.set_title('best approximation amplitude - Fourier constraint satisfied')
ax3.set_title('best approximation amplitude: \nFourier constraint satisfied')
im = ax4.imshow(complex_to_rgb(u2))
# im = ax4.imshow(np.real(u2), cmap='gray')
# f.colorbar(im, ax=ax4)
ax4.set_title('best approximation amplitude - Fourier constraint satisfied')
ax4.set_title('best approximation amplitude: \nFourier constraint satisfied')
f.tight_layout()
g, ((bx1, bx2), (bx3, bx4)) = subplots(2, 2, figsize=(12,8))
g, ((bx1, bx2), (bx3, bx4)) = subplots(2, 2, figsize=(9, 6))
im = bx1.imshow(np.abs(u), cmap='gray')
f.colorbar(im, ax=bx1)
bx1.set_title('best approximation amplitude - physical constraint satisfied')
bx1.set_title('best approximation amplitude: \nphysical constraint satisfied')
im = bx2.imshow(np.real(u), cmap='gray')
f.colorbar(im, ax=bx2)
bx2.set_title('best approximation phase - physical constraint satisfied')
bx2.set_title('best approximation phase: \nphysical constraint satisfied')
bx3.semilogy(change)
bx3.set_xlabel('iteration')
bx3.set_ylabel('$||x^{2k+2}-x^{2k}||$')
......@@ -83,7 +70,7 @@ def Phase_graphics(config, output):
bx4.semilogy(output['gap'])
bx4.set_xlabel('iteration')
bx4.set_ylabel('$||x^{2k+1}-x^{2k}||$')
f.tight_layout()
show()
......@@ -96,13 +83,15 @@ def Phase_graphics_3d(config, output):
phys_constraint_plot = XYZStackViewer(output['u1'])
fourier_constraint_plot = XYZStackViewer(output['u2'])
convergence_plots, ax = subplots(1, 2, figsize=(6, 3.5))
convergence_plots, ax = subplots(1, 2, figsize=(7, 3))
ax[0].semilogy(output['change'])
ax[0].set_title('Change')
ax[0].set_xlabel('iteration')
ax[0].set_ylabel('$||x^{2k+2}-x^{2k}||$')
if "diagnostic" in config:
ax[1].semilogy(output['gap'])
ax[1].set_xlabel('iteration')
ax[1].set_title('Gap')
ax[1].set_ylabel('$||x^{2k+1}-x^{2k}||$')
show()
......
proxtoolbox/Problems/OrbitalTomog/Graphics/stack_viewer.py
View file @ 1939a038
......@@ -52,7 +52,7 @@ class SingleStackViewer:
class XYZStackViewer:
def __init__(self, volume, limit_sliders=(0, 0, 0), cmap: str = 'seismic', clim: tuple = (None, None),
data_transform: callable = None):
data_transform: callable = None, name: str = None):
"""
Dynamic matplotlib plot showing slices out of a 3d dataset
......@@ -62,10 +62,17 @@ class XYZStackViewer:
:param clim: tuple with the min-max values of the colorscale
:param data_transform: simple transformation function to apply to data before plot.
defaults to abs() for complex data
:param name: Set the name of the figure to a string.
"""
if get_backend() not in good_backends:
warn(Warning('Current matplotlib backend may not allow for optimal funcionality! Use, e.g., Qt'))
if volume.ndim == 4 and volume.shape[-1] in [3,4]:
rgb_plot = True
else:
# assert volume.ndim == 3, 'Functionality is designed for 3D arrays'
rgb_plot = False
if cmap is not None:
show_colorbar = True
else:
......@@ -76,6 +83,7 @@ class XYZStackViewer:
elif iscomplexobj(volume):
self.cast_fn = complex_to_rgb
show_colorbar = False
rgb_plot = True
else:
def dummy(array):
return array
......@@ -91,10 +99,13 @@ class XYZStackViewer:
self.maxval = max(abs(volume))
self.minval = -1*self.maxval
else:
self.maxval = max(volume)
self.minval = min([0, min(volume)])
if not rgb_plot:
self.maxval = max(self.cast_fn(volume))
self.minval = min([0, min(self.cast_fn(volume))])
else:
self.minval, self.maxval = 0, 255
self.fig, self.ax = plt.subplots(1, 3, figsize=(10, 4))
self.fig, self.ax = plt.subplots(1, 3, figsize=(10, 4), num=name)
im = self.ax[0].imshow(self.cast_fn(volume[self.indices[0], :, :]), cmap=cmap, vmin=self.minval,
vmax=self.maxval)
im = self.ax[1].imshow(self.cast_fn(volume[:, self.indices[0], :]), cmap=cmap, vmin=self.minval,
......@@ -105,7 +116,7 @@ class XYZStackViewer:
self.ax[2].set_yticks([])
self.fig.subplots_adjust(bottom=0.1, right=0.85)
if show_colorbar:
if show_colorbar and not rgb_plot:
cbar_ax = self.fig.add_axes([0.87, 0.2, 0.01, 0.6])
plt.colorbar(im, cax=cbar_ax)
......@@ -128,7 +139,6 @@ class XYZStackViewer:
self.fig.show()
def update_1(self, n):
"""Update subfigure 1 on change of the slider """
self.indices[0] = int(n)
......
proxtoolbox/Problems/OrbitalTomog/phase.py
View file @ 1939a038
......@@ -7,6 +7,7 @@ from proxtoolbox.ProxOperators.proxoperators import ProxOperator
from proxtoolbox.Problems.OrbitalTomog import Graphics
from numpy.linalg import norm
from numpy import square, sqrt
from numpy.fft import fftshift
from proxtoolbox.Utilities.OrbitalTomog import interpolation, array_tools, binning
......@@ -103,7 +104,7 @@ class Phase(Problem):
# end
else: # i.e. formulation=='product space'
proj_1 = self.config['proxoperators'][0](self.config) # , self.back_end)
u_1 = proj_1.work(self.config['u_0'])
u_1 = proj_1.work(self.config['u_0'].copy())
proj_2 = self.config['proxoperators'][1](self.config) # , self.back_end)
u_2 = proj_2.work(u_1)
......@@ -126,7 +127,7 @@ class Phase(Problem):
# the second tolerance relative to the order of magnitude of the metric
self.config['TOL2'] = self.config['data_ball'] * 1e-15
# self.config['proxoperators']
self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config)#, self.back_end)
self.algorithm = getattr(Algorithms, self.config['algorithm'])(self.config)
def _presolve(self):
"""
......@@ -148,22 +149,23 @@ class Phase(Problem):
"""
Processes the solution and generates the output
"""
# Center the solution (since position is a degree of freedom,
# and if desired, interpolate the results.
# Center the solution (since position is a degree of freedom, and if desired, interpolate the results.
center = tuple([s//2 for s in self.config['u'].shape])
for key in ['u', 'u1', 'u2']:
if 'fourier_shift_arrays' in self.config and self.config['fourier_shift_arrays']:
self.output[key] = fftshift(self.output[key])
self.output[key] = array_tools.roll_to_pos(self.output[key], pos=center, move_maximum=True)
self.output[key] = array_tools.roll_to_pos(self.output[key], pos=center)
# This sequence will work for objects *with a small support* even if they lie over the edge of the array
if 'interpolate_result' in self.config and self.config['interpolate_result']:
self.output[key] = interpolation.fourier_interpolate(self.output[key])
if 'zoomin_on_result' in self.config and self.config['zoomin_on_result']:
# TODO: if some support given, use bounding box of the support as zoom in region
if self.output[key].ndim == 2:
zmy, zmx = self.output[key].shape # self.config['Ny'] // 4, self.config["Nx"] // 4
zmy, zmx = tuple([s//4 for s in self.output[key].shape])
self.output[key] = self.output[key][zmy:-zmy, zmx:-zmx]
elif self.output[key].ndim == 3:
zmy, zmx, zmz = self.output[key].shape
# (self.config['Ny'] // 4, self.config["Nx"] // 4, self.config['Nz'] // 4)
zmy, zmx, zmz = tuple([s//4 for s in self.output[key].shape])
self.output[key] = self.output[key][zmy:-zmy, zmx:-zmx, zmz:-zmz]
def show(self):
......@@ -173,7 +175,7 @@ class Phase(Problem):
print("Calculation time:")
print(self.elapsed_time)
_graphics = getattr(Graphics, self.config['graphics_display'])
_graphics(self.config, self.output)
self.output['plots'] = _graphics(self.config, self.output)
def save(self, save_dir=None):
"""
......
proxtoolbox/Problems/OrbitalTomog/Coronene_demo.py → proxtoolbox/Problems/OrbitalTomog/planar_molecule/Coronene_demo.py
View file @ 1939a038
from proxtoolbox.Problems.OrbitalTomog import coronene_config # base config
from proxtoolbox.Problems.OrbitalTomog import orbitaltomog_data_processor # extends config
from proxtoolbox.Problems.OrbitalTomog.planar_molecule import coronene_config, orbitaltomog_data_processor
from proxtoolbox.Problems.OrbitalTomog.phase import Phase
# sys.path.append('../proxtoolbox/Problems/Phase')
......
proxtoolbox/Problems/OrbitalTomog/coronene_config.py → proxtoolbox/Problems/OrbitalTomog/planar_molecule/coronene_config.py
View file @ 1939a038
......@@ -27,9 +27,8 @@ new_config = {
# Orbital imaging options are: '2D ARPES', '3D ARPES', '2D time', and '3D time''
'experiment': '2D ARPES',
# Next we move to things that most of our users will know
# better than we will. Some of these may be overwritten in the
# data processor file which the user will most likely write.
# Next we move to things that most of our users will know better than we will. Some of these may be overwritten
# in the data processor file which the user will most likely write.
# Are the measurements in the far field or near field?
# Options are: 'far field' or 'near field',
'distance': 'far field',
......@@ -59,32 +58,8 @@ new_config = {
# Algorithm:
'algorithm': 'DRl', # used to be 'Projection',
'numruns': 1, # the only time this parameter will
# be different than 1 is when we are
# benchmarking...not something a normal user
# would be doing.
# # The following are parameters specific to RAAR, HPR, and HAAR that the
# # user should be able to set/modify. Surely
# # there will be other algorithm specific parameters that a user might
# # want to play with. Don't know how best
# # to do this. Thinking of a GUI interface, we could hard code all the
# # parameters the user might encounter and have the menu options change
# # depending on the value of the 'method field.
# # do different things depending on the chosen algorithm:
# if(strcmp('method,'RAAR')||strcmp('method,'AP')||...
# strcmp('method,'HPR')||strcmp('method,'HAAR'))
# the following just points this driver to a patch that communicates
# the parameters defined at this level to the structures used in the
# algorithms developed by Russell.
# 'problem_family' : 'Phase',
# else # moreregularization parameters for Thorsten's algorithms:
# prbl' : complete_itreg_par(prbl),
# This is just a patch to Thorsten's tools. May want to change this
# later
# 'problem_family='Hohage',
# end
# if(strcmp('problem_family,'Phase'))
# be different than 1 is when we are benchmarking... not something a normal user would be doing.
# maximum number of iterations and tolerances
'MAXIT': 500,
'TOL': 1e-10,
......@@ -108,39 +83,12 @@ new_config = {
# need to discuss how/whether the user should
# put in information about the noise
'data_ball': .999826,
# 'data_ball' : .9998261e-0,
# the above is the percentage of the gap
# between the measured data and the
# initial guess satisfying the
# qualitative constraints. For a number
# very close to one, the gap is not expected
# to improve much. For a number closer to 0
# the gap is expected to improve a lot.
# Ultimately the size of the gap depends
# on the inconsistency of the measurement model
# with the qualitative constraints.
# elseif(strcmp('problem_family,'Hohage'))
# 'alpha0' : 1e4,
# 'N_CG' : 70,
# end
# the above is the percentage of the gap between the measured data and the
# initial guess satisfying the qualitative constraints. For a number
# very close to one, the gap is not expected to improve much. For a number closer to 0
# the gap is expected to improve a lot. Ultimately the size of the gap depends
# on the inconsistency of the measurement model with the qualitative constraints.
# # ==========================================
# # parameters for plotting and diagnostics
# # ==========================================
# 'plotWhat.n1=2,
# 'plotWhat.n2=3,
# 'plotWhat.plots' : 'PYWpyw',
# 'verbose' : 1, # options are 0 or 1
# 'graphics' : 1, # whether or not to display figures, options are 0 or 1.
# # default is 1.
# 'anim' : 1, # whether or not to disaply ``real time" reconstructions
# # options are 0=no, 1=yes, 2=make a movie
# # default is 1.
# 'graphics_display' : [], # unless specified, a default
# # plotting subroutine will generate
# # the graphics. Otherwise, the user
# # can write their own plotting subroutine
#
# # ==========================================
# # parameters for plotting and diagnostics
# # ==========================================
......@@ -158,10 +106,3 @@ new_config = {
}
# =====================================================================
# Technical/software specific parameters
# =====================================================================
# Given the parameter values above, the following technical/algorithmic
# parameters are automatically set. The user does not need to know
# about these details, and so probably these parameters should be set in
# a module one level below this one.
proxtoolbox/Problems/OrbitalTomog/orbitaltomog_data_processor.py → proxtoolbox/Problems/OrbitalTomog/planar_molecule/orbitaltomog_data_processor.py
View file @ 1939a038
import numpy as np
import numpy.fft as fft
from scipy.io import loadmat
from scipy.ndimage import binary_dilation
import matplotlib.pyplot as plt
import imageio
from skimage.io import imread
from proxtoolbox.Utilities.OrbitalTomog.binning import bin_2d_array
from proxtoolbox.Utilities.OrbitalTomog.array_tools import shifted_ifft, shifted_fft
......@@ -12,17 +11,18 @@ def data_processor(config):
if config['data_filename'][1] == ':':
path_to_file = config['data_filename']
else:
path_to_file = '../../../InputData/OrbitalTomog/' + config['data_filename']
path_to_file = '../../../../InputData/OrbitalTomog/' + config['data_filename']
try:
if config['data_filename'][-4:] == '.mat': # for the old coronene simulations by W. Bennecke
inp_data = loadmat(path_to_file)
inp = inp_data["I2"]
else:
inp = imageio.imread(path_to_file)
inp = imread(path_to_file)
except FileNotFoundError:
print("Tried path %s, found nothing. "% path_to_file)
path_to_file = input('Please enter the path to the datafile: ')
inp = imageio.imread(path_to_file)
inp = imread(path_to_file)
ny, nx = inp.shape
config['data'] = abs(inp)
......@@ -51,7 +51,7 @@ def data_processor(config):
# Object support determination
try:
config['support'] = imageio.imread(config['support_filename'])
config['support'] = imread(config['support_filename'])
except KeyError: # 'support filename does not exist, so define a support here'
if 'threshold for support' not in config:
config['threshold for support'] = 0.1
......@@ -71,21 +71,21 @@ def data_processor(config):
# Autocorrelation
config['threshold for support'] = 0.1
autocorrelation = abs(fft.fftshift(fft.ifft2(fft.ifftshift(inp))))
autocorrelation = abs(shifted_ifft(inp))
config['abs_illumination'] = np.ones_like(config['support'])
# Initial guess
ph_init = 2 * np.pi * np.random.rand(ny, nx)
config['u_0'] = inp * np.exp(1j * ph_init)
config['u_0'] = np.fft.fftn(config['u_0'])
if ('use_sparsity_with_support' in config
and config['use_sparsity_with_support']
and 'sparsity_support' not in config):
ac = abs(fft.fftshift(fft.ifft2(fft.ifftshift(inp))))
config['sparsity_support'] = binary_dilation(ac > 0.01 * np.amax(ac)).astype(np.uint)
config['sparsity_support'] = binary_dilation(autocorrelation > 0.01 * np.amax(autocorrelation)).astype(np.uint)
if config['dataprocessor_plotting']:
plt.figure(figsize=(15, 4))
plt.figure(figsize=(10, 3.5))
plt.subplot(131)
plt.imshow(inp)
plt.colorbar()
......@@ -110,55 +110,6 @@ def data_processor(config):
return config
#
# def bin_2d_array(arr, new_shape):
# """"
# bins a 2D numpy array
# Args:
# arr: input array to be binned
# new_shape: shape after binning
# Returns:
# binned np array
# """
# shape = (new_shape[0], arr.shape[0] // new_shape[0],
# new_shape[1], arr.shape[1] // new_shape[1])
# if np.any(np.isnan(arr)):
# binfactor = 1
# for i, s in enumerate(arr.shape):
# binfactor *= new_shape[i] / s
# return np.nanmean(arr.reshape(shape), axis=(3, 1)) * binfactor
# else:
# return arr.reshape(shape).sum(-1).sum(1)
# def shifted_fft(arr, axes=None):
# """
# Combined fftshift and fft routine, based on scipy.fftpack
#
# Args:
# arr: numpy array
# axes: identical to argument for scipy.fftpack.fft
#
# Returns:
# transformed array
# """
#
# return fft.ifftshift(fft.fftn(fft.fftshift(arr, axes=axes), axes=axes), axes=axes)
#
#
# def shifted_ifft(arr, axes=None):
# """
# Combined fftshift and fft routine, based on scipy.fftpack
#
# Args:
# arr: numpy array
# axes: identical to argument for scipy.fftpack.fft
#
# Returns:
# transformed array
# """
# return fft.fftshift(fft.ifftn(fft.ifftshift(arr, axes=axes), axes=axes), axes=axes)
def support_from_autocorrelation(input_array: np.ndarray,
threshold: float = 0.1,
......
proxtoolbox/ProxOperators/P_M.py
View file @ 1939a038
......@@ -69,6 +69,7 @@ class P_M_masked(P_M):
-------
array_like - p_M: the projection IN THE PHYSICAL (time) DOMAIN
"""
constrained = super(P_M_masked, self).work(u)
fourier_space_iterate = np.fft.fft2(u, axes=None)
constrained = magproj(fourier_space_iterate, self.M)
update = np.where(self.mask, u, constrained)
return update
return np.fft.ifft2(update, axes=None)
proxtoolbox/ProxOperators/P_Sparsity.py
View file @ 1939a038
from numpy import zeros_like, unravel_index
from numpy import zeros_like, unravel_index, sum
import numpy as np
from .proxoperators import ProxOperator
......@@ -23,6 +23,7 @@ class P_Sparsity(ProxOperator):
if 'sparsity_support' in config:
self.support = config['sparsity_support'].real.astype(np.uint)
assert (sum(self.support) > 0), 'Invalid empty sparsity_support given'
else:
self.support = 1
......
proxtoolbox/Utilities/OrbitalTomog/__init__.py
View file @ 1939a038
from .array_tools import *
from .binning import *
from .interpolation import *
from .padding import *
proxtoolbox/Utilities/OrbitalTomog/array_tools.py
View file @ 1939a038
from numpy import roll, ndarray, floor, iscomplexobj, round
from numpy import roll, ndarray, floor, iscomplexobj, round, any, isnan, nan_to_num
from scipy.ndimage.measurements import maximum_position, center_of_mass
from scipy.fftpack import fftn, fftshift, ifftn, ifftshift
from warnings import warn
__all__ = ["shift_array", 'roll_to_pos', 'shifted_ifft', 'shifted_fft']
......@@ -40,6 +41,9 @@ def roll_to_pos(arr: ndarray, y: int = 0, x: int = 0, pos: tuple = None, move_ma
old = floor(center_of_mass(abs(arr)))
else:
old = floor(center_of_mass(arr))
if any(isnan(old)):
old = nan_to_num(old)
warn(Warning("Unexpected error in the calculation of the center of mass, casting NaNs to num"))
if pos is not None: # dimension-independent method
shifts = tuple([int(round(pos[i] - old[i])) for i in range(len(pos))])
dims = tuple([i for i in range(len(pos))])
......
proxtoolbox/Utilities/OrbitalTomog/padding.py 0 → 100644
View file @ 1939a038
from numpy import min, max, all, array, subtract, pad
__all__ = ['pad_to_square']
def pad_to_square(arr: array, new_shape: tuple = None, **kwargs) -> array:
"""
Given a numpy array of unequal dimensions,
pad with zeros until it all dimensions have equal length
:param arr: numpy array
:param new_shape: desired shape
:param kwargs: are passed on to np.pad
:return:
"""
old_shape = arr.shape
if max(old_shape) == min(old_shape) and new_shape is None:
return arr
if new_shape is None:
new_shape = tuple([max(old_shape) for dim in old_shape])
elif type(new_shape) == int:
new_shape = tuple(new_shape for dim in old_shape)
assert all(array(new_shape) >= array(old_shape)), 'New dimensions must be larger than old dimensions'
padding = subtract(new_shape, old_shape)
left_pad = padding // 2
right_pad = padding - left_pad
cv = kwargs.pop('constant_values', 0)
padmode = kwargs.pop('mode', 'constant')
return pad(arr, tuple([(left_pad[i], right_pad[i]) for i in range(len(left_pad))]),
mode=padmode, constant_values=cv, **kwargs)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment