updated orbital tomography demo

proxtoolbox/Problems/OrbitalTomog/Coronene_demo.py

-import sys
-sys.path.append('../../../')
-# from . 
-import coronene_test # base config
-import orbitaltomog_data_processor # extends config
-from phase import Phase
+from proxtoolbox.Problems.OrbitalTomog import coronene_config  # base config
+from proxtoolbox.Problems.OrbitalTomog import orbitaltomog_data_processor  # extends config
+from proxtoolbox.Problems.OrbitalTomog.phase import Phase
 # sys.path.append('../proxtoolbox/Problems/Phase')
 
-orbit_config = orbitaltomog_data_processor.data_processor(coronene_test.new_config)
+orbit_config = orbitaltomog_data_processor.data_processor(coronene_config.new_config)
 orbit = Phase(orbit_config)
 orbit.solve()
-orbit.show()
\ No newline at end of file
+orbit.show()

proxtoolbox/Problems/OrbitalTomog/coronene_config.py

@@ -8,7 +8,8 @@ new_config = {
     # Problem parameters
     # ==========================================
     # What is the name of the data file?
-    'data_filename': 'coronen_homo1_fourier_noise15.mat',
+    'data_filename': 'coronene_homo1.tif',  # In the directory '../../../InputData/OrbitalTomog/'
+    'from intensity data': False,  # File gives field amplitudes
 
     # What type of object are we working with?
     # Options are: 'phase', 'real', 'nonnegative', 'complex'
@@ -35,8 +36,8 @@ new_config = {
     'distance': 'far field',
 
     # What are the dimensions of the measurements?
-    'Nx': 128,
-    'Ny': 128,
+    'Nx': 64,
+    'Ny': 64,
 
     # 'fresnel_nr' : 0,
 
@@ -99,9 +100,9 @@ new_config = {
     'beta_switch': 30,  # iteration at which beta moves from beta_0 -> beta_max
     'lambda_switch': 100,
 
-    'sparsity_parameter': 16,
+    'sparsity_parameter': 40,
     'use_sparsity_with_support': False,
-    'symmetry_type': 1,  # -1 for antisymmetric functions, 1 for symmetric ones.
+    'symmetry_type': 1,  # -1 for antissymmetric functions, 1 for symmetric ones.
     'symmetry_axis': -1,  # which axis is symmetric. (mirror plane perpendicular to this axis)
 
     # parameter for the data regularization

proxtoolbox/Problems/OrbitalTomog/orbitaltomog_data_processor.py

@@ -13,18 +13,55 @@ def data_processor(config):
     else:
         inp = imageio.imread('../../../InputData/OrbitalTomog/' + config['data_filename'])
     ny, nx = inp.shape
-      
     config['data'] = abs(inp)
-    config['norm_data'] = np.sqrt(np.sum(config['data'] ** 2))
 
     # Keep the same resolution?
-    config['Ny'], config['Nx'] = ny, nx
+    if config['Ny'] is None and config['Nx'] is None:
+        config['Ny'], config['Nx'] = ny, nx
+    elif config['Ny'] == ny and config['Nx'] == nx:
+        pass
+    elif ny % config['Ny'] == 0 and nx % config["Nx"] == 0:
+        # binning must be done for the intensity-data, as that preserves the normalization
+        if not ('from intensity data' in config and config['from intensity data']):
+            config['data'] = np.sqrt(bin_2d_array(config['data'] ** 2, (config['Ny'], config["Nx"])))
+        else:
+            config['data'] = bin_2d_array(config['data'], (config['Ny'], config["Nx"]))
+    else:
+        raise ValueError('Incompatible values for Ny and Nx given in configuration dict')
+
+    # Calculate electric field
+    if 'from intensity data' in config and config['from intensity data']:
+        # avoid sqrt of negative numbers (due to background subtraction)
+        config['data'] = np.where(config['data'] > 0,
+                                  np.sqrt(abs(config['data'])),
+                                  np.zeros_like(config['data']))
+    config['norm_data'] = np.sqrt(np.sum(config['data'] ** 2))
+
+    # Object support determination
+    try:
+        config['support'] = imageio.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
+        config['support'] = support_from_autocorrelation(config['data'],
+                                                         threshold=config['threshold for support'],
+                                                         absolute_autocorrelation=True,
+                                                         binary_dilate_support=1)
+
+    if ('use_sparsity_with_support' in config
+            and config['use_sparsity_with_support']
+            and 'sparsity_support' not in config):
+        if 'threshold for support' not in config:
+            config['threshold for support'] = 0.1
+        config['sparsity_support'] = support_from_autocorrelation(config['data'],
+                                                                  threshold=config['threshold for support'],
+                                                                  binary_dilate_support=1)
 
     # Autocorrelation
-    config['threshold for support'] = 0.2
+    config['threshold for support'] = 0.1
     autocorrelation = abs(fft.fftshift(fft.ifft2(fft.ifftshift(inp))))
-    config['support'] = (autocorrelation > config['threshold for support'] * np.amax(autocorrelation)).astype(np.uint) # [0,0,1,1,0,..]
     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)
@@ -50,13 +87,96 @@ def data_processor(config):
         plt.colorbar()
         plt.title("Initial support")
         plt.show()
-    
+
     # Other settings
     config['fresnel_nr'] = 0
     config['FT_conv_kernel'] = 1
     config['use_farfield_formula'] = True
     config['magn'] = 1
-    config['data_sq'] = abs(config['data'])**2
+    config['data_sq'] = abs(config['data']) ** 2
     config['data_zeros'] = 1
-    
+
     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])
