From d7457e16dbe80fcd68dd84c767748db4e5d85858 Mon Sep 17 00:00:00 2001 From: Eoghan O'Connell Date: Thu, 21 Nov 2024 16:16:10 +0100 Subject: [PATCH] ref: align the 2d qlsi code to 3d --- qpretrieve/interfere/if_qlsi.py | 30 +++++---- tests/test_fourier_base.py | 1 + tests/test_qlsi.py | 108 +++++++++++++++++++++++++++++++- 3 files changed, 126 insertions(+), 13 deletions(-) diff --git a/qpretrieve/interfere/if_qlsi.py b/qpretrieve/interfere/if_qlsi.py index 8000b8b..917a8ad 100644 --- a/qpretrieve/interfere/if_qlsi.py +++ b/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 diff --git a/tests/test_fourier_base.py b/tests/test_fourier_base.py index 345d074..fbfe673 100644 --- a/tests/test_fourier_base.py +++ b/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) diff --git a/tests/test_qlsi.py b/tests/test_qlsi.py index f839114..d151fb9 100644 --- a/tests/test_qlsi.py +++ b/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])