Fixes a bug, which caused slight phase gradients if a ROI was used

docs/tutorial.rst

@@ -425,21 +425,52 @@ The program now should start with the reconstruction and averaging and should ou
 .. code-block:: none
 
     Loading parameters from
-            directory-with-data-files/holoseries-param.json
+            /home/niermann/PycharmProjects/holoaverage/example/holoaverage-a1.json
+    Loading...
+            20 datasets, shape=(2048, 2048), dtype=int16
+            [00] GaN_holographic_focal_series/empty_001.dm3
+            [01] GaN_holographic_focal_series/empty_002.dm3
+            [02] GaN_holographic_focal_series/empty_003.dm3
+            [03] GaN_holographic_focal_series/empty_004.dm3
+            [04] GaN_holographic_focal_series/empty_005.dm3
+            [05] GaN_holographic_focal_series/empty_006.dm3
+            [06] GaN_holographic_focal_series/empty_007.dm3
+            [07] GaN_holographic_focal_series/empty_008.dm3
+            [08] GaN_holographic_focal_series/empty_009.dm3
+            [09] GaN_holographic_focal_series/empty_010.dm3
+            [10] GaN_holographic_focal_series/empty_011.dm3
+            [11] GaN_holographic_focal_series/empty_012.dm3
+            [12] GaN_holographic_focal_series/empty_013.dm3
+            [13] GaN_holographic_focal_series/empty_014.dm3
+            [14] GaN_holographic_focal_series/empty_015.dm3
+            [15] GaN_holographic_focal_series/empty_016.dm3
+            [16] GaN_holographic_focal_series/empty_017.dm3
+            [17] GaN_holographic_focal_series/empty_018.dm3
+            [18] GaN_holographic_focal_series/empty_019.dm3
+            [19] GaN_holographic_focal_series/empty_020.dm3
     Reconstructing...
             . . . . . . . . . . . . . . . . . . . .
     Optimizing after iteration 0
             [NN] sx[px] sy[px] tx[1/px] ty[1/px] def[nm] Ampl   Phase  Error
             [00]  0.000  0.000  0.00000  0.00000   0.000 1.0000 +0.802 2.403186e+09
-            [01]  0.000  0.000  0.00000  0.00000   0.000 1.0000 -0.672 2.365461e+09
-            [02]  0.000  0.000  0.00000  0.00000   0.000 1.0000 -1.478 2.599820e+09
+            [01]  0.000  0.000  0.00000  0.00000   0.000 1.0000 -0.672 2.365460e+09
+            [02]  0.000  0.000  0.00000  0.00000   0.000 1.0000 -1.478 2.599819e+09
+            [03]  0.000  0.000  0.00000  0.00000   0.000 1.0000 -2.849 2.534287e+09
+            [04]  0.000  0.000  0.00000  0.00000   0.000 1.0000 -1.102 2.451658e+09
+            [05]  0.000  0.000  0.00000  0.00000   0.000 1.0000 +0.163 2.365165e+09
+            [06]  0.000  0.000  0.00000  0.00000   0.000 1.0000 +0.439 2.400905e+09
+            [07]  0.000  0.000  0.00000  0.00000   0.000 1.0000 -0.745 2.470458e+09
+            [08]  0.000  0.000  0.00000  0.00000   0.000 1.0000 +0.303 2.364581e+09
+            [09]  0.000  0.000  0.00000  0.00000   0.000 1.0000 +1.830 2.513006e+09
+            [10]  0.000  0.000  0.00000  0.00000   0.000 1.0000 +1.208 2.377789e+09
+            [11]  0.000  0.000  0.00000  0.00000   0.000 1.0000 -0.781 2.302518e+09
     ...
 
 and so on. Eventually it should output something like:
 
 .. code-block:: none
 
-        Iteration   7: total error=7.502063e+07
+        Iteration   7: total error=7.501946e+07
 
 This error number should go down and converge to a stable value within the last iterations.
 

holoaverage/grid.py

@@ -177,13 +177,14 @@ class Grid(object):
     # Private members
     #   _realSampling    np.ndarray[ndim,ndim]: Real space sampling in for 2D: ((YY,YX),(XY,XX)) in nm
     #   _rcprSampling    np.ndarray[ndim,ndim]: Reciprocal space sampling in 2D: ((YY,YX),(XY,XX)) in 1/nm
