Update raven filter and add example

Author: neon60

docs/source/examples/raven_filter_example.py

@@ -0,0 +1,76 @@
+import os
+import numpy as np
+import cupy as cp
+import scipy
+import pyfftw
+import pyfftw.interfaces.numpy_fft as fft
+import httomolibgpu
+from httomolibgpu.misc.raven_filter import raven_filter
+
+import matplotlib.pyplot as plt
+
+# Load the sinogram data
+path_lib = os.path.dirname(httomolibgpu.__file__)
+in_file = os.path.abspath(
+    os.path.join(path_lib, "..", "tests/test_data/", "3600proj_sino.npz")
+)
+l_infile = np.load(in_file)
+sinogram = l_infile["sinogram"]
+angles = l_infile["angles"]
+sinogram = cp.asarray(sinogram)
+
+sino_shape = sinogram.shape
+
+print("The shape of the sinogram is {}".format(cp.shape(sinogram)))
+
+# Make a numpy copy
+sinogram_padded = np.pad(sinogram.get(), 20, "edge")
+
+# GPU filter
+sinogram_gpu_filter = raven_filter(sinogram)
+
+# Size
+width1 = sino_shape[1] + 2 * 20
+height1 = sino_shape[0] + 2 * 20
+
+# Parameters
+v0 = 2
+u0 = 20
+n = 2
+
+# Generate filter function
+centerx = np.ceil(width1 / 2.0) - 1.0
+centery = np.int16(np.ceil(height1 / 2.0) - 1)
+row1 = centery - v0
+row2 = centery + v0 + 1
+listx = np.arange(width1) - centerx
+filtershape = 1.0 / (1.0 + np.power(listx / u0, 2 * n))
+filtershapepad2d = np.zeros((row2 - row1, filtershape.size))
+filtershapepad2d[:] = np.float64(filtershape)
+filtercomplex = filtershapepad2d + filtershapepad2d * 1j
+
+# Generate filter objects
+a = pyfftw.empty_aligned((height1, width1), dtype='complex128', n=16)
+b = pyfftw.empty_aligned((height1, width1), dtype='complex128', n=16)
+c = pyfftw.empty_aligned((height1, width1), dtype='complex128', n=16)
+d = pyfftw.empty_aligned((height1, width1), dtype='complex128', n=16)
+fft_object  = pyfftw.FFTW(a, b, axes=(0, 1))
+ifft_object = pyfftw.FFTW(c, d, axes=(0, 1), direction='FFTW_BACKWARD')
+
+sino = fft.fftshift(fft_object(sinogram_padded))
+sino[row1:row2] = sino[row1:row2] * filtercomplex
+sino = ifft_object(fft.ifftshift(sino))
+
+plt.figure()
+
+#subplot(r,c) provide the no. of rows and columns
+f, axarr = plt.subplots(2,2) 
+
+# use the created array to output your multiple images. In this case I have stacked 4 images vertically
+axarr[0, 0].imshow(sinogram_padded)
+axarr[0, 1].imshow(sinogram_padded - sinogram_gpu_filter.get().real)
+axarr[1, 0].imshow(sinogram_padded - sino.real)
+axarr[1, 1].imshow(sinogram_gpu_filter.get().real - sino.real)
+
+plt.show()
+

httomolibgpu/cuda_kernels/raven_filter.cu

@@ -1,13 +1,36 @@
+#include <cupy/complex.cuh>
 
-template <typename Type, int diameter>
-__global__ void raven__filter_kernel3d(Type *in, int Z, int M, int N) {
+extern "C" __global__ void 
+raven_filter(complex<float> *input, complex<float> *output, int width1, int height1, int u0, int n, int v0) {
 
-  const long i = blockDim.x * blockIdx.x + threadIdx.x;
-  const long j = blockDim.y * blockIdx.y + threadIdx.y;
-  const long k = blockDim.z * blockIdx.z + threadIdx.z;
+  int centerx = (width1 + 1) / 2 - 1;
+  int centery = (height1 + 1) / 2 - 1;
 
-  if (i >= N || j >= M || k >= Z)
+  int px = threadIdx.x + blockIdx.x * blockDim.x;
+  int py = threadIdx.y + blockIdx.y * blockDim.y;
+
+  if (px >= width1)
+    return;
+  if (py >= height1)
     return;
 
-  
+  complex<float> value = input[py * width1 + px];
+  if( py >= (centery - v0) && py <= (centery + v0 + 1) ) {
+    
+    double base = (px - centerx) / u0;
+    double filtered_value = base;
+    for( int i = 1; i < 2 * n; i++ )
+      filtered_value *= base;
+
+    filtered_value = 1.0f / (1.0 + filtered_value);
+    value *= complex<float>(filtered_value, filtered_value);
+  }
+
+  // ifftshifting positions
+  int xshift = (width1 + 1) / 2;
+  int yshift = (height1 + 1) / 2;
+  int outX = (px + xshift) % width1;
+  int outY = (py + yshift) % height1;
+
+  output[outY * width1 + outX] = value;
 }

httomolibgpu/misc/raven_filter.py

@@ -49,12 +49,14 @@
 ]
 
 
