ref: align the 2d qlsi code to 3d

qpretrieve/interfere/if_qlsi.py

@@ -6,7 +6,6 @@
 
 from .base import BaseInterferogram
 from ..fourier import get_best_interface
-from ..data_input import revert_to_data_input_format
 
 
 class QLSInterferogram(BaseInterferogram):
@@ -173,8 +172,14 @@ def run_pipeline(self, **pipeline_kws):
 
         # Obtain the phase gradients in x and y by taking the argument
         # of Hx and Hy.
-        px = unwrap_phase(np.angle(hx))
-        py = unwrap_phase(np.angle(hy))
+        # need to do this along the z axis, as skimage `unwrap_3d` does not
+        # work for our use-case
+        # todo: maybe use np.unwrap for the xy axes instead
+        px = np.zeros_like(hx)
+        py = np.zeros_like(hy)
+        for i, (_hx, _hy) in enumerate(zip(hx, hy)):
+            px[i] = unwrap_phase(np.angle(_hx))
+            py[i] = unwrap_phase(np.angle(_hy))
 
         # Determine the angle by which we have to rotate the gradients in
         # order for them to be aligned with x and y. This angle is defined
@@ -191,8 +196,8 @@ def run_pipeline(self, **pipeline_kws):
                           mode="constant", constant_values=0)
 
         # Perform rotation of the gradients.
-        rotated1 = rotate_noreshape(gradpad1, -angle)
-        rotated2 = rotate_noreshape(gradpad2, -angle)
+        rotated1 = rotate_noreshape(gradpad1, -angle, axes=(-1, -2))
+        rotated2 = rotate_noreshape(gradpad2, -angle, axes=(-1, -2))
 
         # Retrieve the wavefront by integrating the vectorial components
         # (integrate the total differential). This magical approach