+    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,
+                                 relative_threshold: bool = True,
+                                 input_in_fourier_domain: bool = True,
+                                 absolute_autocorrelation: bool = True,
+                                 binary_dilate_support: int = 0) -> np.ndarray:
+    """
+    Determine an initial support from the autocorrelation of an object.
+
+    Args:
+        input_array: either the measured diffraction (arpes pattern) or a guess of the object
+        threshold: support is everywhere where the autocorrelation is higher than the threshold
+        relative_threshold: If true, threshold at threshold*np.amax(autocorrelation)
+        input_in_fourier_domain: False if a guess of the object is given in input_array
+        absolute_autocorrelation: Take the absolute value of the autocorrelation? (Generally a good idea)
+        binary_dilate_support: number of dilation operations to apply to the support.
+
+    Returns:
+        support array (same dimensions as input, dtype=np.int)
+    """
+    if not input_in_fourier_domain:
+        kspace = shifted_fft(input_array)
+    else:
+        kspace = input_array
+
+    autocorrelation = shifted_ifft(abs(kspace))  # Taking absolute value yields autocorrelation by conv. theorem)
+    if absolute_autocorrelation:
+        autocorrelation = abs(autocorrelation)
+    maxval = np.amax(autocorrelation)
+    if relative_threshold:
+        threshold_val = threshold * maxval
+    else:
+        threshold_val = threshold
+    support = (autocorrelation > threshold_val).astype(np.uint)
+    if binary_dilate_support > 0:
+        support = binary_dilation(support, iterations=binary_dilate_support).astype(np.uint)
+
+    return support

proxtoolbox/Problems/Phase/phase.py

@@ -259,7 +259,7 @@ class Phase(Problem):
             f = h5py.File('Phase_test_data/siemens_amplitude_u1_' + str(self.config['MAXIT']) + '.mat')
             u1 = f['u1'].value.view(np.complex)
         elif self.config['data_filename'] == 'Siemens_processor' and self.config[
-            'constraint'] == 'nonnegative and support':
+                'constraint'] == 'nonnegative and support':
             f = h5py.File('Phase_test_data/siemens_nonneg_u1_' + str(self.config['MAXIT']) + '.mat')
             u1 = f['u1'].value.view(np.float64)
         elif self.config['data_filename'] == 'Siemens_processor' and self.config['constraint'] == 'real and support':