savu paganin work6

Author: dkazanc

httomolibgpu/prep/phase.py

@@ -97,20 +97,15 @@ def paganin_filter(
     fft_tomo = fft2(padded_tomo, axes=(-2, -1), overwrite_x=True)
 
     # calculate alpha constant
-    alpha = _calculate_alpha(energy, distance, ratio_delta_beta)
+    alpha = _calculate_alpha(energy, distance / 1e-6, ratio_delta_beta)
 
     # Compute the reciprocal grid
-    # pixel_size *= 1e6  # rescaling to meters
     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)) / (4 * math.pi))
-        )
+        1.0 / (1.0 + alpha * (cp.add.outer(cp.square(indx), cp.square(indy))))
     )
 
     phase_filter = phase_filter / phase_filter.max()  # normalisation
@@ -142,18 +137,10 @@ def paganin_filter(
     return _log_kernel(tomo)
 
 
-def _calculate_alpha(energy, distance, ratio_delta_beta):
-    return _wavelength(energy) * distance * ratio_delta_beta
-
-
-# the scaling is different here and doesn't follow the original formula
-def _calculate_alpha_savu(energy, distance_micron, ratio_delta_beta):
+def _calculate_alpha(energy, distance_micron, ratio_delta_beta):
     return (
-        ((1240.0 / energy) * 10.0 ** (-9))
-        * distance_micron
-        * ratio_delta_beta
-        * math.pi
-    )
+        _wavelength_micron(energy) * distance_micron / (4 * math.pi)
+    ) * ratio_delta_beta
 
 
 def _shift_bit_length(x: int) -> int:
@@ -216,9 +203,10 @@ def _pad_projections_to_second_power(
     return padded_tomo, tuple(pad_list)
 
 
-def _wavelength(energy: float) -> float:
-    # for photons: E = 1keV -> 1.23984193 nm
-    return (1.23984193e-6) / energy
+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_coord(pixel_size: float, num_grid: int) -> cp.ndarray:
@@ -254,7 +242,7 @@ def paganin_filter_savu_legacy(
     """
     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. The ratio_delta_beta parameter here follows implementation in Savu software.
-    The module will be retired in future in favour of paganin_filter. One can rescale parameter ratio_delta_beta / 4*2*pi to achieve the same effect in paganin_filter.
+    The module will be retired in future in favour of paganin_filter. One can rescale parameter ratio_delta_beta / 4 to achieve the same effect in paganin_filter.
 
     Parameters
     ----------
@@ -274,64 +262,5 @@ def paganin_filter_savu_legacy(
     cp.ndarray
         The 3D array of Paganin phase-filtered projection images.
     """
-    # Check the input data is valid
-    if tomo.ndim != 3:
-        raise ValueError(
-            f"Invalid number of dimensions in data: {tomo.ndim},"
-            " please provide a stack of 2D projections."
-        )
-
-    dz_orig, dy_orig, dx_orig = tomo.shape
-
-    # Perform padding to the power of 2 as FFT is O(n*log(n)) complexity
-    # TODO: adding other options of padding?
-    padded_tomo, pad_tup = _pad_projections_to_second_power(tomo)
-
-    dz, dy, dx = padded_tomo.shape
-
-    # 3D FFT of tomo data
-    padded_tomo = cp.asarray(padded_tomo, dtype=cp.complex64)
-    fft_tomo = fft2(padded_tomo, axes=(-2, -1), overwrite_x=True)
-
-    # calculate alpha constant as Savu does
-    micron = 10 ** (-6)
-    KeV = 1000.0
-    alpha = _calculate_alpha_savu(energy * KeV, distance * micron, ratio_delta_beta)
-
-    # Compute the reciprocal grid
-    # NOTE: pixel_size is ALREADY in microns, why it rescaled to micron again in Savu?
-    indx = _reciprocal_coord(pixel_size * micron, dy)
-    indy = _reciprocal_coord(pixel_size * micron, dx)
-
-    # Build Lorentzian-type filter
-    phase_filter = fftshift(
-        1.0 / (1.0 + alpha * (cp.add.outer(cp.square(indx), cp.square(indy))))
-    )
-
-    phase_filter = phase_filter / phase_filter.max()  # normalisation
 
-    # Filter projections
-    fft_tomo *= phase_filter
-
-    # Apply filter and take inverse FFT
-    ifft_filtered_tomo = ifft2(fft_tomo, axes=(-2, -1), overwrite_x=True).real
-
-    # slicing indices for cropping
-    slc_indices = (
-        slice(pad_tup[0][0], pad_tup[0][0] + dz_orig, 1),
-        slice(pad_tup[1][0], pad_tup[1][0] + dy_orig, 1),
-        slice(pad_tup[2][0], pad_tup[2][0] + dx_orig, 1),
-    )
-
-    # crop the padded filtered data:
-    tomo = ifft_filtered_tomo[slc_indices].astype(cp.float32)
-
-    # taking the negative log
-    _log_kernel = cp.ElementwiseKernel(
-        "C tomo",
-        "C out",
-        "out = -log(tomo)",
-        name="log_kernel",
-    )
-
-    return _log_kernel(tomo)
+    return paganin_filter(tomo, pixel_size, distance, energy, ratio_delta_beta / 4)