@@ -217,10 +222,10 @@ def run_pipeline(self, **pipeline_kws):
         # Rotate the wavefront back and crop it so that the FOV matches
         # the input data.
         raw_wavefront = rotate_noreshape(
-            wfr, angle)[:, sx // 2:-sx // 2, sy // 2:-sy // 2]
+            wfr, angle, axes=(-1, -2))[:, sx // 2:-sx // 2, sy // 2:-sy // 2]
         # Multiply by qlsi pitch term and the scaling factor to get
         # the quantitative wavefront.
-        scaling_factor = self.fft_origin.shape[0] / wfr.shape[0]
+        scaling_factor = self.fft_origin.shape[-2] / wfr.shape[-2]
         raw_wavefront *= qlsi_pitch_term * scaling_factor
 
         self._phase = raw_wavefront / wavelength * 2 * np.pi
@@ -231,8 +236,8 @@ def run_pipeline(self, **pipeline_kws):
 
         self.pipeline_kws.update(pipeline_kws)
 
-        raw_wavefront = revert_to_data_input_format(
-            self.fft.data_format, raw_wavefront)
+        # raw_wavefront = revert_to_data_input_format(
+        #     self.fft.data_format, raw_wavefront)
         self.wavefront = raw_wavefront
 
         return raw_wavefront
@@ -288,8 +293,8 @@ def find_peaks_qlsi(ft_data, periodicity=4, copy=True):
     ft_data[:, cy - 3:cy + 3] = 0
 
     # circular bandpass according to periodicity
-    fx = np.fft.fftshift(np.fft.fftfreq(ft_data.shape[0])).reshape(-1, 1)
-    fy = np.fft.fftshift(np.fft.fftfreq(ft_data.shape[1])).reshape(1, -1)
+    fx = np.fft.fftshift(np.fft.fftfreq(ft_data.shape[-2])).reshape(-1, 1)
+    fy = np.fft.fftshift(np.fft.fftfreq(ft_data.shape[-1])).reshape(1, -1)
     frmask1 = np.sqrt(fx ** 2 + fy ** 2) > 1 / (periodicity * .8)
     frmask2 = np.sqrt(fx ** 2 + fy ** 2) < 1 / (periodicity * 1.2)
     ft_data[np.logical_or(frmask1, frmask2)] = 0
@@ -302,10 +307,11 @@ def find_peaks_qlsi(ft_data, periodicity=4, copy=True):
     return fx[i1x, 0], fy[0, i1y]
 
 
-def rotate_noreshape(arr, angle, mode="mirror", reshape=False):
+def rotate_noreshape(arr, angle, axes, mode="mirror", reshape=False):
     return scipy.ndimage.rotate(
         arr,  # input
         angle=np.rad2deg(angle),  # angle
+        axes=axes,
         reshape=reshape,  # reshape
         order=0,  # order
         mode=mode,  # mode

tests/test_fourier_base.py

@@ -152,6 +152,7 @@ def test_scale_to_filter_qlsi():
     ifr_phase = ifr.phase[0]
 
     phase = unwrap_phase(ifh_phase - ifr_phase)
+
     assert phase.shape == (720, 720)
     assert np.allclose(phase.mean(), 0.12434563269684816, atol=1e-6)
 

tests/test_qlsi.py

@@ -1,9 +1,11 @@
 import pathlib
+import pytest
 
 import h5py
 import numpy as np
-import qpretrieve
+from skimage.restoration import unwrap_phase
 
+import qpretrieve
 
 data_path = pathlib.Path(__file__).parent / "data"
 
@@ -29,3 +31,107 @@ def test_qlsi_phase():
         assert qlsi.phase.argmax() == 242294
         assert np.allclose(qlsi.phase.max(), 0.9343997734657971,
                            atol=0, rtol=1e-12)
+
+
+def test_qlsi_fftfreq_reshape_2d_3d(hologram):
+    data_2d = hologram
+    data_3d, _ = qpretrieve.data_input._convert_2d_to_3d(data_2d)
+
+    fx_2d = np.fft.fftfreq(data_2d.shape[-1]).reshape(-1, 1)
+    fx_3d = np.fft.fftfreq(data_3d.shape[-1]).reshape(data_3d.shape[0], -1, 1)
+
+    fy_2d = np.fft.fftfreq(data_2d.shape[-2]).reshape(1, -1)
+    fy_3d = np.fft.fftfreq(data_3d.shape[-2]).reshape(data_3d.shape[0], 1, -1)
+
+    assert np.array_equal(fx_2d, fx_3d[0])
+    assert np.array_equal(fy_2d, fy_3d[0])
+
+
+@pytest.mark.xfail
+def test_qlsi_unwrap_phase_2d_3d():
+    """
+    Check whether skimage unwrap_2d and unwrap_3d give the same result.
+    In other words, does unwrap_3d apply th unwrapping along the z axis.
+
+    Answer is no, they are different. unwrap_3d is designed for 3D data that
+    is to be unwrapped on all axes at once.
+
+    """
+    with h5py.File(data_path / "qlsi_paa_bead.h5") as h5:
+        image = h5["0"][:]
+
+    # Standard analysis pipeline
+    pipeline_kws = {
+        'wavelength': 550e-9,
+        'qlsi_pitch_term': 1.87711e-08,
+        'filter_name': "disk",
+        'filter_size': 180,
+        'filter_size_interpretation': "frequency index",
+        'scale_to_filter': False,
+        'invert_phase': False
+    }
+
+    data_2d = image
+    data_3d, _ = qpretrieve.data_input._convert_2d_to_3d(data_2d)
+
+    ft_2d = qpretrieve.fourier.FFTFilterNumpy(data_2d, subtract_mean=False)
+    ft_3d = qpretrieve.fourier.FFTFilterNumpy(data_3d, subtract_mean=False)
+
+    pipeline_kws["sideband_freq"] = qpretrieve.interfere. \
+        if_qlsi.find_peaks_qlsi(ft_2d.fft_origin[0])
+
+    hx_2d = ft_2d.filter(filter_name=pipeline_kws["filter_name"],
+                         filter_size=pipeline_kws["filter_size"],
+                         scale_to_filter=pipeline_kws["scale_to_filter"],
+                         freq_pos=pipeline_kws["sideband_freq"])
+    hx_3d = ft_3d.filter(filter_name=pipeline_kws["filter_name"],
+                         filter_size=pipeline_kws["filter_size"],
+                         scale_to_filter=pipeline_kws["scale_to_filter"],
+                         freq_pos=pipeline_kws["sideband_freq"])
+
+    assert np.array_equal(hx_2d, hx_3d)
+
+    px_3d_loop = np.zeros_like(hx_3d)
+    for i, _hx in enumerate(hx_3d):
+        px_3d_loop[i] = unwrap_phase(np.angle(_hx))
+
+    px_2d = unwrap_phase(np.angle(hx_2d[0]))
+    px_3d = unwrap_phase(np.angle(hx_3d))
+
+    assert np.array_equal(px_2d, px_3d_loop[0])  # this passes
+    assert np.array_equal(px_2d, px_3d[0])  # this fails
+
+
+def test_qlsi_rotate_2d_3d(hologram):
+    data_2d = hologram
+    data_3d, _ = qpretrieve.data_input._convert_2d_to_3d(data_2d)
+
+    rot_2d = qpretrieve.interfere.if_qlsi.rotate_noreshape(
+        data_2d,
+        angle=2,
+        axes=(1, 0),  # this was the default used before
+        reshape=False,
+    )
+    rot_3d = qpretrieve.interfere.if_qlsi.rotate_noreshape(
+        data_3d,
+        angle=2,
+        axes=(-1, -2),  # the y and x axes
+        reshape=False,
+    )
+
+    assert np.array_equal(rot_2d, rot_3d[0])
+
+
+def test_qlsi_pad_2d_3d(hologram):
+    data_2d = hologram
+    data_3d, _ = qpretrieve.data_input._convert_2d_to_3d(data_2d)
+
+    sx, sy = data_2d.shape[-2:]
+    gradpad_2d = np.pad(
+        data_2d, ((sx // 2, sx // 2), (sy // 2, sy // 2)),
+        mode="constant", constant_values=0)
+    gradpad_3d = np.pad(
+        data_3d, ((0, 0), (sx // 2, sx // 2), (sy // 2, sy // 2)),
+        mode="constant", constant_values=0)
+
+    assert np.array_equal(gradpad_2d, gradpad_3d[0])