Add 'physical radius' as disk filter option

CHANGELOG

@@ -1,4 +1,5 @@
 0.6.2
+ - enh: add 'physical radius' as disk filter option (#3, #26)
  - ref: use get_available_interfaces in get_best_interface (#14, #25)
 0.6.1
  - ref: reorder best fft interface (#23, #24)

qpretrieve/interfere/base.py

@@ -1,3 +1,4 @@
+from __future__ import annotations
 import warnings
 from abc import ABC, abstractmethod
 from typing import Type
@@ -175,7 +176,10 @@ def compute_filter_size(
             self,
             filter_size: float,
             filter_size_interpretation: str,
-            sideband_freq: tuple[float, float] = None) -> float:
+            sideband_freq: tuple[float, float] = None,
+            pixel_size: float | None = None,
+            numerical_aperture: float | None = None,
+            wavelength: float | None = None) -> float:
         """Compute the actual filter size in Fourier space"""
         if filter_size_interpretation == "frequency":
             # convert frequency to frequency index
@@ -201,6 +205,28 @@ def compute_filter_size(
                                  + f"'{filter_size}'!")
             # convert to frequencies (compatible with fx and fy)
             fsize = filter_size / self.fft.shape[-2]
+        elif filter_size_interpretation == "physical radius":
+            if not all(v is not None and v > 0 for v in (
+                    pixel_size, numerical_aperture, wavelength)):
+                raise ValueError(
+                    "For `filter_size_interpretation='physical radius'`, "
+                    "`pixel_size`, `numerical_aperture`, and `wavelength` "
+                    "must be set and must be positive.")
+            if filter_size <= 0:
+                raise ValueError(
+                    "For `filter_size_interpretation='physical radius'`, "
+                    "`filter_size` must be positive.")
+            # base radius in Fourier pixels: # r ~= n * dx * NA / lambda
+            # use detector-space size n from unpadded input stack.
+            n = float(max(self.fft.origin.shape[-2:]))
+            radius_px = n * pixel_size * numerical_aperture / wavelength
+            # `filter_size` acts as scaling factor
+            radius_px *= filter_size
+            if radius_px >= self.fft.shape[-2] / 2:
+                raise ValueError(
+                    "Physical-radius-derived filter size exceeds Fourier "
+                    f"limit {self.fft.shape[-2] / 2}; got '{radius_px}'.")
+            fsize = radius_px / self.fft.shape[-2]
         else:
             raise ValueError("Invalid value for `filter_size_interpretation`: "
                              + f"'{filter_size_interpretation}'")
@@ -221,7 +247,7 @@ def get_pipeline_kw(self, key):
 
     @abstractmethod
     def run_pipeline(self, **pipeline_kws):
-        """Perform pipeline analysis, populating `self.field`
+        r"""Perform pipeline analysis, populating `self.field`
 
         Parameters
         ----------
@@ -237,6 +263,12 @@ def run_pipeline(self, **pipeline_kws):
             (this is the default). If set to "frequency index", the filter
             size is interpreted as a Fourier frequency index ("pixel size")
             and must be between 0 and `max(hologram.shape)/2`.
+            If set to "physical radius", the radius is derived from
+            :math:`r \approx n dx NA / \lambda` (in Fourier pixels), where
+            `n` is the input image size in pixels, `dx` is `pixel_size`,
+            `NA` is `numerical_aperture`, and :math:`\lambda` is
+            `wavelength`. In this mode, `filter_size` is a scaling factor
+            (use `filter_size=1.0` for the direct formula).
         scale_to_filter: bool or float
             Crop the image in Fourier space after applying the filter,
             effectively removing surplus (zero-padding) data and
@@ -252,6 +284,15 @@ def run_pipeline(self, **pipeline_kws):
         sideband_freq: tuple of floats
             Frequency coordinates of the sideband to use. By default,
             a heuristic search for the sideband is done.
+        pixel_size: float
+            Sensor pixel size `dx` in meters, used when
+            `filter_size_interpretation="physical radius"`.
+        numerical_aperture: float
+            Collection NA, used when
+            `filter_size_interpretation="physical radius"`.
+        wavelength: float
+            Illumination wavelength in meters, used when
+            `filter_size_interpretation="physical radius"`.
         invert_phase: bool
             Invert the phase data.
         """

qpretrieve/interfere/if_oah.py

@@ -49,6 +49,10 @@ def run_pipeline(self, **pipeline_kws) -> xp.ndarray:
             (this is the default). If set to "frequency index", the filter
             size is interpreted as a Fourier frequency index ("pixel size")
             and must be between 0 and `max(hologram.shape)/2`.
+            If set to "physical radius", the base radius is
+            :math:`r \approx n dx NA / \lambda` in Fourier pixels. In this
+            mode `filter_size` is a scaling factor (`1.0` means direct
+            physical radius).
         scale_to_filter: bool or float
             Crop the image in Fourier space after applying the filter,
             effectively removing surplus (zero-padding) data and
@@ -66,6 +70,12 @@ def run_pipeline(self, **pipeline_kws) -> xp.ndarray:
             a heuristic search for the sideband is done.
             If you pass a 3D array, the first hologram is used to
             determine the sideband frequencies.
+        pixel_size: float
+            Sensor pixel size `dx` in meters for physical-radius mode.
+        numerical_aperture: float
+            Collection NA for physical-radius mode.
+        wavelength: float
+            Illumination wavelength in meters for physical-radius mode.
         invert_phase: bool
             Invert the phase data.
         """
@@ -82,7 +92,10 @@ def run_pipeline(self, **pipeline_kws) -> xp.ndarray:
             filter_size=pipeline_kws["filter_size"],
             filter_size_interpretation=(
                 pipeline_kws["filter_size_interpretation"]),
-            sideband_freq=pipeline_kws["sideband_freq"])
+            sideband_freq=pipeline_kws["sideband_freq"],
+            pixel_size=pipeline_kws.get("pixel_size"),
+            numerical_aperture=pipeline_kws.get("numerical_aperture"),
+            wavelength=pipeline_kws.get("wavelength"))
 
         # perform filtering
         filter_size = float(fsize)

