Merge branch 'main' into phase_cross_correlation-radway-62

Author: dkazanc

httomolibgpu/__init__.py

@@ -5,7 +5,7 @@
 from httomolibgpu.misc.morph import sino_360_to_180, data_resampler
 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.normalize import dark_flat_field_correction, minus_log
 from httomolibgpu.prep.phase import paganin_filter, paganin_filter_savu_legacy
 from httomolibgpu.prep.stripe import (
     remove_stripe_based_sorting,

httomolibgpu/prep/normalize.py

@@ -20,7 +20,6 @@
 # ---------------------------------------------------------------------------
 """Modules for raw projection data normalization"""
 
-import numpy as np
 from httomolibgpu import cupywrapper
 
 cp = cupywrapper.cp
@@ -34,25 +33,21 @@
     mean = Mock()
 
 from numpy import float32
-from typing import Tuple
 
-__all__ = ["normalize"]
 
+__all__ = ["dark_flat_field_correction", "minus_log"]
 
-def normalize(
+
+def dark_flat_field_correction(
     data: cp.ndarray,
     flats: cp.ndarray,
     darks: cp.ndarray,
     flats_multiplier: float = 1.0,
     darks_multiplier: float = 1.0,
     cutoff: float = 10.0,
-    minus_log: bool = True,
-    nonnegativity: bool = False,
-    remove_nans: bool = False,
 ) -> cp.ndarray:
     """
     Normalize raw projection data using the flat and dark field projections.
-    This is a raw CUDA kernel implementation with CuPy wrappers.
 
     Parameters
     ----------
@@ -68,17 +63,11 @@ def normalize(
         A multiplier to apply to darks, can work as an intensity compensation constant.
     cutoff : float
         Permitted maximum value for the normalised data.
-    minus_log : bool
-        Apply negative log to the normalised data.
-    nonnegativity : bool
-        Remove negative values in the normalised data.
-    remove_nans : bool
-        Remove NaN and Inf values in the normalised data.
 
-    Returns
+            Returns
     -------
     cp.ndarray
-        Normalised 3D tomographic data as a CuPy array.
+        Normalised by dark/flat fields 3D tomographic data as a CuPy array.
     """
     _check_valid_input_normalise(data, flats, darks)
 
@@ -99,16 +88,6 @@ def normalize(
         }
         float v = (float(data) - float(darks))/denom;
         """
-    if minus_log:
-        kernel += "v = -log(v);\n"
-        kernel_name += "_mlog"
-    if nonnegativity:
-        kernel += "if (v < 0.0f) v = 0.0f;\n"
-        kernel_name += "_nneg"
-    if remove_nans:
-        kernel += "if (isnan(v)) v = 0.0f;\n"
-        kernel += "if (isinf(v)) v = 0.0f;\n"
-        kernel_name += "_remnan"
     kernel += "if (v > cutoff) v = cutoff;\n"
     kernel += "if (v < -cutoff) v = cutoff;\n"
     kernel += "out = v;\n"
@@ -128,6 +107,24 @@ def normalize(
     return out
 
 
+def minus_log(data: cp.ndarray) -> cp.ndarray:
+    """
+    Apply -log(data) operation
+
+    Parameters
+    ----------
+    data : cp.ndarray
+        Data as a CuPy array.
+
+    Returns
+    -------
+    cp.ndarray
+        data after -log(data)
+    """
+
+    return -cp.log(data)
+
+
 def _check_valid_input_normalise(data, flats, darks) -> None:
     """Helper function to check the validity of inputs to normalisation functions"""
     if data.ndim != 3:

httomolibgpu/recon/algorithm.py

@@ -64,7 +64,6 @@ def FBP2d_astra(
     filter_d: Optional[float] = None,
     recon_size: Optional[int] = None,
     recon_mask_radius: float = 0.95,
-    neglog: bool = False,
     gpu_id: int = 0,
 ) -> np.ndarray:
     """
@@ -96,9 +95,6 @@ def FBP2d_astra(
         The radius of the circular mask that applies to the reconstructed slice in order to crop
         out some undesirable artifacts. The values outside the given diameter will be set to zero.
         To implement the cropping one can use the range [0.7-1.0] or set to 2.0 when no cropping required.
-    neglog: bool
-        Take negative logarithm on input data to convert to attenuation coefficient or a density of the scanned object. Defaults to False,
-        assuming that the negative log is taken either in normalisation procedure on with Paganin filter application.
     gpu_id : int
         A GPU device index to perform operation on.
 
@@ -120,8 +116,6 @@ def FBP2d_astra(
     reconstruction = np.empty(
         (recon_size, detY_size, recon_size), dtype=float32, order="C"
     )
-    _take_neg_log_np(data) if neglog else data
-
     # loop over detY slices
     for slice_index in range(0, detY_size):
         reconstruction[:, slice_index, :] = np.flipud(
@@ -145,7 +139,6 @@ def FBP3d_tomobar(
     filter_freq_cutoff: float = 0.35,
     recon_size: Optional[int] = None,
     recon_mask_radius: Optional[float] = 0.95,
-    neglog: bool = False,
     gpu_id: int = 0,
 ) -> cp.ndarray:
     """
@@ -174,9 +167,6 @@ def FBP3d_tomobar(
         The radius of the circular mask that applies to the reconstructed slice in order to crop
         out some undesirable artifacts. The values outside the given diameter will be set to zero.
         To implement the cropping one can use the range [0.7-1.0] or set to 2.0 when no cropping required.
-    neglog: bool
-        Take negative logarithm on input data to convert to attenuation coefficient or a density of the scanned object. Defaults to False,
-        assuming that the negative log is taken either in normalisation procedure on with Paganin filter application.
     gpu_id : int
         A GPU device index to perform operation on.
 
@@ -191,7 +181,7 @@ def FBP3d_tomobar(
     )
 
     reconstruction = RecToolsCP.FBP(
-        _take_neg_log(data) if neglog else data,
+        data,
         cutoff_freq=filter_freq_cutoff,
         recon_mask_radius=recon_mask_radius,
         data_axes_labels_order=input_data_axis_labels,
@@ -214,7 +204,6 @@ def LPRec3d_tomobar(
     power_of_2_cropping: Optional[bool] = False,
     min_mem_usage_filter: Optional[bool] = True,
     min_mem_usage_ifft2: Optional[bool] = True,
-    neglog: bool = False,
 ) -> cp.ndarray:
     """
     Fourier direct inversion in 3D on unequally spaced (also called as Log-Polar) grids using
@@ -243,9 +232,6 @@ def LPRec3d_tomobar(
         The radius of the circular mask that applies to the reconstructed slice in order to crop
         out some undesirable artifacts. The values outside the given diameter will be set to zero.
         To implement the cropping one can use the range [0.7-1.0] or set to 2.0 when no cropping required.
-    neglog: bool
-        Take negative logarithm on input data to convert to attenuation coefficient or a density of the scanned object. Defaults to False,
-        assuming that the negative log is taken either in normalisation procedure on with Paganin filter application.
 
     Returns
     -------
@@ -258,7 +244,7 @@ def LPRec3d_tomobar(
     )
 
     reconstruction = RecToolsCP.FOURIER_INV(
-        _take_neg_log(data) if neglog else data,
+        data,
         recon_mask_radius=recon_mask_radius,
         data_axes_labels_order=input_data_axis_labels,
         filter_type=filter_type,
@@ -282,7 +268,6 @@ def SIRT3d_tomobar(
     recon_mask_radius: float = 0.95,
     iterations: int = 300,
     nonnegativity: bool = True,
-    neglog: bool = False,
     gpu_id: int = 0,
 ) -> cp.ndarray:
     """
@@ -312,10 +297,7 @@ def SIRT3d_tomobar(
     iterations : int
         The number of SIRT iterations.
     nonnegativity : bool
-        Impose nonnegativity constraint on reconstructed image.
-    neglog: bool
-        Take negative logarithm on input data to convert to attenuation coefficient or a density of the scanned object. Defaults to False,
-        assuming that the negative log is taken either in normalisation procedure on with Paganin filter application.
+        Impose nonnegativity constraint on the reconstructed image.
     gpu_id : int
         A GPU device index to perform operation on.
 
@@ -336,7 +318,7 @@ def SIRT3d_tomobar(
     )
 
     _data_ = {
-        "projection_norm_data": _take_neg_log(data) if neglog else data,
+        "projection_norm_data": data,
         "data_axes_labels_order": input_data_axis_labels,
     }  # data dictionary
     _algorithm_ = {
@@ -359,7 +341,6 @@ def CGLS3d_tomobar(
     recon_mask_radius: float = 0.95,
     iterations: int = 20,
     nonnegativity: bool = True,
-    neglog: bool = False,
     gpu_id: int = 0,
 ) -> cp.ndarray:
     """
@@ -390,9 +371,6 @@ def CGLS3d_tomobar(
         The number of CGLS iterations.
     nonnegativity : bool
         Impose nonnegativity constraint on reconstructed image.
-    neglog: bool
-        Take negative logarithm on input data to convert to attenuation coefficient or a density of the scanned object. Defaults to False,
-        assuming that the negative log is taken either in normalisation procedure on with Paganin filter application.
     gpu_id : int, optional
         A GPU device index to perform operation on.
 
@@ -407,7 +385,7 @@ def CGLS3d_tomobar(
     )
 
     _data_ = {
-        "projection_norm_data": _take_neg_log(data) if neglog else data,
+        "projection_norm_data": data,
         "data_axes_labels_order": input_data_axis_labels,
     }  # data dictionary
     _algorithm_ = {
@@ -435,7 +413,6 @@ def FISTA3d_tomobar(
     regularisation_iterations: int = 50,
     regularisation_half_precision: bool = True,
     nonnegativity: bool = True,
-    neglog: bool = False,
     gpu_id: int = 0,
 ) -> cp.ndarray:
     """
@@ -474,9 +451,6 @@ def FISTA3d_tomobar(
         Perform faster regularisation computation in half-precision with a very minimal sacrifice in quality.
     nonnegativity : bool
         Impose nonnegativity constraint on the reconstructed image.
-    neglog: bool
-        Take negative logarithm on input data to convert to attenuation coefficient or a density of the scanned object. Defaults to False,
-        assuming that the negative log is taken either in normalisation procedure on with Paganin filter application.
     gpu_id : int
         A GPU device index to perform operation on.
 
@@ -490,7 +464,7 @@ def FISTA3d_tomobar(
     )
 
     _data_ = {
-        "projection_norm_data": _take_neg_log(data) if neglog else data,
+        "projection_norm_data": data,
         "OS_number": subsets_number,
         "data_axes_labels_order": input_data_axis_labels,
     }
@@ -649,24 +623,6 @@ def _instantiate_iterative_recon_class(
     return RecToolsCP
 
 
-def _take_neg_log(data: cp.ndarray) -> cp.ndarray:
-    """Taking negative log"""
-    data[data <= 0] = 1
-    data = -cp.log(data)
-    data[cp.isnan(data)] = 6.0
-    data[cp.isinf(data)] = 0
-    return data
-
-
-def _take_neg_log_np(data: np.ndarray) -> np.ndarray:
-    """Taking negative log"""
-    data[data <= 0] = 1
-    data = -np.log(data)
-    data[np.isnan(data)] = 6.0
-    data[np.isinf(data)] = 0
-    return data
-
-
 def __estimate_detectorHoriz_padding(detX_size) -> int:
     det_half = detX_size // 2
     padded_value_exact = int(np.sqrt(2 * (det_half**2))) - det_half