added documentation and non regression test

nabu/preproc/phase.py

@@ -51,8 +51,12 @@ class PhaseRetrieval(object):
         use_R2C=True,
     ):
         """
-        Paganin Phase Retrieval for an infinitely distant point source.
-        Formula (10) in [1].
+        Base class for Paganin Phase Retrieval for an infinitely distant point source.
+        This class must not be used as it is, but through a derived class which implements
+        the compute_filter method.
+        Derived class that you can use are
+        PaganinPhaseRetrieval which implements   Formula (10) in [1].
+        Paganin2020PhaseRetrieval which implements   Formula (17) in [2].
 
         Parameters
         ----------
@@ -134,32 +138,22 @@ class PhaseRetrieval(object):
              to its edges, including the edges. See ``numpy.pad()``.
 
 
-        **Formulas**
-
-        The radio is divided, in the Fourier domain, by the original "Paganin filter" `[1]`.
-
-        .. math::
-
-          F = 1 + \\frac{\\delta}{\\beta} \\lambda D \\pi |k|^2
-
-        where k is the wave vector.
-
 
         References
         -----------
         [1] D. Paganin Et Al, "Simultaneous phase and amplitude extraction
             from a single defocused image of a homogeneous object",
             Journal of Microscopy, Vol 206, Part 1, 2002
+        [2] D. Paganin et Al. "Boosting spatial resolution by incorporating
+            periodic boundary conditions into single-distance hard-x-ray
+            phase retrieval ",
+            Journal of Optics, Volume 22, Number 11, 2020
         """
-        self._init_parameters(
-            distance, energy, pixel_size, delta_beta, padding, use_R2C
-        )
+        self._init_parameters(distance, energy, pixel_size, delta_beta, padding, use_R2C)
         self._calc_shape(shape, margin)
         self.compute_filter()
 
-    def _init_parameters(
-        self, distance, energy, pixel_size, delta_beta, padding, use_R2C
-    ):
+    def _init_parameters(self, distance, energy, pixel_size, delta_beta, padding, use_R2C):
         self.distance_cm = distance
         self.distance_micron = distance * 1e4
         self.energy_kev = energy
@@ -254,7 +248,6 @@ class PhaseRetrieval(object):
             return n
         return self.powers[idx]
 
-
     def pad_with_values(self, data, top_val=0, bottom_val=0, left_val=0, right_val=0):
         """
         Pad the data into `self.padded_data` with values.
@@ -284,9 +277,7 @@ class PhaseRetrieval(object):
         self.data_padded = np.roll(self.data_padded, (-Pu, -Pl), axis=(0, 1))
 
     def _pad_zeros(self, data):
-        return self.pad_with_values(
-            data, top_val=0, bottom_val=0, left_val=0, right_val=0
-        )
+        return self.pad_with_values(data, top_val=0, bottom_val=0, left_val=0, right_val=0)
 
     def _pad_mean(self, data):
         """
@@ -332,8 +323,7 @@ class PhaseRetrieval(object):
         padding_method = PaddingMode.from_value(padding_method)
         if padding_method not in self.padding_methods:
             raise ValueError(
-                "Unknown padding method %s. Available are: %s"
-                % (padding_method, str(list(self.padding_methods.keys())))
+                "Unknown padding method %s. Available are: %s" % (padding_method, str(list(self.padding_methods.keys())))
             )
         pad_func = self.padding_methods[padding_method]
         pad_func(data)
