centering and interpolation of the result

proxtoolbox/Problems/OrbitalTomog/Graphics/Phase_graphics.py

@@ -47,10 +47,6 @@ def Phase_graphics(config, output):
     # iter = output['iter']
     change = output['change']
 
-    #    if 'time' in output:
-    #        time = output['time']
-    #    else:
-    #        time=999
 
     f, ((ax1, ax2), (ax3, ax4)) = subplots(2, 2, figsize=(12,8))
 

proxtoolbox/Problems/OrbitalTomog/Graphics/__init__.py

 from .Phase_graphics import Phase_graphics
+from .array_tools import *
+from .complex_field_visualization import complex_to_rgb
 
-__all__ = ["Phase_graphics"]
+__all__ = ["Phase_graphics", "roll_to_pos", "complex_to_rgb",'fourier_interpolate']

proxtoolbox/Problems/OrbitalTomog/Graphics/array_tools.py

+import numpy as np
+from scipy.ndimage import measurements
+import scipy.fftpack as fft
+
+
+def shift_array(arr, dy, dx):
+    temp = np.roll(arr, (dy, dx), (0, 1))
+    return temp
+
+
+def roll_to_pos(arr, y=0, x=0, move_maximum=False, by_abs_val=True):
+    if move_maximum:
+        if by_abs_val or arr.dtype in [np.complex64, np.complex128]:
+            old = np.floor(measurements.maximum_position(abs(arr)))
+        else:
+            old = np.floor(measurements.maximum_position(arr))
+    else:
+        if by_abs_val or arr.dtype in [np.complex64, np.complex128]:
+            old = np.floor(measurements.center_of_mass(abs(arr)))
+        else:
+            old = np.floor(measurements.center_of_mass(arr))
+    temp = shift_array(arr, int(y - old[0]), int(x - old[1]))
+    return temp
+
+
+def shifted_fft(arr):
+    return fft.ifftshift(fft.fft2(fft.fftshift(arr)))
+
+
+def shifted_ifft(arr):
+    return fft.fftshift(fft.ifft2(fft.ifftshift(arr)))
+
+
+def fourier_interpolate(arr):
+    """
+    Interpolate 2d complex array by a factor 2
+
+    :param arr: numpy array
+    :return: interpolated 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)

proxtoolbox/Problems/OrbitalTomog/coronene_config.py

@@ -143,14 +143,13 @@ new_config = {
     'verbose': 1,  # options are 0 or 1
     'graphics': 1,  # whether or not to display figures, options are 0 or 1.
     # default is 1.
-    'anim': 2,  # whether or not to disaply ``real time" reconstructions
+    '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',  # unless specified, a default
-    'dataprocessor_plotting': True
-    # plotting subroutine will generate
-    # the graphics.  Otherwise, the user
-    # can write their own plotting subroutine
+    'graphics_display': 'Phase_graphics',  # name of the plotting routine
+    'dataprocessor_plotting': True,
+    'interpolate_result': True
+
 
 }
 

proxtoolbox/Problems/OrbitalTomog/phase.py

@@ -156,15 +156,21 @@ class Phase(Problem):
         """
         Processes the solution and generates the output
         """
-        pass
+        # Center the solution (since position is a degree of freedom,
+        # and if desired, interpolate the results.
+        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
+            self.output[key] = Graphics.roll_to_pos(self.output[key], yc, xc) # then move center of mass.
+            # 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] = Graphics.fourier_interpolate(self.output[key])
 
     def show(self):
         """
         Generates graphical output from the solution
         """
-
         print("Calculation time:")
         print(self.elapsed_time)
         _graphics = getattr(Graphics, self.config['graphics_display'])
-        # print("Retrieved graphics function")
         _graphics(self.config, self.output)