Merge pull request #200 from DiamondLightSource/zenodo_stripes

Zenodo stripes removal tests

Author: namannimmo10

.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

@@ -205,40 +205,12 @@ def remove_all_stripe(
     for m in range(data.shape[1]):
         sino = data[:, m, :]
         sino = _rs_dead(sino, snr, la_size, matindex)
-        sino = _rs_sort2(sino, sm_size, matindex, dim)
+        sino = _rs_sort(sino, sm_size, dim)
         sino = cp.nan_to_num(sino)
         data[:, m, :] = sino
     return data
 
 
-def _rs_sort2(sinogram, size, matindex, dim):
-    """
-    Remove stripes using the sorting technique.
-    """
-    sinogram = cp.transpose(sinogram)
-    matcomb = cp.asarray(cp.dstack((matindex, sinogram)))
-
-    # matsort = cp.asarray([row[row[:, 1].argsort()] for row in matcomb])
-    ids = cp.argsort(matcomb[:, :, 1], axis=1)
-    matsort = matcomb.copy()
-    matsort[:, :, 0] = cp.take_along_axis(matsort[:, :, 0], ids, axis=1)
-    matsort[:, :, 1] = cp.take_along_axis(matsort[:, :, 1], ids, axis=1)
-    if dim == 1:
-        matsort[:, :, 1] = median_filter(matsort[:, :, 1], (size, 1))
-    else:
-        matsort[:, :, 1] = median_filter(matsort[:, :, 1], (size, size))
-
-    # matsortback = cp.asarray([row[row[:, 0].argsort()] for row in matsort])
-
-    ids = cp.argsort(matsort[:, :, 0], axis=1)
-    matsortback = matsort.copy()
-    matsortback[:, :, 0] = cp.take_along_axis(matsortback[:, :, 0], ids, axis=1)
-    matsortback[:, :, 1] = cp.take_along_axis(matsortback[:, :, 1], ids, axis=1)
-
-    sino_corrected = matsortback[:, :, 1]
-    return cp.transpose(sino_corrected)
-
-
 def _mpolyfit(x, y):
     n = len(x)
     x_mean = cp.mean(x)
@@ -261,8 +233,6 @@ def _detect_stripe(listdata, snr):
     listsorted = cp.sort(listdata)[::-1]
     xlist = cp.arange(0, numdata, 1.0)
     ndrop = cp.int16(0.25 * numdata)
-    # (_slope, _intercept) = cp.polyfit(xlist[ndrop:-ndrop - 1],
-    #   listsorted[ndrop:-ndrop - 1], 1)
     (_slope, _intercept) = _mpolyfit(
         xlist[ndrop : -ndrop - 1], listsorted[ndrop : -ndrop - 1]
     )
@@ -293,11 +263,6 @@ def _rs_large(sinogram, snr, size, matindex, drop_ratio=0.1, norm=True):
     sinosmooth = median_filter(sinosort, (1, size))
     list1 = cp.mean(sinosort[ndrop : nrow - ndrop], axis=0)
     list2 = cp.mean(sinosmooth[ndrop : nrow - ndrop], axis=0)
-    # listfact = cp.divide(list1,
-    #                      list2,
-    #                      out=cp.ones_like(list1),
-    #                      where=list2 != 0)
-
     listfact = list1 / list2
 
     # Locate stripes
@@ -310,14 +275,12 @@ def _rs_large(sinogram, snr, size, matindex, drop_ratio=0.1, norm=True):
     sinogram1 = cp.transpose(sinogram)
     matcombine = cp.asarray(cp.dstack((matindex, sinogram1)))
 
-    # matsort = cp.asarray([row[row[:, 1].argsort()] for row in matcombine])
     ids = cp.argsort(matcombine[:, :, 1], axis=1)
     matsort = matcombine.copy()
     matsort[:, :, 0] = cp.take_along_axis(matsort[:, :, 0], ids, axis=1)
     matsort[:, :, 1] = cp.take_along_axis(matsort[:, :, 1], ids, axis=1)
 
     matsort[:, :, 1] = cp.transpose(sinosmooth)
-    # matsortback = cp.asarray([row[row[:, 0].argsort()] for row in matsort])
     ids = cp.argsort(matsort[:, :, 0], axis=1)
     matsortback = matsort.copy()
     matsortback[:, :, 0] = cp.take_along_axis(matsortback[:, :, 0], ids, axis=1)
@@ -330,12 +293,9 @@ def _rs_large(sinogram, snr, size, matindex, drop_ratio=0.1, norm=True):
 
 
 def _rs_dead(sinogram, snr, size, matindex, norm=True):
-    """
-    Remove unresponsive and fluctuating stripes.
-    """
+    """remove unresponsive and fluctuating stripes"""
     sinogram = cp.copy(sinogram)  # Make it mutable
     (nrow, _) = sinogram.shape
-    # sinosmooth = cp.apply_along_axis(uniform_filter1d, 0, sinogram, 10)
     sinosmooth = uniform_filter1d(sinogram, 10, axis=0)
 
     listdiff = cp.sum(cp.abs(sinogram - sinosmooth), axis=0)
@@ -344,22 +304,22 @@ def _rs_dead(sinogram, snr, size, matindex, norm=True):
     listfact = listdiff / listdiffbck
 
     listmask = _detect_stripe(listfact, snr)
+    del listfact
     listmask = binary_dilation(listmask, iterations=1).astype(listmask.dtype)
     listmask[0:2] = 0.0
     listmask[-2:] = 0.0
-    listx = cp.where(listmask < 1.0)[0]
-    listy = cp.arange(nrow)
-    matz = sinogram[:, listx]
 
+    listx = cp.where(listmask < 1.0)[0]
     listxmiss = cp.where(listmask > 0.0)[0]
+    del listmask
 
-    # finter = interpolate.interp2d(listx.get(), listy.get(), matz.get(), kind='linear')
     if len(listxmiss) > 0:
-        # sinogram_c[:, listxmiss.get()] = finter(listxmiss.get(), listy.get())
         ids = cp.searchsorted(listx, listxmiss)
-        sinogram[:, listxmiss] = matz[:, ids - 1] + (listxmiss - listx[ids - 1]) * (
-            matz[:, ids] - matz[:, ids - 1]
-        ) / (listx[ids] - listx[ids - 1])
+        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]]
+        )
 
     # Remove residual stripes
     if norm is True:
@@ -455,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