+    #   _realOffset      np.ndarray[ndim]: Real space offset in nm
     #   _shape           tuple: Size of simulation grid, len=ndim
     #   _complexType     np.dtype: Complex type
     #   _floatType       np.dtype: Float type
     #   _planFFT         FFT plan (destroys input)
     #   _planIFFT        IFFT plan (destroys input)
     #   _fft_nthreads    Number of threads to use in FFT plans
-    def __init__(self, shape, realSampling=None, rcprSampling=None, dtype=np.complex64):
+    def __init__(self, shape, realSampling=None, rcprSampling=None, realOffset=None, dtype=np.complex64):
         """
         Initializes grid. Either realSampling or rcprSampling must be given.
 
@@ -198,7 +199,9 @@ def __init__(self, shape, realSampling=None, rcprSampling=None, dtype=np.complex
                 Float: Isotropic sampling in 1/nm
                 1D-array: Diagonal sampling in 1/nm
                 2D-array: Anisotropic sampling 1/nm
-            dtype
+            _realOffset:
+                2D-array: Offset 1/nm
+            d type
                 np.dtype: Data type to use
         """
         self._shape = tuple(shape)
@@ -219,6 +222,9 @@ def __init__(self, shape, realSampling=None, rcprSampling=None, dtype=np.complex
             self._setRcprSampling(rcprSampling)
         else:
             raise ValueError("Either 'realSampling' or 'rcprSampling' must be given.")
+        self._realOffset = np.zeros(2, dtype=float)
+        if realOffset is not None:
+            self._realOffset[...] = realOffset
         self._planFFT = None
         self._planIFFT = None
         self._planRFFT = None
@@ -229,14 +235,16 @@ def __getstate__(self):
                 'complexType': self._complexType,
                 'floatType': self._floatType,
                 'realSampling': self._realSampling,
-                'rcprSampling': self._rcprSampling}
+                'rcprSampling': self._rcprSampling,
+                'realOffset': self._realOffset}
 
     def __setstate__(self, state):
         self._shape = state["shape"]
         self._complexType = state["complexType"]
         self._floatType = state["floatType"]
         self._realSampling = state["realSampling"]
         self._rcprSampling = state["rcprSampling"]
+        self._realOffset = state["realOffset"]
         self._planFFT = None
         self._planIFFT = None
         self._planRFFT = None
@@ -262,7 +270,7 @@ def fromDataSet(dataset):
             dim_unit = None
         dim_scale = ScaleMatrix(len(dataset.shape))
         if 'dim_scale' in dataset.attrs:
-            dim_scale.set(dataset.attrs['dim_scale']).getMatrix()
+            dim_scale.setReverse(dataset.attrs['dim_scale']).getMatrix()
         else:
             dim_scale.set(np.ones(len(dataset.shape)))
 
@@ -281,7 +289,12 @@ def fromDataSet(dataset):
                     warnings.warn("Expected unit to be '1/nm' not '%s'" % dim_unit[0])
 
         if space >= 0:
-            return Grid(shape, realSampling=dim_scale, dtype=dtype)
+            if 'dim_offset' in dataset.attrs:
+                dim_offset = np.zeros(2, dtype=float)
+                dim_offset[...] = dataset.attrs["dim_offset"][::-1]
+            else:
+                dim_offset = None
+            return Grid(shape, realSampling=dim_scale, realOffset=dim_offset, dtype=dtype)
         else:
             dim_scale.set(dim_scale.getMatrix() * np.atleast_1d(shape))
             return Grid(shape, rcprSampling=dim_scale, dtype=dtype)
@@ -314,6 +327,10 @@ def _setRcprSampling(self, value):
         self._rcprSampling.set(value)
         self._realSampling.set(np.linalg.inv(self._rcprSampling.getMatrix()))
 
+    @property
+    def realOffset(self):
+        return self._realOffset.copy()
+
     def hasDiagonalSampling(self):
         """Returns whether the sampling matrix is diagonal (sampling along grid axes)."""
         return self._realSampling.isDiagonal()
@@ -483,9 +500,11 @@ def getRealGrid(self):
             Td = np.diag(T)
             for i in range(ndim):
                 result[i] *= Td[i]
+                result[i] += self._realOffset[i]
         else:
             result = np.zeros((ndim,) + self._shape, dtype=self.floatType)
             for i in range(ndim):
                 for j in range(ndim):