@@ -360,11 +350,30 @@ class PhaseRetrieval(object):
         """
         return lmicron_to_db(Lmicron, self.energy_kev, self.distance_micron)
 
+    def compute_filter(self):
+        raise NotImplementedError("This class must be used as a derived class, which impelemnts the filter formula ")
+
     __call__ = apply_filter
 
 
 class PaganinPhaseRetrieval(PhaseRetrieval):
-    """The original Paganin method
+    """The original Paganin method it redefine the compute_filter method
+
+    **Formulas**
+
+    The radio is divided, in the Fourier domain, by the original "Paganin filter" `[1]`, formula (10).
+
+    .. math::
+
+      F = 1 + \\frac{\\delta}{\\beta} \\lambda D \\pi |k|^2
+
+    where k is the wave vector.
+
+    References
+    -----------
+    [1] D. Paganin Et Al, "Simultaneous phase and amplitude extraction
+        from a single defocused image of a homogeneous object",
+        Journal of Microscopy, Vol 206, Part 1, 2002
     """
 
     def compute_filter(self):
@@ -386,8 +395,15 @@ class PaganinPhaseRetrieval(PhaseRetrieval):
 
 
 class Paganin2020PhaseRetrieval(PhaseRetrieval):
-    """The 2020 Paganin method
-    TODO: add reference to the paper
+    """The 2020 Paganin method, it redefine the compute_filter method according
+    to formula (17) in [1]
+
+    References
+    -----------
+    [1] D. Paganin et Al. "Boosting spatial resolution by incorporating
+        periodic boundary conditions into single-distance hard-x-ray
+        phase retrieval ",
+        Journal of Optics, Volume 22, Number 11, 2020
     """
 
     def compute_filter(self):
@@ -401,12 +417,7 @@ class Paganin2020PhaseRetrieval(PhaseRetrieval):
 
         self._coords_grid = np.add.outer(fy ** 2, fx ** 2)
 
-        k2_modified = (
-            (2.0 - np.add.outer(np.cos(2 * np.pi * fy), np.cos(2 * np.pi * fx)))
-            * 2.0
-            / (2 * np.pi)
-            / (2 * np.pi)
-        )
+        k2_modified = (2.0 - np.add.outer(np.cos(2 * np.pi * fy), np.cos(2 * np.pi * fx))) * 2.0 / (2 * np.pi) / (2 * np.pi)
 
         k2 = self._coords_grid
         D = self.distance_micron

nabu/preproc/tests/test_paganin.py

 import pytest
 import numpy as np
-from nabu.preproc.phase import PaganinPhaseRetrieval, PaddingMode
+from nabu.preproc.phase import PaganinPhaseRetrieval, Paganin2020PhaseRetrieval, PaddingMode
 from nabu.testutils import get_data
 from nabu.thirdparty.tomopy_phase import retrieve_phase
 from nabu.cuda.utils import get_cuda_context, __has_pycuda__, __has_cufft__
+
 if __has_pycuda__:
     import pycuda.gpuarray as garray
-    from nabu.preproc.phase_cuda import CudaPaganinPhaseRetrieval
+    from nabu.preproc.phase_cuda import CudaPaganinPhaseRetrieval, CudaPaganin2020PhaseRetrieval
 
 scenarios = [
     {
@@ -17,7 +18,8 @@ scenarios = [
     }
 ]
 
-@pytest.fixture(scope='class', params=scenarios)
+
+@pytest.fixture(scope="class", params=scenarios)
 def bootstrap(request):
     cls = request.cls
     cls.paganin_config = request.param
@@ -28,7 +30,18 @@ def bootstrap(request):
     cls.paganin = PaganinPhaseRetrieval(cls.data.shape, **cls.paganin_config)
 
 
-@pytest.mark.usefixtures('bootstrap')
+@pytest.fixture(scope="class", params=scenarios)
+def bootstrap2020(request):
+    cls = request.cls
+    cls.paganin_config = request.param
+
+    cls.data = get_data("mri_proj_astra.npz")["data"]
+    cls.rtol = 3e-3
+    cls.rtol_pag = 5e-3
+    cls.paganin = Paganin2020PhaseRetrieval(cls.data.shape, **cls.paganin_config)
+
+
+@pytest.mark.usefixtures("bootstrap")
 class TestPaganin:
     """
     Test the Paganin phase retrieval.
