Match with the CPU implementation

Author: neon60

docs/source/examples/raven_filter_example.py

@@ -23,32 +23,44 @@
 
 print("The shape of the sinogram is {}".format(cp.shape(sinogram)))
 
+# Parameters
+v0 = 2
+n = 4
+u0 = 20
+
 # Make a numpy copy
 sinogram_padded = np.pad(sinogram.get(), 20, "edge")
 
 # GPU filter
-sinogram_gpu_filter = raven_filter(sinogram)
+sinogram_gpu_filter = raven_filter(sinogram, u0, n, v0)
 
 # 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
+
+print(np.arange(width1))
+print(listx)
+print(listx / u0)
+
 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
 
+print("Centery: ", centery, "Row1: ", row1, "Row2: ", row2)
+
+# plt.figure()
+# plt.imshow(filtercomplex.real)
+# plt.title("Uncorrected image")
+# plt.show()
+
 # Generate filter objects
 a = pyfftw.empty_aligned((height1, width1), dtype='complex128', n=16)
 b = pyfftw.empty_aligned((height1, width1), dtype='complex128', n=16)

httomolibgpu/cuda_kernels/raven_filter.cu

@@ -3,8 +3,8 @@
 extern "C" __global__ void 
 raven_filter(complex<float> *input, complex<float> *output, int width1, int height1, int u0, int n, int v0) {
 
-  int centerx = (width1 + 1) / 2 - 1;
-  int centery = (height1 + 1) / 2 - 1;
+  int centerx = width1 / 2;
+  int centery = height1 / 2;
 
   int px = threadIdx.x + blockIdx.x * blockDim.x;
   int py = threadIdx.y + blockIdx.y * blockDim.y;
@@ -15,14 +15,15 @@ raven_filter(complex<float> *input, complex<float> *output, int width1, int heig
     return;
 
   complex<float> value = input[py * width1 + px];
-  if( py >= (centery - v0) && py <= (centery + v0 + 1) ) {
+  if( py >= (centery - v0) && py < (centery + v0 + 1) ) {
 
-    double base = (px - centerx) / u0;
-    double filtered_value = base;
+    // +1 needed to match with CPU implementation
+    float base = float(px - centerx + 1) / u0;
+    float power = base;
     for( int i = 1; i < 2 * n; i++ )
-      filtered_value *= base;
+      power *= base;
 
-    filtered_value = 1.0f / (1.0 + filtered_value);
+    float filtered_value = 1.f / (1.f + power);
     value *= complex<float>(filtered_value, filtered_value);
   }
 

httomolibgpu/misc/raven_filter.py

@@ -51,12 +51,12 @@
 
 def raven_filter(
     data: cp.ndarray,
-    pad_y: int = 20,
-    pad_x: int = 20,
-    pad_method: str = "edge",
     uvalue: int = 20,
     nvalue: int = 4,
     vvalue: int = 2,
+    pad_y: int = 20,
+    pad_x: int = 20,
+    pad_method: str = "edge",
 ) -> cp.ndarray:
     """
     Applies raven filter to a 3D CuPy array. For more detailed information, see :ref:`method_raven_filter`.