-                    result[i] = result[i] + T[i, j] * idx[j]
+                    result[i] += T[i, j] * idx[j]
+                result[i] += self._realOffset[i]
         return result

holoaverage/main.py

@@ -140,9 +140,9 @@ def rescale_fourier(data, shape=None, out=None):
     elif out.shape != shape:
         raise ValueError("Output array shape must fit given shape.")
     out.fill(0)
-    out[shape[0] // 2 - half[0]:shape[0] // 2 + half[0], shape[1] // 2 - half[1]:shape[1] // 2 + half[1]] = \
-        np.fft.fftshift(np.fft.fft2(data))[data.shape[0] // 2 - half[0]:data.shape[0] // 2 + half[0],
-                                           data.shape[1] // 2 - half[1]:data.shape[1] // 2 + half[1]]
+    out[shape[0] // 2 - half[0]:shape[0] // 2 + half[0] + 1, shape[1] // 2 - half[1]:shape[1] // 2 + half[1] + 1] = \
+        np.fft.fftshift(np.fft.fft2(data))[data.shape[0] // 2 - half[0]:data.shape[0] // 2 + half[0] + 1,
+                                           data.shape[1] // 2 - half[1]:data.shape[1] // 2 + half[1] + 1]
     out[...] = np.fft.ifft2(np.fft.ifftshift(out))
     return out
 

holoaverage/rawalign.py

@@ -101,7 +101,7 @@ def rawAlign(series, qMax=None, verbose=0, roi=None, filter=None):
     return series
 
 
-def extractROI(series, roi, binning=1, verbose=0, dtype=None):
+def extractROI(series, roi, verbose=0, dtype=None):
     """
     Extracts ROI from a series the result as a new series.
 
@@ -112,17 +112,12 @@ def extractROI(series, roi, binning=1, verbose=0, dtype=None):
             The series to be reconstructed
         roi
             Region of interest in [px] (left, top, right, bottom)
-        binning
-            Binning (X,Y) to use
         dtype
             Type of ROI series, defaults to series type
     """
     if len(series.shape) != 2:
         raise ValueError("2D series expected.")
     shape = (roi[3] - roi[1], roi[2] - roi[0])
-    binning = np.diag(ScaleMatrix(2).set(binning).getMatrix()).astype(int)
-    if (shape[0] % binning[1]) != 0 or (shape[1] % binning[0]) != 0:
-        raise ValueError("ROI shape must be a multiple of binning.")
     if dtype is None:
         dtype = series.dtype
     raw_shift = series.attrs.get("raw_shift", np.zeros(series.indexShape + (2,), dtype=int))
@@ -131,14 +126,11 @@ def extractROI(series, roi, binning=1, verbose=0, dtype=None):
     dim_offset = series.attrs.get('dim_offset', np.zeros(2, dtype=float))
     if dim_scale.get() is not None:
         dim_offset = np.dot(dim_scale.getMatrix(), roi[0:2]) + dim_offset
-        dim_scale.set(dim_scale.getMatrix() * binning.astype(float))
-    old_binning = np.diag(ScaleMatrix(2).set(series.attrs.get('binning', 1)).getMatrix())
     overrides = {"dim_scale": dim_scale.getCompact(minRank=1),
-                 "dim_offset": dim_offset,
-                 "binning": old_binning * binning}
+                 "dim_offset": dim_offset}
 
     # Extracting
-    result = Series(series.indexShape, (shape[0] / binning[1], shape[1] / binning[0]), dtype=dtype)
+    result = Series(series.indexShape, shape, dtype=dtype)
     result.attrs.update(series.attrs)
     result.attrs.update(overrides)
     result.attrs['roi'] = roi
@@ -170,13 +162,11 @@ def extractROI(series, roi, binning=1, verbose=0, dtype=None):
         tmp[...] = np.mean(series[index].array)
         tmp[dy:dy + ny, dx:dx + nx] = series[index].array[sy:sy + ny, sx:sx + nx]
         data = DataSet(result.shape, dtype=dtype)
-        data.array[...] = 0
-        for i in range(binning[1]):
-            for j in range(binning[0]):
-                data.array[...] += tmp[i::binning[1], j::binning[0]]
+        data.array[...] = tmp
         data.attrs.update(series[index].attrs)
         data.attrs.update(overrides)
         data.attrs['roi'] = (_x + roi[0], _y + roi[1], _x + roi[2], _y + roi[3])
+        data.attrs['dim_offset'] = np.dot(dim_scale.getMatrix(), data.attrs['roi'][0:2]) + dim_offset
         result[index] = data
     if verbose > 0:
         print()

holoaverage/reconstruction.py

@@ -75,12 +75,13 @@ def __init__(self, grid, mask=None, carrier=None, shape=None, mtf=None):
         src_scale = grid.realSampling
         # Test for pixels out of desired region
         carrier_px = np.dot(src_scale, np.array((carrier[1], carrier[0]), dtype=float)) * np.array(self._grid.shape, dtype=float)
+        icarrier_px = np.floor(carrier_px).astype(int)
         ry, rx = self._grid.shape[0] // 2, self._grid.shape[1] // 2
         ny, nx = shape
         hy, hx = ny // 2, nx // 2
         dy, dx = 0, 0
-        sy = ry + int(carrier_px[0]) - hy
-        sx = rx + int(carrier_px[1]) - hx
+        sy = ry + icarrier_px[0] - hy
+        sx = rx + icarrier_px[1] - hx
         # print sx, sy, rx, ry, carrier_px, hx, hy
         if sy < 0:
             ny += sy
@@ -117,6 +118,7 @@ def __init__(self, grid, mask=None, carrier=None, shape=None, mtf=None):
         # Prepare FFT
         self._reco_grid = Grid(shape, realSampling=self._destScale, dtype=self._grid.complexType)
         self._reco_grid.prepareFFT(forwardFFT=False, backwardFFT=True, forwardRFFT=False, backwardRFFT=False)
+        self._reco_phase_wedge = np.dot(np.linalg.inv(src_scale), (carrier_px - icarrier_px) / np.array(self._grid.shape, dtype=float))
 
     def apply(self, data):
         if (data.shape != self._grid.shape):
@@ -128,12 +130,19 @@ def apply(self, data):
             tmp = np.fft.fftshift(fromHermite(self._grid.forwardRFFT(sub), sub.shape))
         else:
             tmp = sub
+
         # Extract and iFFT
         dx0, dy0, dx1, dy1 = self._destRect
         sx0, sy0, sx1, sy1 = self._srcRect
         side[dy0:dy1, dx0:dx1] = tmp[sy0:sy1, sx0:sx1] * self._mask
         result = DataSet(self._destShape, dtype=self._grid.complexType)
         result.array[...] = self._reco_grid.backwardFFT(np.fft.ifftshift(side))
+
+        # Correct for non-pixel aligned carrier
+        if not np.allclose(self._reco_phase_wedge, 0.0):
+            ry, rx = self._reco_grid.getRealGrid()
+            result.array[...] *= np.exp(-2.0j * np.pi * self._reco_phase_wedge[0] * ry) * np.exp(-2.0j * np.pi * self._reco_phase_wedge[1] * rx)
+
         result.attrs.update(data.attrs)
         result.attrs.update(self._reco_attrs)
         return result

tests/README.txt

@@ -1,4 +1,4 @@
-The tests require the example data to be present in the
+The tests in test_example.py require the example data to be present in the
 
     example/GaN_holographic_focal_series
 

tests/test_example.py

@@ -4,6 +4,10 @@
 import h5py
 import numpy as np
 
+import warnings
+warnings.filterwarnings("error",category=DeprecationWarning)
+warnings.filterwarnings("error",category=PendingDeprecationWarning)
+
 from holoaverage.main import holoaverage, rescale_fourier
 
 # Setup path to example data
@@ -18,9 +22,9 @@
 class TestExamples(unittest.TestCase):
     """Test whether the examples run."""
 
-    GAN_DATA_CONVERGENCE_256 = 1.337721e8           # Expected error for 256px object reconstruction
-    GAN_DATA_CONVERGENCE_256_FOCUS = 1.249292e8     # Expected error for 256px object reconstruction with defocus adjustment
-    GAN_EMPTY_CONVERGENCE_384 = 5.681450e10         # Expected error for 384px empty reconstruction
+    GAN_DATA_CONVERGENCE_256 = 58546370.75          # Expected error for 256px object reconstruction
+    GAN_DATA_CONVERGENCE_256_FOCUS = 5.467907e+07   # Expected error for 256px object reconstruction with defocus adjustment
+    GAN_EMPTY_CONVERGENCE_384 = 3.20882540e+10      # Expected error for 384px empty reconstruction
 
     # Test parameters
     VERBOSE = 0
@@ -48,11 +52,11 @@ def default_GaN_param(self):
             "object_names": os.path.join(EXAMPLE_PATH, "a1_%03d.dm3"),
             "object_first": 1,
             "object_last": 20,
-            "object_size": 256,
+            "object_size": 384,
             "empty_names": os.path.join(EXAMPLE_PATH, "empty_%03d.dm3"),
             "empty_first": 1,
             "empty_last": 20,
-            "empty_size": 384,
+            "empty_size": 512,
             "sampling": 0.00519824,
             "defocus_first": 20.0,
             "defocus_step": -2.0,
@@ -92,12 +96,12 @@ def test_GaN_example_other_sizes(self):
             empty_convergence = output["empty"].attrs["convergence"]
 
         # Test for empty convergence
-        self.assertTrue(np.allclose(empty_convergence[-1], self.GAN_EMPTY_CONVERGENCE_384 * (384.0/512.0)**2, rtol=0.1))
+        self.assertTrue(np.allclose(empty_convergence[-1], self.GAN_EMPTY_CONVERGENCE_384, rtol=0.1))
         delta = (empty_convergence[1:] - empty_convergence[:-1]) / empty_convergence[:-1]
         self.assertTrue(np.all(delta < +1e-4))
 
         # Test for data convergence
-        self.assertTrue(np.allclose(data_convergence[-1], self.GAN_DATA_CONVERGENCE_256 * (256.0/384.0)**2, rtol=0.1))
+        self.assertTrue(np.allclose(data_convergence[-1], self.GAN_DATA_CONVERGENCE_256, rtol=0.1))
         delta = (data_convergence[1:] - data_convergence[:-1]) / data_convergence[:-1]
         self.assertTrue(np.all(delta < +1e-4))
 
@@ -119,7 +123,7 @@ def test_GaN_example_disabled_pyfftw(self):
         self.assertTrue(np.all(delta < +1e-4))
 
         # Test for data convergence
-        self.assertTrue(np.allclose(data_convergence[-1], 1.337719e8, rtol=1e-4))
+        self.assertTrue(np.allclose(data_convergence[-1], self.GAN_DATA_CONVERGENCE_256, rtol=1e-4))
         delta = (data_convergence[1:] - data_convergence[:-1]) / data_convergence[:-1]
         self.assertTrue(np.all(delta < +1e-4))
 
@@ -184,7 +188,7 @@ def test_GaN_no_empty(self):
             data_convergence = output["data"].attrs["convergence"]
 
         # Test for data convergence
-        self.assertTrue(np.allclose(data_convergence[-1], 5.6592e14, rtol=1e-4))
+        self.assertTrue(np.allclose(data_convergence[-1], 4.88520723e+13, rtol=1e-4))
         delta = (data_convergence[1:] - data_convergence[:-1]) / data_convergence[:-1]
         self.assertTrue(np.all(delta < +1e-4))
 
@@ -203,7 +207,7 @@ def test_GaN_empty_from_raw(self):
             data_convergence = output["data"].attrs["convergence"]
 
         # Test for data convergence
-        self.assertTrue(np.allclose(data_convergence[-1], 5.6592e14, rtol=1e-4))
+        self.assertTrue(np.allclose(data_convergence[-1], 2.47313565e+14, rtol=1e-4))
         delta = (data_convergence[1:] - data_convergence[:-1]) / data_convergence[:-1]
         self.assertTrue(np.all(delta < +1e-4))
 
@@ -310,7 +314,7 @@ def test_GaN_synthetic_empty(self):
         np.testing.assert_allclose(np.angle(empty / new_empty)[10:-10, 10:-10], 0.0, atol=0.1)
 
         # Test for data convergence
-        self.assertTrue(np.allclose(data_convergence[-1], 5.1516e+11, rtol=1e-4))
+        self.assertTrue(np.allclose(data_convergence[-1], 2.25341306e+11, rtol=1e-4))
         delta = (data_convergence[1:] - data_convergence[:-1]) / data_convergence[:-1]
         self.assertTrue(np.all(delta < +1e-4))