adding parameterised tests for all stripe algorithms using zenodo data

Author: dkazanc

.scripts/download_zenodo.py

@@ -18,15 +18,15 @@ def calculate_md5(filename):
 
 def download_zenodo_files(output_dir: Path):
     """
-    Download all files from Zenodo record 14833590 and verify their checksums.
+    Download all files from Zenodo record 14938787 and verify their checksums.
 
     Args:
         output_dir: Directory where files should be downloaded
     """
     try:
-        print("Fetching files from Zenodo record 14833590...")
+        print("Fetching files from Zenodo record 14938787...")
         with urllib.request.urlopen(
-            "https://zenodo.org/api/records/14833590"
+            "https://zenodo.org/api/records/14938787"
         ) as response:
             data = json.loads(response.read())
 

httomolibgpu/misc/morph.py

@@ -112,7 +112,7 @@ def data_resampler(
 
     Parameters
     ----------
-    data : cp.ndarray 
+    data : cp.ndarray
         3d cupy array.
     newshape : list
         2d list that defines the 2D slice shape of new shape data.

httomolibgpu/prep/stripe.py

@@ -317,9 +317,8 @@ def _rs_dead(sinogram, snr, size, matindex, norm=True):
         ids = cp.searchsorted(listx, listxmiss)
         weights = (listxmiss - listx[ids - 1]) / (listx[ids] - listx[ids - 1])
         # direct interpolation without making an extra copy
-        sinogram[:, listxmiss] = (
-            sinogram[:, listx[ids - 1]] +
-            weights * (sinogram[:, listx[ids]] - sinogram[:, listx[ids - 1]])
+        sinogram[:, listxmiss] = sinogram[:, listx[ids - 1]] + weights * (
+            sinogram[:, listx[ids]] - sinogram[:, listx[ids - 1]]
         )
 
     # Remove residual stripes
@@ -416,7 +415,7 @@ def raven_filter(
     # Removing padding
     data = data[pad_y : height - pad_y, :, pad_x : width - pad_x].real
 
-    return data
+    return cp.require(data, requirements="C")
 
 
 def _create_matindex(nrow, ncol):

httomolibgpu/recon/algorithm.py

@@ -83,8 +83,8 @@ def FBP(
         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. 
+        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.
 
@@ -137,8 +137,8 @@ def LPRec(
         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.         
+        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
     -------
@@ -189,8 +189,8 @@ def SIRT(
     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.              
+        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.
 
@@ -254,8 +254,8 @@ def CGLS(
     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.              
+        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.
 
@@ -277,6 +277,7 @@ def CGLS(
     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,
@@ -352,10 +353,11 @@ def _instantiate_iterative_recon_class(
     )
     return RecToolsCP
 
+
 def _take_neg_log(data: cp.ndarray) -> cp.ndarray:
     """Taking negative log"""
-    data[data<=0] = 1
+    data[data <= 0] = 1
     data = -cp.log(data)
     data[cp.isnan(data)] = 6.0
     data[cp.isinf(data)] = 0
-    return data
+    return data

tests/test_misc/test_morph.py

@@ -6,6 +6,7 @@
 from numpy.testing import assert_allclose
 from httomolibgpu.misc.morph import sino_360_to_180, data_resampler
 
+
 @pytest.mark.parametrize(
     "overlap, rotation",
     [

tests/test_prep/stripe_cpu_reference.py

@@ -2,26 +2,29 @@
 import pyfftw
 import pyfftw.interfaces.numpy_fft as fft
 
+
 def raven_filter_numpy(
-        sinogram,
-        uvalue: int = 20,
-        nvalue: int = 4,
-        vvalue: int = 2,
-        pad_y: int = 20,
-        pad_x: int = 20,
-        pad_method: str = "edge"):
-    
+    sinogram,
+    uvalue: int = 20,
+    nvalue: int = 4,
+    vvalue: int = 2,
+    pad_y: int = 20,
+    pad_x: int = 20,
+    pad_method: str = "edge",
+):
     # Parameters
     v0 = vvalue
     n = nvalue
     u0 = uvalue
 
     # Make a padded copy
-    sinogram_padded = np.pad(sinogram, ((pad_y,pad_y), (0, 0), (pad_x,pad_x)), pad_method)
-    
+    sinogram_padded = np.pad(
+        sinogram, ((pad_y, pad_y), (0, 0), (pad_x, pad_x)), pad_method
+    )
+
     # Size
     height, images, width = sinogram_padded.shape
-    
+
     # Generate filter function
     centerx = np.ceil(width / 2.0) - 1.0
     centery = np.int16(np.ceil(height / 2.0) - 1)
@@ -32,19 +35,19 @@ def raven_filter_numpy(
     filtershapepad2d = np.zeros((row2 - row1, filtershape.size))
     filtershapepad2d[:] = np.float64(filtershape)
     filtercomplex = filtershapepad2d + filtershapepad2d * 1j
-    
+
     # Generate filter objects
-    a = pyfftw.empty_aligned((height, images, width), dtype='complex128', n=16)
-    b = pyfftw.empty_aligned((height, images, width), dtype='complex128', n=16)
-    c = pyfftw.empty_aligned((height, images, width), dtype='complex128', n=16)
-    d = pyfftw.empty_aligned((height, images, width), dtype='complex128', n=16)
-    fft_object  = pyfftw.FFTW(a, b, axes=(0, 2))
-    ifft_object = pyfftw.FFTW(c, d, axes=(0, 2), direction='FFTW_BACKWARD')
-    
+    a = pyfftw.empty_aligned((height, images, width), dtype="complex128", n=16)
+    b = pyfftw.empty_aligned((height, images, width), dtype="complex128", n=16)
+    c = pyfftw.empty_aligned((height, images, width), dtype="complex128", n=16)
+    d = pyfftw.empty_aligned((height, images, width), dtype="complex128", n=16)
+    fft_object = pyfftw.FFTW(a, b, axes=(0, 2))
+    ifft_object = pyfftw.FFTW(c, d, axes=(0, 2), direction="FFTW_BACKWARD")
+
     sino = fft.fftshift(fft_object(sinogram_padded), axes=(0, 2))
     for m in range(sino.shape[1]):
         sino[row1:row2, m] = sino[row1:row2, m] * filtercomplex
     sino = ifft_object(fft.ifftshift(sino, axes=(0, 2)))
-    sinogram = sino[pad_y:height-pad_y, :, pad_x:width-pad_x]
+    sinogram = sino[pad_y : height - pad_y, :, pad_x : width - pad_x]
 
-    return sinogram.real
+    return sinogram.real

tests/test_prep/test_stripe.py

@@ -14,6 +14,7 @@
 from numpy.testing import assert_allclose
 from .stripe_cpu_reference import raven_filter_numpy
 
+
 def test_remove_stripe_ti_on_data(data, flats, darks):
     # --- testing the CuPy implementation from TomoCupy ---#
     data = normalize(data, flats, darks, cutoff=10, minus_log=True)
@@ -53,6 +54,7 @@ def test_remove_stripe_ti_on_data(data, flats, darks):
 #         np.median(corrected_data), np.median(corrected_host_data), rtol=1e-6
 #     )
 
+
 def test_stripe_removal_sorting_cupy(data, flats, darks):
     # --- testing the CuPy port of TomoPy's implementation ---#
     data = normalize(data, flats, darks, cutoff=10, minus_log=True)
@@ -67,6 +69,7 @@ def test_stripe_removal_sorting_cupy(data, flats, darks):
     assert corrected_data.dtype == np.float32
     assert corrected_data.flags.c_contiguous
 
+
 def test_stripe_raven_cupy(data, flats, darks):
     # --- testing the CuPy port of TomoPy's implementation ---#
 
@@ -82,13 +85,16 @@ def test_stripe_raven_cupy(data, flats, darks):
     assert data_after_raven_gpu.dtype == np.float32
     assert data_after_raven_gpu.shape == data_after_raven_cpu.shape
 
+
 @pytest.mark.parametrize("uvalue", [20, 50, 100])
 @pytest.mark.parametrize("nvalue", [2, 4, 6])
 @pytest.mark.parametrize("vvalue", [2, 4])
 @pytest.mark.parametrize("pad_x", [0, 10, 20])
 @pytest.mark.parametrize("pad_y", [0, 10, 20])
 @cp.testing.numpy_cupy_allclose(rtol=0, atol=3e-01)
-def test_stripe_raven_parameters_cupy(ensure_clean_memory, xp, uvalue, nvalue, vvalue, pad_x, pad_y):
+def test_stripe_raven_parameters_cupy(
+    ensure_clean_memory, xp, uvalue, nvalue, vvalue, pad_x, pad_y
+):
     # because it's random, we explicitly seed and use numpy only, to match the data
     np.random.seed(12345)
     data = np.random.random_sample(size=(256, 5, 512)).astype(np.float32) * 2.0 + 0.001
@@ -97,14 +103,15 @@ def test_stripe_raven_parameters_cupy(ensure_clean_memory, xp, uvalue, nvalue, v
     if xp.__name__ == "numpy":
         results = raven_filter_numpy(
             data, uvalue=uvalue, nvalue=nvalue, vvalue=vvalue, pad_x=pad_x, pad_y=pad_y
-        ).astype(np.float32) 
+        ).astype(np.float32)
     else:
         results = raven_filter(
             data, uvalue=uvalue, nvalue=nvalue, vvalue=vvalue, pad_x=pad_x, pad_y=pad_y
         ).get()
 
     return xp.asarray(results)
 
+
 @pytest.mark.perf
 def test_stripe_removal_sorting_cupy_performance(ensure_clean_memory):
     data_host = (
@@ -154,6 +161,7 @@ def test_remove_stripe_ti_performance(ensure_clean_memory):
 
     assert "performance in ms" == duration_ms
 
+
 @pytest.mark.perf
 def test_raven_filter_performance(ensure_clean_memory):
     data_host = (
@@ -178,6 +186,7 @@ def test_raven_filter_performance(ensure_clean_memory):
 
     assert "performance in ms" == duration_ms
 
+
 def test_remove_all_stripe_on_data(data, flats, darks):
     # --- testing the CuPy implementation from TomoCupy ---#
     data = normalize(data, flats, darks, cutoff=10, minus_log=True)
@@ -190,8 +199,12 @@ def test_remove_all_stripe_on_data(data, flats, darks):
     )
     assert_allclose(np.median(data_after_stripe_removal), 0.015338, rtol=1e-04)
     assert_allclose(np.max(data_after_stripe_removal), 2.298123, rtol=1e-05)
-    assert_allclose(np.median(data_after_stripe_removal, axis=(1, 2)).sum(), 2.788046, rtol=1e-6)
-    assert_allclose(np.median(data_after_stripe_removal, axis=(0, 1)).sum(), 28.661312, rtol=1e-6)
+    assert_allclose(
+        np.median(data_after_stripe_removal, axis=(1, 2)).sum(), 2.788046, rtol=1e-6
+    )
+    assert_allclose(
+        np.median(data_after_stripe_removal, axis=(0, 1)).sum(), 28.661312, rtol=1e-6
+    )
 
     data = None  #: free up GPU memory
     # make sure the output is float32

zenodo-tests/conftest.py

@@ -89,6 +89,22 @@ def i12_dataset3(i12_dataset3_file):
     )
 
 
+@pytest.fixture(scope="session")
+def i12_dataset4_file(test_data_path):
+    in_file = os.path.join(test_data_path, "i12_dataset4.npz")
+    return np.load(in_file)
+
+
+@pytest.fixture
+def i12_dataset4(i12_dataset4_file):
+    return (
+        cp.asarray(i12_dataset4_file["projdata"]),
+        i12_dataset4_file["angles"],
+        cp.asarray(i12_dataset4_file["flats"]),
+        cp.asarray(i12_dataset4_file["darks"]),
+    )
+
+
 @pytest.fixture(scope="session")
 def i13_dataset1_file(test_data_path):
     in_file = os.path.join(test_data_path, "i13_dataset1.npz")