Merge pull request #196 from DiamondLightSource/neglog

adding neglog as a parameter of each reconstruction module

Author: dkazanc

docs/source/doc-conda-requirements.yml

@@ -3,7 +3,7 @@ channels:
   - conda-forge
 dependencies:
   - python>=3.10
-  - sphinx
+  - sphinx>=8.0,<8.2.0
   - sphinx-book-theme
   - nbsphinx
   - pandoc
@@ -13,4 +13,4 @@ dependencies:
   - sphinxcontrib-bibtex
   - pyyaml
   - ipython
-  - ghp-import
+  - ghp-import

httomolibgpu/recon/algorithm.py

@@ -54,9 +54,10 @@ def FBP(
     data: cp.ndarray,
     angles: np.ndarray,
     center: Optional[float] = None,
-    filter_freq_cutoff: Optional[float] = 0.35,
+    filter_freq_cutoff: float = 0.35,
     recon_size: Optional[int] = None,
-    recon_mask_radius: Optional[float] = 0.95,
+    recon_mask_radius: float = 0.95,
+    neglog: bool = False,
     gpu_id: int = 0,
 ) -> cp.ndarray:
     """
@@ -81,6 +82,9 @@ def FBP(
         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.
         It is recommended to keep the value in the range [0.7-1.0].
+    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.
 
@@ -94,7 +98,7 @@ def FBP(
     )
 
     reconstruction = RecToolsCP.FBP(
-        data,
+        _take_neg_log(data) if neglog else data,
         cutoff_freq=filter_freq_cutoff,
         recon_mask_radius=recon_mask_radius,
         data_axes_labels_order=input_data_axis_labels,
@@ -110,6 +114,7 @@ def LPRec(
     center: Optional[float] = None,
     recon_size: Optional[int] = None,
     recon_mask_radius: Optional[float] = 0.95,
+    neglog: bool = False,
 ) -> cp.ndarray:
     """
     Fourier direct inversion in 3D on unequally spaced (also called as Log-Polar) grids using
@@ -131,6 +136,9 @@ def LPRec(
         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.
         It is recommended to keep the value in the range [0.7-1.0].
+    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
     -------
@@ -140,7 +148,7 @@ def LPRec(
     RecToolsCP = _instantiate_direct_recon_class(data, angles, center, recon_size, 0)
 
     reconstruction = RecToolsCP.FOURIER_INV(
-        data,
+        _take_neg_log(data) if neglog else data,
         recon_mask_radius=recon_mask_radius,
         data_axes_labels_order=input_data_axis_labels,
     )
@@ -156,12 +164,14 @@ def SIRT(
     recon_size: Optional[int] = None,
     iterations: Optional[int] = 300,
     nonnegativity: Optional[bool] = True,
+    neglog: bool = False,
     gpu_id: int = 0,
 ) -> cp.ndarray:
     """
     Perform Simultaneous Iterative Recostruction Technique (SIRT) using ASTRA toolbox :cite:`van2016fast` and
     ToMoBAR :cite:`kazantsev2020tomographic` wrappers.
-    This is 3D recon directly from a CuPy array while using ASTRA GPUlink capability.
+    This is 3D recon directly from a CuPy array while using ASTRA GPUlink capability to avoid host-device
+    transactions for projection and backprojection.
 
     Parameters
     ----------
@@ -178,6 +188,9 @@ def SIRT(
         The number of SIRT iterations.
     nonnegativity : bool, optional
         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.
 
@@ -187,11 +200,16 @@ def SIRT(
         The SIRT reconstructed volume as a CuPy array.
     """
     RecToolsCP = _instantiate_iterative_recon_class(
-        data, angles, center, recon_size, gpu_id, datafidelity="LS"
+        data,
+        angles,
+        center,
+        recon_size,
+        gpu_id,
+        datafidelity="LS",
     )
 
     _data_ = {
-        "projection_norm_data": data,
+        "projection_norm_data": _take_neg_log(data) if neglog else data,
         "data_axes_labels_order": input_data_axis_labels,
     }  # data dictionary
     _algorithm_ = {
@@ -211,12 +229,14 @@ def CGLS(
     recon_size: Optional[int] = None,
     iterations: Optional[int] = 20,
     nonnegativity: Optional[bool] = True,
+    neglog: bool = False,
     gpu_id: int = 0,
 ) -> cp.ndarray:
     """
-    Perform Congugate Gradient Least Squares (CGLS) using ASTRA toolbox :cite:`van2016fast` and
+    Perform Conjugate Gradient Least Squares (CGLS) using ASTRA toolbox :cite:`van2016fast` and
     ToMoBAR :cite:`kazantsev2020tomographic` wrappers.
-    This is 3D recon directly from a CuPy array while using ASTRA GPUlink capability.
+    This is 3D recon directly from a CuPy array while using ASTRA GPUlink capability to avoid host-device
+    transactions for projection and backprojection.
 
     Parameters
     ----------
@@ -233,6 +253,9 @@ def CGLS(
         The number of CGLS iterations.
     nonnegativity : bool, optional
         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.
 
@@ -246,15 +269,14 @@ def CGLS(
     )
 
     _data_ = {
-        "projection_norm_data": data,
+        "projection_norm_data": _take_neg_log(data) if neglog else data,
         "data_axes_labels_order": input_data_axis_labels,
     }  # data dictionary
     _algorithm_ = {"iterations": iterations, "nonnegativity": nonnegativity}
     reconstruction = RecToolsCP.CGLS(_data_, _algorithm_)
     cp._default_memory_pool.free_all_blocks()
     return cp.require(cp.swapaxes(reconstruction, 0, 1), requirements="C")
 
-
 ## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  ##
 def _instantiate_direct_recon_class(
     data: cp.ndarray,
@@ -329,3 +351,11 @@ def _instantiate_iterative_recon_class(
         device_projector=gpu_id,
     )
     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

tests/test_recon/test_algorithm.py

@@ -31,6 +31,25 @@ def test_reconstruct_FBP_1(data, flats, darks, ensure_clean_memory):
     assert recon_data.shape == (160, 128, 160)
 
 
+def test_reconstruct_FBP_1_neglog(data, flats, darks, ensure_clean_memory):
+    recon_data = FBP(
+        normalize_cupy(data, flats, darks, cutoff=10, minus_log=False),
+        np.linspace(0.0 * np.pi / 180.0, 180.0 * np.pi / 180.0, data.shape[0]),
+        79.5,
+        filter_freq_cutoff=1.1,
+        recon_mask_radius=None,
+        neglog=True,
+    )
+    assert recon_data.flags.c_contiguous
+    recon_data = recon_data.get()
+    assert_allclose(np.mean(recon_data), 0.000798, rtol=1e-07, atol=1e-6)
+    assert_allclose(np.mean(recon_data, axis=(0, 2)).sum(), 0.102106, rtol=1e-05)
+    assert_allclose(np.std(recon_data), 0.006293, rtol=1e-07, atol=1e-6)
+    assert_allclose(np.median(recon_data), -0.000555, rtol=1e-07, atol=1e-6)
+    assert recon_data.dtype == np.float32
+    assert recon_data.shape == (160, 128, 160)
+
+
 def test_reconstruct_FBP_2(data, flats, darks, ensure_clean_memory):
     recon_data = FBP(
         normalize_cupy(data, flats, darks, cutoff=20.5, minus_log=False),