Merge pull request #215 from DiamondLightSource/remove_paganin_savu

retiring paganin savu method

Author: dkazanc

httomolibgpu/__init__.py

@@ -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_savu, paganin_filter_tomopy
+from httomolibgpu.prep.phase import paganin_filter_tomopy
 from httomolibgpu.prep.stripe import (
     remove_stripe_based_sorting,
     remove_stripe_ti,

httomolibgpu/cuda_kernels/paganin_filter_gen.cu

@@ -29,10 +29,8 @@
 from unittest.mock import Mock
 
 if cupy_run:
-    from httomolibgpu.cuda_kernels import load_cuda_module
     from cupyx.scipy.fft import fft2, ifft2, fftshift
 else:
-    load_cuda_module = Mock()
     fft2 = Mock()
     ifft2 = Mock()
     fftshift = Mock()
@@ -44,225 +42,10 @@
 from httomolibgpu.misc.supp_func import data_checker
 
 __all__ = [
-    "paganin_filter_savu",
     "paganin_filter_tomopy",
 ]
 
 
-## %%%%%%%%%%%%%%%%%%%%%%% paganin_filter %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  ##
-#: CuPy implementation of Paganin filter from Savu
-def paganin_filter_savu(
-    data: cp.ndarray,
-    ratio: float = 250.0,
-    energy: float = 53.0,
-    distance: float = 1.0,
-    resolution: float = 1.28,
-    pad_y: int = 100,
-    pad_x: int = 100,
-    pad_method: str = "edge",
-    increment: float = 0.0,
-) -> cp.ndarray:
-    """
-    Apply Paganin filter (for denoising or contrast enhancement) to
-    projections.
-
-    Parameters
-    ----------
-    data : cp.ndarray
-        The stack of projections to filter.
-
-    ratio : float, optional
-        Ratio of delta/beta.
-
-    energy : float, optional
-        Beam energy in keV.
-
-    distance : float, optional
-        Distance from sample to detector in metres.
-
-    resolution : float, optional
-        Pixel size in microns.
-
-    pad_y : int, optional
-        Pad the top and bottom of projections.
-
-    pad_x : int, optional
-        Pad the left and right of projections.
-
-    pad_method : str, optional
-        Numpy pad method to use.
-
-    increment : float, optional
-        Increment all values by this amount before taking the log.
-
-    Returns
-    -------
-    cp.ndarray
-        The stack of filtered projections.
-    """
-    # Check the input data is valid
-    if data.ndim != 3:
-        raise ValueError(
-            f"Invalid number of dimensions in data: {data.ndim},"
-            " please provide a stack of 2D projections."
-        )
-
-    data = data_checker(data, verbosity=True, method_name="paganin_filter_savu")
-
-    # Setup various values for the filter
-    _, height, width = data.shape
-    micron = 1e-6
-    keV = 1000.0
-    energy *= keV
-    resolution *= micron
-    wavelength = (1240.0 / energy) * 1e-9
-
-    height1 = height + 2 * pad_y
-    width1 = width + 2 * pad_x
-
-    # Define the paganin filter, taking into account the padding that will be
-    # applied to the projections (if any)
-
-    # Using raw kernel her as indexing is direct and it avoids a lot of temporaries
-    # and tiny kernels
-    module = load_cuda_module("paganin_filter_gen")
-    kernel = module.get_function("paganin_filter_gen")
-
-    # Apply padding to all the 2D projections
-    # Note: this takes considerable time on GPU...
-    data = cp.pad(data, ((0, 0), (pad_y, pad_y), (pad_x, pad_x)), mode=pad_method)
-
-    precond_kernel_float = cp.ElementwiseKernel(
-        "T data",
-        "T out",
-        """
-        if (isnan(data)) {
-            out = T(0);
-        } else if (isinf(data)) {
-            out = data < 0.0 ? -3.402823e38f : 3.402823e38f;  // FLT_MAX, not available in cupy
-        } else if (data == 0.0) {
-            out = 1.0;
-        } else {
-            out = data;
-        }
-        """,
-        name="paganin_precond_float",
-        no_return=True,
-    )
-    precond_kernel_int = cp.ElementwiseKernel(
-        "T data",
-        "T out",
-        """out = data == 0 ? 1 : data""",
-        name="paganin_precond_int",
-        no_return=True,
-    )
-
-    if data.dtype in (cp.float32, cp.float64):
-        precond_kernel_float(data, data)
-    else:
-        precond_kernel_int(data, data)
-
-    # avoid normalising in both directions - we include multiplier in the post_kernel
-    data = cp.asarray(data, dtype=cp.complex64)
-    data = fft2(data, axes=(-2, -1), overwrite_x=True, norm="backward")
-
-    # prepare filter here, while the GPU is busy with the FFT
-    filtercomplex = cp.empty((height1, width1), dtype=cp.complex64)
-    bx = 16
-    by = 8
-    gx = (width1 + bx - 1) // bx
-    gy = (height1 + by - 1) // by
-    kernel(
-        grid=(gx, gy, 1),
-        block=(bx, by, 1),
-        args=(
-            cp.int32(width1),
-            cp.int32(height1),
-            cp.float32(resolution),
-            cp.float32(wavelength),
-            cp.float32(distance),
-            cp.float32(ratio),
-            filtercomplex,
-        ),
-    )
-    data *= filtercomplex
-
-    data = ifft2(data, axes=(-2, -1), overwrite_x=True, norm="forward")
-
-    post_kernel = cp.ElementwiseKernel(
-        "C pci1, raw float32 increment, raw float32 ratio, raw float32 fft_scale",
-        "T out",
-        "out = -0.5 * ratio * log(abs(pci1) * fft_scale + increment)",
-        name="paganin_post_proc",
-        no_return=True,
-    )
-    fft_scale = 1.0 / (data.shape[1] * data.shape[2])
-    res = cp.empty((data.shape[0], height, width), dtype=cp.float32)
-    post_kernel(
-        data[:, pad_y : pad_y + height, pad_x : pad_x + width],
-        np.float32(increment),
-        np.float32(ratio),
-        np.float32(fft_scale),
-        res,
-    )
-    return res
-
-
-def _wavelength(energy: float) -> float:
-    SPEED_OF_LIGHT = 299792458e2  # [cm/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
-    and discretization.
-
-    Parameters
-    ----------
-    pixel_size : float
-        Detector pixel size in cm.
-    num_grid : int
-        Size of the reciprocal grid.
-
-    Returns
-    -------
-    ndarray
-        Grid coordinates.
-    """
-    n = num_grid - 1
-    rc = cp.arange(-n, num_grid, 2, dtype=cp.float32)
-    rc *= 2 * math.pi / (n * pixel_size)
-    return rc
-
-
-##-------------------------------------------------------------##
 ##-------------------------------------------------------------##
 # Adaptation of retrieve_phase (Paganin filter) from TomoPy
 def paganin_filter_tomopy(
@@ -412,3 +195,57 @@ def _pad_projections_to_second_power(
 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]
+    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
+    and discretization.
+
+    Parameters
+    ----------
+    pixel_size : float
+        Detector pixel size in cm.
+    num_grid : int
+        Size of the reciprocal grid.
+
+    Returns
+    -------
+    ndarray
+        Grid coordinates.
+    """
+    n = num_grid - 1
+    rc = cp.arange(-n, num_grid, 2, dtype=cp.float32)
+    rc *= 2 * math.pi / (n * pixel_size)
+    return rc

httomolibgpu/prep/phase.py

@@ -152,7 +152,7 @@ def FBP3d_tomobar(
     """
     Perform Filtered Backprojection (FBP) reconstruction using ASTRA toolbox :cite:`van2016fast` and
     ToMoBAR :cite:`kazantsev2020tomographic` wrappers.
-    This is a 3D recon from the CuPy array directly and using a custom built SINC filter for filtration in Fourier space, 
+    This is a 3D recon from the CuPy array directly and using a custom built SINC filter for filtration in Fourier space,
     see more in :ref:`method_FBP3d_tomobar`.
 
     Parameters
@@ -166,7 +166,7 @@ def FBP3d_tomobar(
     detector_pad : int
         Detector width padding with edge values to remove circle/arc type artifacts in the reconstruction.
     filter_freq_cutoff : float
-        Cutoff frequency parameter for the SINC filter, the lower values may produce better contrast but noisy reconstruction. The filter change will also affect the dynamic range of the reconstructed image. 
+        Cutoff frequency parameter for the SINC filter, the lower values may produce better contrast but noisy reconstruction. The filter change will also affect the dynamic range of the reconstructed image.
     recon_size : int, optional
         The [recon_size, recon_size] shape of the reconstructed slice in pixels.
         By default (None), the reconstructed size will be the dimension of the horizontal detector.