adding more zenodo tests and update default FBP parameters

Author: dkazanc

httomolibgpu/recon/algorithm.py

@@ -54,15 +54,15 @@ def FBP(
     data: cp.ndarray,
     angles: np.ndarray,
     center: Optional[float] = None,
-    filter_freq_cutoff: Optional[float] = 1.1,
+    filter_freq_cutoff: Optional[float] = 0.35,
     recon_size: Optional[int] = None,
-    recon_mask_radius: Optional[float] = None,
+    recon_mask_radius: Optional[float] = 0.95,
     gpu_id: int = 0,
 ) -> cp.ndarray:
     """
     Perform Filtered Backprojection (FBP) reconstruction using ASTRA toolbox :cite:`van2016fast` and
     ToMoBAR :cite:`kazantsev2020tomographic` wrappers.
-    This is a 3D recon from a CuPy array directly and a custom built filter.
+    This is a 3D recon from the CuPy array directly and using a custom built SINC filter for filtration in Fourier space.
 
     Parameters
     ----------
@@ -72,16 +72,16 @@ def FBP(
         An array of angles given in radians.
     center : float, optional
         The center of rotation (CoR).
-    filter_freq_cutoff : float, optional
-        Cutoff frequency parameter for the sinc filter, the lowest values produce more crispy but noisy reconstruction.
+    filter_freq_cutoff : float
+        Cutoff frequency parameter for the SINC filter, the lower values produce better contrast but noisy reconstruction.
     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.
-    recon_mask_radius: float, optional
+    recon_mask_radius: float
         The radius of the circular mask that applies to the reconstructed slice in order to crop
-        out some undesirable artifacts. The values outside the diameter will be set to zero.
-        None by default, to see the effect of the mask try setting the value in the range [0.7-1.0].
-    gpu_id : int, optional
+        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].
+    gpu_id : int
         A GPU device index to perform operation on.
 
     Returns
@@ -109,7 +109,7 @@ def LPRec(
     angles: np.ndarray,
     center: Optional[float] = None,
     recon_size: Optional[int] = None,
-    recon_mask_radius: Optional[float] = None,
+    recon_mask_radius: Optional[float] = 0.95,
 ) -> cp.ndarray:
     """
     Fourier direct inversion in 3D on unequally spaced (also called as Log-Polar) grids using
@@ -127,10 +127,10 @@ def LPRec(
     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.
-    recon_mask_radius: float, optional
+    recon_mask_radius: float
         The radius of the circular mask that applies to the reconstructed slice in order to crop
-        out some undesirable artifacts. The values outside the diameter will be set to zero.
-        None by default, to see the effect of the mask try setting the value in the range [0.7-1.0].
+        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].
 
     Returns
     -------

tests/test_misc/test_morph.py

@@ -6,22 +6,6 @@
 from numpy.testing import assert_allclose
 from httomolibgpu.misc.morph import sino_360_to_180, data_resampler
 
-
-# @pytest.mark.parametrize("rotation", ["left", "right"])
-# def test_sino_360_to_180_realdata(ensure_clean_memory, sino3600, rotation):
-#     shape_new = (3601, 3, 2560)
-#     data3d = cp.zeros(shape_new, dtype=np.float32)
-#     data3d[:, 0, :] = sino3600
-#     data3d[:, 1, :] = sino3600
-#     data3d[:, 2, :] = sino3600
-
-#     sino_360_to_180(data3d, overlap=900, rotation=rotation)
-
-#     assert data3d.shape == shape_new
-#     assert data3d.dtype == np.float32
-#     assert data3d.flags.c_contiguous
-
-
 @pytest.mark.parametrize(
     "overlap, rotation",
     [

tests/test_recon/test_algorithm.py

@@ -19,6 +19,7 @@ def test_reconstruct_FBP_1(data, flats, darks, ensure_clean_memory):
         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,
     )
     assert recon_data.flags.c_contiguous
     recon_data = recon_data.get()
@@ -36,6 +37,7 @@ def test_reconstruct_FBP_2(data, flats, darks, ensure_clean_memory):
         np.linspace(5.0 * np.pi / 360.0, 180.0 * np.pi / 360.0, data.shape[0]),
         15.5,
         filter_freq_cutoff=1.1,
+        recon_mask_radius=None,
     )
 
     recon_data = recon_data.get()

