Commit ddd71f05 authored by Matthijs's avatar Matthijs
Browse files

move tools for orbital imaging code to utilities folder

parent aaa0521d
......@@ -4,3 +4,4 @@ docs/Manual/autotoc/
_build
build
dist
proxtoolbox/Algorithms/samsara
\ No newline at end of file
......@@ -31,6 +31,7 @@
from matplotlib.pyplot import subplots, show
import numpy as np
from .complex_field_visualization import complex_to_rgb
from .stack_viewer import XYZStackViewer, SingleStackViewer
def Phase_graphics(config, output):
......@@ -84,3 +85,8 @@ def Phase_graphics(config, output):
bx4.set_ylabel('$||x^{2k+1}-x^{2k}||$')
show()
def Phase_graphics_3d(config, output):
# TODO: do something smart
return 0
from .Phase_graphics import Phase_graphics
from .array_tools import *
from proxtoolbox.Utilities.OrbitalTomog.array_tools import *
from .complex_field_visualization import complex_to_rgb
from .stack_viewer import XYZStackViewer, SingleStackViewer
__all__ = ["Phase_graphics", "roll_to_pos", "complex_to_rgb",'fourier_interpolate']
__all__ = ["Phase_graphics", "roll_to_pos", "complex_to_rgb", 'fourier_interpolate',
'XYZStackViewer', 'SingleStackViewer']
......@@ -4,7 +4,7 @@ from numpy import max, min
from matplotlib.widgets import Slider
from warnings import warn
good_backends = ['Qt5Agg','TkAgg','Qt4Agg']
good_backends = ['Qt5Agg', 'TkAgg', 'Qt4Agg']
class SingleStackViewer:
......
from .orbitaltomog_data_processor import shifted_fft, shifted_ifft, support_from_autocorrelation, bin_2d_array
from .orbitaltomog_data_processor import support_from_autocorrelation
from proxtoolbox.Utilities.OrbitalTomog.array_tools import shifted_fft, shifted_ifft
from proxtoolbox.Utilities.OrbitalTomog.binning import bin_2d_array, bin_array
import numpy as np
from skimage.io import imread
from .Graphics.stack_viewer import XYZStackViewer
......
......@@ -4,6 +4,8 @@ from scipy.io import loadmat
from scipy.ndimage import binary_dilation
import matplotlib.pyplot as plt
import imageio
from proxtoolbox.Utilities.OrbitalTomog.binning import bin_2d_array
from proxtoolbox.Utilities.OrbitalTomog.array_tools import shifted_ifft, shifted_fft
def data_processor(config):
......@@ -108,54 +110,54 @@ 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 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,
......
......@@ -102,7 +102,7 @@ new_config = {
'beta_switch': 30, # iteration at which beta moves from beta_0 -> beta_max
'lambda_switch': 100,
'sparsity_parameter': 40,
'sparsity_parameter': 100,
'use_sparsity_with_support': False,
'symmetry_type': 1, # -1 for antissymmetric functions, 1 for symmetric ones.
'symmetry_axis': -1, # which axis is symmetric. (mirror plane perpendicular to this axis)
......@@ -154,17 +154,10 @@ new_config = {
'anim': 2, # whether or not to display ``real time" reconstructions
# options are 0=no, 1=yes, 2=make a movie
# default is 1.
'graphics_display': 'Phase_graphics', # name of the plotting routine
'graphics_display': 'Phase_graphics_3d', # name of the plotting routine
'dataprocessor_plotting': True,
'interpolate_result': True,
'zoomin_on_result': True
}
# =====================================================================
# 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.
......@@ -160,6 +160,7 @@ class Phase(Problem):
"""
# Center the solution (since position is a degree of freedom,
# and if desired, interpolate the results.
# TODO: make compatible with non-2D data (e.g. 3D)
yc, xc = int(self.config['Ny'] / 2), int(self.config["Nx"] / 2)
for key in ['u', 'u1', 'u2']:
self.output[key] = Graphics.roll_to_pos(self.output[key], yc, xc, move_maximum=True) # first move maximum
......
......@@ -23,26 +23,30 @@ def roll_to_pos(arr, y=0, x=0, move_maximum=False, by_abs_val=True):
return temp
def shifted_fft(arr):
return fft.ifftshift(fft.fft2(fft.fftshift(arr)))
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
"""
def shifted_ifft(arr):
return fft.fftshift(fft.ifft2(fft.ifftshift(arr)))
return fft.ifftshift(fft.fftn(fft.fftshift(arr, axes=axes), axes=axes), axes=axes)
def fourier_interpolate(arr):
def shifted_ifft(arr, axes=None):
"""
Interpolate 2d complex array by a factor 2
Combined fftshift and fft routine, based on scipy.fftpack
Args:
arr: numpy array
axes: identical to argument for scipy.fftpack.fft
:param arr: numpy array
:return: interpolated array
Returns:
transformed array
"""
ny, nx = arr.shape
assert ny == nx
pd = int(ny / 2)
_arr = roll_to_pos(arr, pd, pd, move_maximum=True)
_arr = roll_to_pos(_arr, pd, pd)
fd = shifted_fft(arr)
tmp = np.pad(fd, ((pd, pd), (pd, pd)), mode='constant')
return shifted_ifft(tmp)
return fft.fftshift(fft.ifftn(fft.ifftshift(arr, axes=axes), axes=axes), axes=axes)
import numpy as np
def bin_array(arr: np.ndarray, new_shape: any, pad_zeros=True) -> np.ndarray:
"""
Reduce the size of an array by binning
:param arr: original
:param new_shape: tuple which must be an integer divisor of the original shape, or integer to bin by that factor
:return: new array
"""
# make tuple with new shape
if type(new_shape) == int: # binning factor is given
_shape = tuple([i // new_shape for i in arr.shape])
binfactor = tuple([new_shape for i in _shape])
else:
_shape = new_shape
binfactor = tuple([s // _shape[i] for i, s in enumerate(arr.shape)])
# determine if padding is needed
padding = tuple([(0, (binfactor[i] - s % binfactor[i]) % binfactor[i]) for i, s in enumerate(arr.shape)])
if pad_zeros and np.any(np.array(padding) != 0):
_arr = np.pad(arr, padding, mode='constant', constant_values=0) # pad array
_shape = tuple([s//binfactor[i] for i, s in enumerate(_arr.shape)]) # update binned size due to padding
else:
_arr = arr # expected to fail if padding has non-zeros
# send to 2d or 3d padding functions
try:
if len(arr.shape) == 2:
out = bin_2d_array(_arr, _shape)
elif len(arr.shape) == 3:
out = bin_3d_array(_arr, _shape)
else:
raise NotImplementedError('Cannot only bin 3d or 2d arrays')
return out
except ValueError:
raise ValueError("Cannot bin data with this shape. Try setting pad_zeros=True, or change the binning.")
def bin_2d_array(arr: np.ndarray, new_shape: tuple) -> np.ndarray:
"""
bins a 2D numpy array
Args:
arr: input array to be binned
new_shape: shape after binning, must be an integer divisor of the original shape
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 bin_3d_array(arr: np.ndarray, new_shape: tuple) -> np.ndarray:
""""
bins a 3D numpy array
Args:
arr: input array to be binned
new_shape: shape after binning, must be an integer divisor of the original shape
Returns:
binned np array
"""
shape = (new_shape[0], arr.shape[0] // new_shape[0],
new_shape[1], arr.shape[1] // new_shape[1],
new_shape[2], arr.shape[2] // new_shape[2])
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=(5, 3, 1)) * binfactor
else:
return np.sum(arr.reshape(shape), axis=(5, 3, 1))
import numpy as np
from .array_tools import shifted_fft, shifted_ifft, roll_to_pos
def fourier_interpolate(arr: np.ndarray, factor: any = 2., **kwargs) -> np.ndarray:
"""
Interpolation by padding in the Fourier domain.
Keyword arguments get passed to sub-functionalities. Options are for example:
- center_data: whether to move the center of mass to the center of the array before interpolation,
this is generally a good idea
NOTE: testing by simple 2x interpolation and subsequent 2x2 binning shows that errors are introduced in the procedure
:param arr: numpy array, 2 or 3-dimensional
:param factor: tuple or float to indicate by which factor to interpolate, default to factor 2
:returns: interpolated array, cast to identical dtype as arr
"""
try:
if len(factor) == 2:
out = fourier_interpolate_2d_nonsquare(arr, factor=factor, **kwargs)
elif len(factor) == 3:
out = fourier_interpolate_3d_nonsquare(arr, factor=factor, **kwargs)
else:
raise Exception("This error should never raise")
except TypeError:
if len(arr.shape) == 2:
out = fourier_interpolate_2d(arr, factor=factor, **kwargs)
elif len(arr.shape) == 3:
out = fourier_interpolate_3d(arr, factor=factor, **kwargs)
else:
raise NotImplementedError('Can only interpolate 2d or 3d arrays as of now')
return out.astype(arr.dtype, casting='unsafe')
def fourier_interpolate_2d_nonsquare(arr: np.ndarray, factor: tuple = (1., 2.),
center_data: bool = False) -> np.ndarray:
"""
Interpolate 2d array (can be complex-valued, should be well-sampled)
:param arr: numpy array of rectangular size
:param factor: tuple of interpolation factor (2 for doubling of the number of sampling points)
:param center_data: automatic centering of the data
:return: interpolated array
"""
ny, nx = arr.shape
pdy, pdx = int((ny / 2) * (factor[0] - 1)), int((nx / 2) * (factor[1] - 1))
if center_data:
_arr = roll_to_pos(arr, ny // 2, nx // 2, move_maximum=True)
_arr = roll_to_pos(_arr, ny // 2, nx // 2, )
fd = shifted_fft(arr)
tmp = np.pad(fd, ((pdy, pdy), (pdx, pdx)), mode='constant')
return shifted_ifft(tmp)
def fourier_interpolate_3d_nonsquare(arr: np.ndarray, factor: tuple = (1., 2., 3.)) -> np.ndarray:
"""
Interpolate 3d data (can be complex-valued, should be oversampled)
:param arr: dataset to be interpolated.
:param factor: interpolation factor (2 for doubling of the number of sampling points)
:return: interpolated data
"""
ny, nx, nz = arr.shape
pdy, pdx, pdz = int((ny / 2) * (factor[0] - 1)), int((nx / 2) * (factor[1] - 1)), int((nz / 2) * (factor[2] - 1))
fd = shifted_fft(arr)
tmp = np.pad(fd, ((pdy, pdy), (pdx, pdx), (pdz, pdz)), mode='constant')
return shifted_ifft(tmp)
def fourier_interpolate_2d(arr: np.ndarray, factor: float = 2., center_data=False) -> np.ndarray:
"""
Interpolate 2d array (can be complex-valued, should be oversampled)
:param arr: numpy array of rectangular size
:param factor: interpolation factor (2 for doubling of the number of sampling points)
:param center_data: automatic centering of the data
:return: interpolated array
"""
ny, nx = arr.shape
assert ny == nx, "Accepts only rectangular arrays currently"
pd = int((ny / 2) * (factor - 1))
if center_data:
_arr = roll_to_pos(arr, ny // 2, nx // 2, move_maximum=True)
_arr = roll_to_pos(_arr, ny // 2, nx // 2, )
fd = shifted_fft(arr)
tmp = np.pad(fd, ((pd, pd), (pd, pd)), mode='constant')
return shifted_ifft(tmp)
def fourier_interpolate_3d(arr: np.ndarray, center_data: bool = False, factor: float = 2.) -> np.ndarray:
"""
Interpolate 3d data (can be complex-valued, should be oversampled)
:param arr: dataset to be interpolated.
:param center_data: not implemented: would enable automatic centering of the data
:param factor: interpolation factor (2 for doubling of the number of sampling points)
:return: interpolated data
"""
ny, nx, nz = arr.shape
pd = int((ny / 2) * (factor - 1))
if center_data:
raise NotImplementedError
else:
_arr = arr
fd = shifted_fft(arr)
tmp = np.pad(fd, ((pd, pd), (pd, pd), (pd, pd)), mode='constant')
return shifted_ifft(tmp)
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