@@ -38,34 +51,67 @@ class TestPaganin:
     def crop_to_margin(self, data):
         s0, s1 = self.paganin.shape_inner
         ((U, _), (L, _)) = self.paganin.margin
-        return data[U:U+s0, L:L+s1]
+        return data[U : U + s0, L : L + s1]
 
     def test_paganin(self):
         data_tomopy = np.atleast_3d(np.copy(self.data)).T
 
         res_tomopy = retrieve_phase(
             data_tomopy,
-            pixel_size=self.paganin.pixel_size_micron*1e-4,
+            pixel_size=self.paganin.pixel_size_micron * 1e-4,
             dist=self.paganin.distance_cm,
             energy=self.paganin.energy_kev,
-            alpha=1./(4*3.141592**2*self.paganin.delta_beta),
+            alpha=1.0 / (4 * 3.141592 ** 2 * self.paganin.delta_beta),
         )
         res_tomopy = self.crop_to_margin(res_tomopy[0].T)
 
         res = self.paganin.apply_filter(self.data)
 
-        errmax = np.max(np.abs(res - res_tomopy)/np.max(res_tomopy))
+        errmax = np.max(np.abs(res - res_tomopy) / np.max(res_tomopy))
         assert errmax < self.rtol_pag, "Max error is too high"
 
-
-    @pytest.mark.skipif(not(__has_pycuda__ and __has_cufft__), reason="Need pycuda and scikit-cuda for this test")
+    @pytest.mark.skipif(not (__has_pycuda__ and __has_cufft__), reason="Need pycuda and scikit-cuda for this test")
     def test_gpu_paganin(self):
-        gpu_paganin = CudaPaganinPhaseRetrieval(
-            self.data.shape,
-            **self.paganin_config
-        )
+        gpu_paganin = CudaPaganinPhaseRetrieval(self.data.shape, **self.paganin_config)
         ref = self.paganin.apply_filter(self.data)
         res = gpu_paganin.apply_filter(self.data)
-        errmax = np.max(np.abs((res - ref)/np.max(ref)))
+        errmax = np.max(np.abs((res - ref) / np.max(ref)))
         assert errmax < self.rtol, "Max error is too high"
 
+
+@pytest.mark.usefixtures("bootstrap2020")
+class TestPaganin2020:
+    """
+    Test the Paganin2020 phase retrieval.
+    The reference implementation is tomopy.
+    """
+
+    def crop_to_margin(self, data):
+        s0, s1 = self.paganin.shape_inner
+        ((U, _), (L, _)) = self.paganin.margin
+        return data[U : U + s0, L : L + s1]
+
+    def test_paganin(self):
+        data_tomopy = np.atleast_3d(np.copy(self.data)).T
+
+        res_tomopy = retrieve_phase(
+            data_tomopy,
+            pixel_size=self.paganin.pixel_size_micron * 1e-4,
+            dist=self.paganin.distance_cm,
+            energy=self.paganin.energy_kev,
+            alpha=1.0 / (4 * 3.141592 ** 2 * self.paganin.delta_beta),
+        )
+        res_tomopy = self.crop_to_margin(res_tomopy[0].T)
+
+        res = self.paganin.apply_filter(self.data)
+
+        errmax = np.max(np.abs(res - res_tomopy) / np.max(res_tomopy))
+        assert errmax < self.rtol_pag, "Max error is too high"
+
+    @pytest.mark.skipif(not (__has_pycuda__ and __has_cufft__), reason="Need pycuda and scikit-cuda for this test")
+    def test_gpu_paganin(self):
+        gpu_paganin = CudaPaganinPhaseRetrieval(self.data.shape, **self.paganin_config)
+        ref = self.paganin.apply_filter(self.data)
+        res = gpu_paganin.apply_filter(self.data)
+        errmax = np.max(np.abs((res - ref) / np.max(ref)))
+        assert errmax < self.rtol, "Max error is too high"