tests/test_recon/test_rotation.py

@@ -93,17 +93,6 @@ def test_find_center_vo_performance():
     assert "performance in ms" == duration_ms
 
 
-def test_find_center_360_ones():
-    mat = cp.ones(shape=(100, 100, 100), dtype=cp.float32)
-
-    (cor, overlap, side, overlap_position) = find_center_360(mat)
-
-    assert_allclose(cor, 5.0)
-    assert_allclose(overlap, 12.0)
-    assert side == 0
-    assert_allclose(overlap_position, 7.0)
-
-
 def test_find_center_360_data(data):
     eps = 1e-5
     (cor, overlap, side, overlap_pos) = find_center_360(data, norm=True, denoise=False)

zenodo-tests/conftest.py

@@ -3,9 +3,6 @@
 import numpy as np
 import pytest
 
-CUR_DIR = os.path.abspath(os.path.dirname(__file__))
-
-
 def pytest_addoption(parser):
     parser.addoption(
         "--performance",
@@ -34,6 +31,9 @@ def pytest_collection_modifyitems(config, items):
                 item.add_marker(skip_perf)
 
 
+#CUR_DIR = os.path.abspath(os.path.dirname(__file__))
+CUR_DIR = "/dls/science/users/kjy41806/zenodo-tests/"
+
 @pytest.fixture(scope="session")
 def test_data_path():
     return os.path.join(CUR_DIR, "large_data_archive")

zenodo-tests/test_misc/__init__.py

@@ -0,0 +1,31 @@
+import time
+import cupy as cp
+import numpy as np
+from cupy.cuda import nvtx
+import pytest
+from numpy.testing import assert_allclose
+from httomolibgpu.misc.morph import sino_360_to_180
+from httomolibgpu.prep.normalize import normalize
+
+
+def test_sino_360_to_180_i13_dataset1(i13_dataset1, ensure_clean_memory):
+
+    projdata = i13_dataset1[0]
+    flats = i13_dataset1[2]
+    darks = i13_dataset1[3]
+    del i13_dataset1
+    
+    data_normalised = normalize(projdata, flats, darks, minus_log=True)
+    del flats, darks, projdata
+    ensure_clean_memory
+
+    stiched_data_180degrees = sino_360_to_180(data_normalised, overlap = 473.822265625, rotation = 'right')
+    stiched_data_180degrees = stiched_data_180degrees.get()
+
+    assert_allclose(np.sum(stiched_data_180degrees), 28512826.0, rtol=1e-07, atol=1e-6)
+    assert stiched_data_180degrees.shape == (3000, 10, 4646)
+    assert stiched_data_180degrees.dtype == np.float32
+    assert stiched_data_180degrees.flags.c_contiguous
+
+
+

zenodo-tests/test_misc/test_morph.py

@@ -0,0 +1,138 @@
+import cupy as cp
+import numpy as np
+import pytest
+from cupy.cuda import nvtx
+import time
+
+from httomolibgpu.prep.normalize import normalize
+from httomolibgpu.recon.algorithm import (
+    FBP,
+    LPRec,
+)
+from numpy.testing import assert_allclose
+import time
+import pytest
+from cupy.cuda import nvtx
+
+def test_reconstruct_FBP_i12_dataset1(i12_dataset1, ensure_clean_memory):
+
+    projdata = i12_dataset1[0]
+    angles = i12_dataset1[1]
+    flats = i12_dataset1[2]
+    darks = i12_dataset1[3]
+    del i12_dataset1
+
+    data_normalised = normalize(projdata, flats, darks, minus_log=True)
+    del flats, darks, projdata
+    ensure_clean_memory
+
+    recon_data = FBP(
+        data_normalised,
+        np.deg2rad(angles),
+        center=1253.75,
+        filter_freq_cutoff=0.35,     
+    )
+    assert recon_data.flags.c_contiguous
+    recon_data = recon_data.get()
+    assert_allclose(np.sum(recon_data), 46569.395, rtol=1e-07, atol=1e-6)
+    assert recon_data.dtype == np.float32
+    assert recon_data.shape == (2560, 50, 2560)
+
+@pytest.mark.perf
+def test_FBP_performance_i13_dataset2(i13_dataset2, ensure_clean_memory):
+    dev = cp.cuda.Device()
+    projdata = i13_dataset2[0]
+    angles = i13_dataset2[1]
+    flats = i13_dataset2[2]
+    darks = i13_dataset2[3]
+    del i13_dataset2
+
+    data_normalised = normalize(projdata, flats, darks, minus_log=True)
+    del flats, darks, projdata
+    ensure_clean_memory
+
+    # cold run first
+    FBP(
+        data_normalised,
+        np.deg2rad(angles),
+        center=1253.75,
+        filter_freq_cutoff=0.35,
+        )
+    dev.synchronize()
+
+    start = time.perf_counter_ns()
+    nvtx.RangePush("Core")
+    for _ in range(10):
+        FBP(
+        data_normalised,
+        np.deg2rad(angles),
+        center=1286.25,
+        filter_freq_cutoff=0.35,
+        )
+    nvtx.RangePop()
+    dev.synchronize()
+    duration_ms = float(time.perf_counter_ns() - start) * 1e-6 / 10
+
+    assert "performance in ms" == duration_ms
+
+
+def test_reconstruct_LPREC_i13_dataset2(i13_dataset2, ensure_clean_memory):
+    
+    projdata = i13_dataset2[0]
+    angles = i13_dataset2[1]
+    flats = i13_dataset2[2]
+    darks = i13_dataset2[3]
+    del i13_dataset2
+
+    data_normalised = normalize(projdata, flats, darks, minus_log=True)
+    del flats, darks, projdata
+    ensure_clean_memory
+
+
+    recon_data = LPRec(
+        data=data_normalised,
+        angles=np.deg2rad(angles),
+        center=1286.25,
+    )
+    assert recon_data.flags.c_contiguous
+    recon_data = recon_data.get()
+
+    assert_allclose(np.sum(recon_data), 2044.9531, rtol=1e-07, atol=1e-6)
+    assert recon_data.dtype == np.float32
+    assert recon_data.shape == (2560, 10, 2560)
+
+
+@pytest.mark.perf
+def test_LPREC_performance_i13_dataset2(i13_dataset2, ensure_clean_memory):
+    dev = cp.cuda.Device()
+    projdata = i13_dataset2[0]
+    angles = i13_dataset2[1]
+    flats = i13_dataset2[2]
+    darks = i13_dataset2[3]
+    del i13_dataset2
+
+    data_normalised = normalize(projdata, flats, darks, minus_log=True)
+    del flats, darks, projdata
+    ensure_clean_memory
+
+    # cold run first
+    LPRec(
+        data=data_normalised,
+        angles=np.deg2rad(angles),
+        center=1286.25,
+    )
+    dev.synchronize()
+
+    start = time.perf_counter_ns()
+    nvtx.RangePush("Core")
+    for _ in range(10):
+        LPRec(
+        data=data_normalised,
+        angles=np.deg2rad(angles),
+        center=1286.25,
+        )
+    nvtx.RangePop()
+    dev.synchronize()
+    duration_ms = float(time.perf_counter_ns() - start) * 1e-6 / 10
+
+    assert "performance in ms" == duration_ms

zenodo-tests/test_recon/test_algorithm.py

@@ -133,6 +133,7 @@ def test_center_360_i13_dataset1(i13_dataset1, ensure_clean_memory):
     cor, overlap, side, overlap_position = find_center_360(data_normalised)
 
     assert int(cor) == 2322
+    assert int(overlap) == 473 # actual 473.822265625
     assert cor.dtype == np.float64