Merge pull request #238 from DiamondLightSource/paganin_changes

Paganin filter changes

Author: dkazanc

docs/source/_static/figures/paganin/lprec_no_paganin.jpg

@@ -8,6 +8,16 @@ @article{Paganin02
   year={2002},
   publisher={Wiley Online Library}
 }
+@article{paganin2020boosting,
+  title={Boosting spatial resolution by incorporating periodic boundary conditions into single-distance hard-x-ray phase retrieval},
+  author={Paganin, David M and Favre-Nicolin, Vincent and Mirone, Alessandro and Rack, Alexander and Villanova, Julie and Olbinado, Margie P and Fernandez, Vincent and da Silva, Julio C and Pelliccia, Daniele},
+  journal={Journal of Optics},
+  volume={22},
+  number={11},
+  pages={115607},
+  year={2020},
+  publisher={IOP Publishing}
+}
 @article{vo2021data,
   title={Data processing methods and data acquisition for samples larger than the field of view in parallel-beam tomography},
   author={Vo, Nghia T and Atwood, Robert C and Drakopoulos, Michael and Connolley, Thomas},

docs/source/_static/figures/paganin/paganin_10.jpg

@@ -10,8 +10,8 @@ Here we present a list of methods of the HTTomolibGPU library with more detailed
 
    methods_list/correction_methods
    methods_list/stripe_removal_methods
+   methods_list/phase_contrast_methods
    methods_list/reconstruction_methods
    methods_list/denoising_methods
    methods_list/rescale_methods
-   
-   
+

docs/source/_static/figures/paganin/paganin_100.jpg

@@ -0,0 +1,110 @@
+.. _phase_contrast_module:
+
+Phase-contrast enhancement
+**************************
+
+Methods from the :mod:`httomolibgpu.prep.phase` module needed when the data requires phase contrast enhancement.
+
+**Description**
+
+In conventional X-ray tomography, image contrast comes mainly from absorption, i.e., how much a material attenuates X-rays.
+For some samples, especially biological, the absorption difference between structures can be insignificant and therefore the edge contrast
+can be lost. X-rays also undergo phase shifts when passing through materials with different refractive indices and
+these phase shifts carry valuable structural information that can be converted into contrast if properly reconstructed.
+
+The phase contrast enhances edges and interfaces, revealing fine internal structure that pure absorption tomography would miss.
+When we record propagation-based phase contrast images (e.g., at a certain distance from the sample), the detected intensity includes
+a mixture of absorption and phase effects. To reconstruct a clean 3D volume, we must separate phase from intensity.
+This process is called the phase retrieval. The Paganin method :cite:`Paganin02`  provides a simple, robust way to do this.
+
+The Paganin formula is the following, see eq.10 in :cite:`paganin2020boosting`:
+
+.. math:: T(x,y) = - \frac{1}{\mu}\ln\left (\mathcal{F}^{-1}\left
+    (\frac{\mathcal{F}\left [ I(x, y, z = \Delta) / I_{0} \right ]}{1 +
+        \frac{\delta\Delta}{\mu}\left ( k_x^2 + k_y^2 \right )}  \right )\right ),
+
+where:
+
+* :math:`I(x, y, z = \Delta)` - measured X-ray beam intensity at the propagation distance :math:`\Delta`.
+
+* :math:`I_{0}` - measured X-ray beam intensity at the zero-distance from the X-ray source (incident beam).
+
+* :math:`\Delta > 0` - propagation distance of the wavefront from sample to detector.
+
+* :math:`\mu` - linear attenuation coefficient of the single-material object defined as :math:`\mu = 2k\beta`, where :math:`k=\frac{2\pi}{\lambda}` is the wave-number corresponding to the vacuum wavelength :math:`\lambda`.
+
+* :math:`\delta` - the phase decrement, related to the phase shift of X-rays. It is the real part of the complex refractive index:  :math:`n = (1 - \delta) + i \beta`.
+
+* :math:`\beta` - the absorption index, related to the attenuation. It the complex part of the material refractive index: :math:`n = (1 - \delta) + i \beta`.
+
+One can re-write the formula above as:
+
+.. math:: T(x,y) = - \frac{1}{\mu}\ln\left (\mathcal{F}^{-1}\left
+    (\frac{\mathcal{F}\left [ I(x, y, z = \Delta) / I_{0} \right ]}{1 +
+        \alpha \left ( k_x^2 + k_y^2 \right )}  \right )\right ),
+
+where:
+
+.. math:: \alpha = \frac{\lambda \Delta \delta}{4 \pi \beta}.
+
+**Where and how to use it:**
+
+Works well for single-material and relatively homogeneous or weakly heterogeneous samples (e.g., biological tissues, polymers, soft matter). The filter is particularly useful for weakly absorbing samples imaged with hard X-rays, where traditional absorption contrast is poor.
+
+**What are the adjustable parameters:**
+
+* :code:`input data` is the flat/dark normalised raw data before the negative log. Note that the :math:`-ln()` operation is the part of the Paganin filter.
+
+* :code:`pixel_size` Detector pixel size (resolution) in MICRON units.
+
+* :code:`distance` Propagation distance (:math:`\Delta`) of the wavefront from sample to detector in METRE units.
+
+* :code:`energy` Incident beam energy in keV.
+
+* :code:`ratio_delta_beta` The ratio of :math:`\frac{\delta}{\beta}` is a critical parameter as it defines how the algorithm balances phase contrast versus absorption contrast. Higher ratio values lead to stronger smoothing, more phase contrast recovered, but potential loss of edge sharpness. It is recommended to keep :math:`\frac{\delta}{\beta} > 250` for weakly absorbing materials (like biological tissue, polymers, and light materials). Lower values :math:`\frac{\delta}{\beta} << 250` lead to less smoothing, more edges preserved.
+
+**Real data examples:**
+
+In this section we will be applying Paganin filter to data obtained at Diamond Light Source I12 beamline (beamtime NT41036-2, PI: G. Burca).
+The sample is a set of fixed bovine liver sections in a centrifuge tube. This is a biological sample and a good candidate for Paganin filter demonstration.
+
+We used :ref:`method_remove_all_stripe` to remove ring artifacts and :ref:`method_LPRec3d_tomobar` for reconstruction of the filtered projections.
+We also used the following parameters for Paganin filter, while varying only :code:`ratio_delta_beta`.
+
+.. code-block:: yaml
+
+    pixel_size: 32.4 # microns
+    distance: 3.2    # meters
+    energy: 53.0     # KeV
+
+.. list-table::
+
+
+    * - .. figure:: ../../_static/figures/paganin/lprec_no_paganin.jpg
+           :width: 335px
+
+           No Paganin filter applied. The reconstruction is too noisy and without any features.
+
+      - .. figure:: ../../_static/figures/paganin/paganin_10.jpg
+           :width: 335px
+
+           Paganin filter with :math:`\frac{\delta}{\beta} = 10`. Some features are recongnisable, but the image is still too noisy.
+
+
+.. list-table::
+
+
+    * - .. figure:: ../../_static/figures/paganin/paganin_100.jpg
+           :width: 200px
+
+           Paganin filter with :math:`\frac{\delta}{\beta} = 100`. Still a bit noisy.
+
+      - .. figure:: ../../_static/figures/paganin/paganin_500.jpg
+           :width: 200px
+
+           Paganin filter with :math:`\frac{\delta}{\beta} = 500`. Close to the optimal value.
+
+      - .. figure:: ../../_static/figures/paganin/paganin_2000.jpg
+           :width: 200px
+
+           Paganin filter with :math:`\frac{\delta}{\beta} = 2000`. Oversmoothed.

docs/source/_static/figures/paganin/paganin_2000.jpg

@@ -4,7 +4,7 @@
 from httomolibgpu.misc.rescale import rescale_to_int
 from httomolibgpu.prep.alignment import distortion_correction_proj_discorpy
 from httomolibgpu.prep.normalize import normalize
-from httomolibgpu.prep.phase import paganin_filter_tomopy
+from httomolibgpu.prep.phase import paganin_filter
 from httomolibgpu.prep.stripe import (
     remove_stripe_based_sorting,
     remove_stripe_ti,

docs/source/_static/figures/paganin/paganin_500.jpg

@@ -18,7 +18,7 @@
 # Created By  : Tomography Team at DLS <[email protected]>
 # Created Date: 01 November 2022
 # ---------------------------------------------------------------------------
-"""Modules for phase retrieval and phase-contrast enhancement"""
+"""Modules for phase retrieval and phase-contrast enhancement. For more detailed information, see :ref:`phase_contrast_module`."""
 
 import numpy as np
 from httomolibgpu import cupywrapper
@@ -42,35 +42,36 @@
 from httomolibgpu.misc.supp_func import data_checker
 
 __all__ = [
-    "paganin_filter_tomopy",
+    "paganin_filter",
 ]
 
 
-##-------------------------------------------------------------##
-# Adaptation of retrieve_phase (Paganin filter) from TomoPy
-def paganin_filter_tomopy(
+# This implementation originated from the TomoPy version. It has been modified to conform
+# different unit standards and also control of the filter driven by 'delta/beta' ratio
+# as opposed to 'alpha' in the TomoPy's implementation.
+def paganin_filter(
     tomo: cp.ndarray,
-    pixel_size: float = 1e-4,
-    dist: float = 50.0,
+    pixel_size: float = 1.28,
+    distance: float = 1.0,
     energy: float = 53.0,
-    alpha: float = 1e-3,
+    ratio_delta_beta: float = 250,
 ) -> cp.ndarray:
     """
-    Perform single-material phase retrieval from flats/darks corrected tomographic measurements. See
-    :cite:`Paganin02` for a reference.
+    Perform single-material phase retrieval from flats/darks corrected tomographic measurements. For more detailed information, see :ref:`phase_contrast_module`.
+    Also see :cite:`Paganin02` and :cite:`paganin2020boosting` for references.
 
     Parameters
     ----------
     tomo : cp.ndarray
         3D array of f/d corrected tomographic projections.
-    pixel_size : float, optional
-        Detector pixel size in cm.
-    dist : float, optional
-        Propagation distance of the wavefront in cm.
-    energy : float, optional
-        Energy of incident wave in keV.
-    alpha : float, optional
-        Regularization parameter, the ratio of delta/beta. Smaller values lead to less noise and more blur.
+    pixel_size : float
+        Detector pixel size (resolution) in micron units.
+    distance : float
+        Propagation distance of the wavefront from sample to detector in metre units.
+    energy : float
+        Beam energy in keV.
+    ratio_delta_beta : float
+        The ratio of delta/beta, where delta is the phase shift and real part of the complex material refractive index and beta is the absorption.
 
     Returns
     -------
@@ -84,7 +85,7 @@ def paganin_filter_tomopy(
             " please provide a stack of 2D projections."
         )
 
-    tomo = data_checker(tomo, verbosity=True, method_name="paganin_filter_tomopy")
+    tomo = data_checker(tomo, verbosity=True, method_name="paganin_filter")
 
     dz_orig, dy_orig, dx_orig = tomo.shape
 
@@ -98,11 +99,18 @@ def paganin_filter_tomopy(
     padded_tomo = cp.asarray(padded_tomo, dtype=cp.complex64)
     fft_tomo = fft2(padded_tomo, axes=(-2, -1), overwrite_x=True)
 
-    # Compute the reciprocal grid.
-    w2 = _reciprocal_grid(pixel_size, (dy, dx))
+    # calculate alpha constant
+    alpha = _calculate_alpha(energy, distance / 1e-6, ratio_delta_beta)
+
+    # Compute the reciprocal grid
+    indx = _reciprocal_coord(pixel_size, dy)
+    indy = _reciprocal_coord(pixel_size, dx)
+
+    # Build Lorentzian-type filter
+    phase_filter = fftshift(
+        1.0 / (1.0 + alpha * (cp.add.outer(cp.square(indx), cp.square(indy))))
+    )
 
-    # Build filter in the Fourier space.
-    phase_filter = fftshift(_paganin_filter_factor2(energy, dist, alpha, w2))
     phase_filter = phase_filter / phase_filter.max()  # normalisation
 
     # Filter projections
@@ -132,6 +140,12 @@ def paganin_filter_tomopy(
     return _log_kernel(tomo)
 
 
+def _calculate_alpha(energy, distance_micron, ratio_delta_beta):
+    return (
+        _wavelength_micron(energy) * distance_micron / (4 * math.pi)
+    ) * ratio_delta_beta
+
+
 def _shift_bit_length(x: int) -> int:
     return 1 << (x - 1).bit_length()
 
@@ -192,42 +206,12 @@ def _pad_projections_to_second_power(
     return padded_tomo, tuple(pad_list)
 
 
-def _paganin_filter_factor2(energy, dist, alpha, w2):
-    # Alpha represents the ratio of delta/beta.
-    return 1 / (_wavelength(energy) * dist * w2 / (4 * math.pi) + alpha)
-
-
-def _wavelength(energy: float) -> float:
-    SPEED_OF_LIGHT = 299792458e2  # [cm/s]
+def _wavelength_micron(energy: float) -> float:
+    SPEED_OF_LIGHT = 299792458e2 * 10000.0  # [microns/s]
     PLANCK_CONSTANT = 6.58211928e-19  # [keV*s]
     return 2 * math.pi * PLANCK_CONSTANT * SPEED_OF_LIGHT / energy
 
 
-def _reciprocal_grid(pixel_size: float, shape_proj: tuple) -> cp.ndarray:
-    """
-    Calculate reciprocal grid.
-
-    Parameters
-    ----------
-    pixel_size : float
-        Detector pixel size in cm.
-    shape_proj : tuple
-        Shape of the reciprocal grid along x and y axes.
-
-    Returns
-    -------
-    ndarray
-        Grid coordinates.
-    """
-    # Sampling in reciprocal space.
-    indx = _reciprocal_coord(pixel_size, shape_proj[0])
-    indy = _reciprocal_coord(pixel_size, shape_proj[1])
-    indx_sq = cp.square(indx)
-    indy_sq = cp.square(indy)
-
-    return cp.add.outer(indx_sq, indy_sq)
-
-
 def _reciprocal_coord(pixel_size: float, num_grid: int) -> cp.ndarray:
     """
     Calculate reciprocal grid coordinates for a given pixel size
@@ -236,7 +220,7 @@ def _reciprocal_coord(pixel_size: float, num_grid: int) -> cp.ndarray:
     Parameters
     ----------
     pixel_size : float
-        Detector pixel size in cm.
+        Detector pixel size in microns.
     num_grid : int
         Size of the reciprocal grid.