qpretrieve/interfere/if_qlsi.py

@@ -73,6 +73,10 @@ def run_pipeline(self, **pipeline_kws) -> xp.ndarray:
             (this is the default). If set to "frequency index", the filter
             size is interpreted as a Fourier frequency index ("pixel size")
             and must be between 0 and `max(hologram.shape)/2`.
+            If set to "physical radius", the base radius is
+            :math:`r \approx n dx NA / \lambda` in Fourier pixels. In this
+            mode `filter_size` is a scaling factor (`1.0` means direct
+            physical radius).
         scale_to_filter: bool or float
             Crop the image in Fourier space after applying the filter,
             effectively removing surplus (zero-padding) data and
@@ -90,6 +94,10 @@ def run_pipeline(self, **pipeline_kws) -> xp.ndarray:
             a heuristic search for the sideband is done.
             If you pass a 3D array, the first hologram is used to
             determine the sideband frequencies.
+        pixel_size: float
+            Sensor pixel size `dx` in meters for physical-radius mode.
+        numerical_aperture: float
+            Collection NA for physical-radius mode.
         invert_phase: bool
             Invert the phase data.
         wavelength: float
@@ -129,7 +137,10 @@ def run_pipeline(self, **pipeline_kws) -> xp.ndarray:
             filter_size=pipeline_kws["filter_size"],
             filter_size_interpretation=(
                 pipeline_kws["filter_size_interpretation"]),
-            sideband_freq=pipeline_kws["sideband_freq"])
+            sideband_freq=pipeline_kws["sideband_freq"],
+            pixel_size=pipeline_kws.get("pixel_size"),
+            numerical_aperture=pipeline_kws.get("numerical_aperture"),
+            wavelength=pipeline_kws.get("wavelength"))
 
         # get pitch ratio
         qlsi_pitch_term = pipeline_kws["qlsi_pitch_term"]

tests/test_interfere_base.py

@@ -76,7 +76,10 @@ def test_interfere_base_get_data_with_input_layout_fft_warning():
     assert holo.field.shape == (1, 200, 210)
 
     fft_orig1 = holo.get_data_with_input_layout(data="fft_filtered")
-    fft_orig2 = holo.get_data_with_input_layout(data="fft")
+    with pytest.warns(
+            UserWarning, match="You have asked for 'fft' which is a class. "
+                               "Returning 'fft_filtered'. "):
+        fft_orig2 = holo.get_data_with_input_layout(data="fft")
     assert fft_orig1.shape == fft_orig2.shape
 
 

tests/test_oah.py

@@ -64,6 +64,68 @@ def test_get_field_error_invalid_interpretation(hologram):
         holo.run_pipeline(filter_size_interpretation="blequency")
 
 
+def test_get_field_interpretation_physical_radius_matches_frequency_index(
+        hologram):
+    """Verify 'physical radius' with matching 'frequency index' input"""
+    data = hologram
+    holo = qpretrieve.OffAxisHologram(data)
+
+    pixel_size = 6.5e-6
+    numerical_aperture = 0.03
+    wavelength = 532e-9
+
+    n = float(max(holo.fft.origin.shape[-2:]))
+    radius_px = n * pixel_size * numerical_aperture / wavelength
+
+    res_phys = holo.run_pipeline(
+        filter_name="disk",
+        filter_size=1.0,
+        filter_size_interpretation="physical radius",
+        pixel_size=pixel_size,
+        numerical_aperture=numerical_aperture,
+        wavelength=wavelength,
+    )
+    res_freq_idx = holo.run_pipeline(
+        filter_name="disk",
+        filter_size=radius_px,
+        filter_size_interpretation="frequency index",
+    )
+    assert np.all(res_phys == res_freq_idx)
+
+
+def test_get_field_interpretation_physical_radius_requires_parameters(
+        hologram):
+    data = hologram
+    holo = qpretrieve.OffAxisHologram(data)
+    with pytest.raises(ValueError, match="must be set and must be positive"):
+        holo.run_pipeline(
+            filter_size_interpretation="physical radius",
+            filter_size=1.0,
+        )
+
+    pixel_size = -6.5e-6  # negative will create error
+    numerical_aperture = 0.03
+    wavelength = 532e-9
+    with pytest.raises(ValueError, match="must be set and must be positive"):
+        holo.run_pipeline(
+            filter_size_interpretation="physical radius",
+            pixel_size=pixel_size,
+            numerical_aperture=numerical_aperture,
+            wavelength=wavelength,
+            filter_size=1.0,
+        )
+
+    pixel_size = 6.5e-6
+    with pytest.raises(ValueError, match="`filter_size` must be positive"):
+        holo.run_pipeline(
+            filter_size_interpretation="physical radius",
+            pixel_size=pixel_size,
+            numerical_aperture=numerical_aperture,
+            wavelength=wavelength,
+            filter_size=-1.0,
+        )
+
+
 def test_get_field_filter_names(hologram):
     data_2d = hologram
     data_3d, _ = convert_data_to_3d_array_layout(data_2d)