-def raven_filter_savu(
+def raven_filter(
     data: cp.ndarray,
-    kernel_size: int = 3,
-    pad_y: int = 100,
-    pad_x: int = 100,
+    pad_y: int = 20,
+    pad_x: int = 20,
     pad_method: str = "edge",
+    uvalue: int = 20,
+    nvalue: int = 4,
+    vvalue: int = 2,
 ) -> cp.ndarray:
     """
     Applies raven filter to a 3D CuPy array. For more detailed information, see :ref:`method_raven_filter`.
@@ -63,8 +65,6 @@ def raven_filter_savu(
     ----------
     data : cp.ndarray
         Input CuPy 3D array either float32 or uint16 data type.
-    kernel_size : int, optional
-        The size of the filter's kernel (a diameter).
 
     pad_y : int, optional
         Pad the top and bottom of projections.
@@ -75,10 +75,19 @@ def raven_filter_savu(
     pad_method : str, optional
         Numpy pad method to use.
 
+    uvalue : int, optional
+        The shape of filter.
+
+    nvalue : int, optional
+        The shape of filter.
+
+    vvalue : int, optional
+        The number of rows to be applied the filter
+
     Returns
     -------
     ndarray
-        Raven filtered 3D CuPy array either float32 or uint16 data type.
+        Raven filtered 3D CuPy array in float32 data type.
 
     Raises
     ------
@@ -87,55 +96,40 @@ def raven_filter_savu(
     """
     input_type = data.dtype
 
-    if input_type not in ["float32", "uint16"]:
-        raise ValueError("The input data should be either float32 or uint16 data type")
-
-    if data.ndim == 3:
-        if 0 in data.shape:
-            raise ValueError("The length of one of dimensions is equal to zero")
-    else:
-        raise ValueError("The input array must be a 3D array")
-
-    if kernel_size not in [3, 5, 7, 9, 11, 13]:
-        raise ValueError("Please select a correct kernel size: 3, 5, 7, 9, 11, 13")
-
+    if input_type not in ["float32"]:
+        raise ValueError("The input data should be float32")
 
-    dz_orig, dy_orig, dx_orig = data.shape
+    # if data.ndim != 3:
+    #     raise ValueError("only 3D data is supported")
 
-    padded_data, pad_tup = cp.pad(data, fftpad, "edge")
-    dz, dy, dx = padded_data.shape
+    # Padding of the data
+    padded_data = cp.pad(data, ((pad_y, pad_y), (pad_x, pad_x)), mode=pad_method)
 
-    # 3D FFT of data
-    padded_data = cp.pad(data, ((0, 0), (pad_y, pad_y), (pad_x, pad_x)), mode=pad_method)
+    # FFT and shift of data
     fft_data = fft2(padded_data, axes=(-2, -1), overwrite_x=True)
     fft_data_shifted = fftshift(fft_data)
 
     # Setup various values for the filter
-    _, height, width = data.shape
+    height, width = data.shape
 
     height1 = height + 2 * pad_y
     width1 = width + 2 * pad_x
 
-    # raven
-    kernel_args = "raven_general_kernel3d<{0}, {1}>".format(
-        "float" if input_type == "float32" else "unsigned short", kernel_size
-    )
-    block_x = 128
     # setting grid/block parameters
+    block_x = 128
     block_dims = (block_x, 1, 1)
-    grid_x = (dx + block_x - 1) // block_x
-    grid_y = dy
-    grid_z = dz
-    grid_dims = (grid_x, grid_y, grid_z)
-    params = (fft_data_shifted, dz, dy, dx)
-
-    raven_module = load_cuda_module("raven_kernel", name_expressions=[kernel_args])
-    raven_filt = raven_module.get_function(kernel_args)
-
+    grid_x = (width1 + block_x - 1) // block_x
+    grid_y = height1
+    grid_dims = (grid_x, grid_y, 1)
+    params = (fft_data_shifted, fft_data, width1, height1, uvalue, nvalue, vvalue)
+
+    raven_module = load_cuda_module("raven_filter")
+    raven_filt = raven_module.get_function("raven_filter")
+    
     raven_filt(grid_dims, block_dims, params)
-
-    fft_data = fftshift(fft_data_shifted)
-
-    data = ifft2(fft_data, axes=(-2, -1), overwrite_x=True, norm="forward")
+    
+    # raven_fil already doing ifftshifting
+    # fft_data = ifftshift(fft_data_shifted)
+    data = ifft2(fft_data, axes=(-2, -1), overwrite_x=True)
 
     return data

httomolibgpu/prep/stripe.py

@@ -29,7 +29,10 @@
 from unittest.mock import Mock
 
 if cupy_run:
-    from cupyx.scipy.ndimage import median_filter, binary_dilation, raven_filter, uniform_filter1d
+    from cupyx.scipy.ndimage import median_filter, binary_dilation, uniform_filter1d
+    from httomolibgpu.misc.raven_filter import (
+        raven_filter,
+    )
 else:
     median_filter = Mock()
     binary_dilation = Mock()