Update CudaSinoNormalization with division op

nabu/reconstruction/sinogram_cuda.py

@@ -188,17 +188,28 @@ class CudaSinoNormalization(SinoNormalization):
             }
         elif self.normalization_array is not None:
             normalization_array = self.normalization_array
-            n_angles, n_x = self.sinos_shape[-2:]
             # If normalization_array is 1D, make a 2D array by repeating the line
             if normalization_array.ndim == 1:
-                normalization_array = np.tile(normalization_array, (n_angles, 1))
-            if normalization_array.shape != (n_angles, n_x):
-                raise ValueError(
-                    "Expected normalization array to be of shape (n_angles, n_x) = (%d, %d)"
-                    % (n_angles, n_x)
-                )
+                normalization_array = np.tile(normalization_array, (self.n_angles, 1))
             self._d_normalization_array = garray.to_gpu(normalization_array.astype("f"))
-
+            if self.normalization_kind == "subtraction":
+                generic_op_val = 1
+            elif self.normalization_kind == "division":
+                generic_op_val = 3
+            self._norm_kernel = CudaKernel(
+                "inplace_generic_op_2Dby2D",
+                filename=get_cuda_srcfile("ElementOp.cu"),
+                signature="PPii",
+                options=["-DGENERIC_OP=%d" % generic_op_val]
+            )
+            self._norm_kernel_args = [
+                self._d_normalization_array, np.int32(self.n_angles), np.int32(self.n_x)
+            ]
+            blk = (32, 32, 1)
+            self._norm_kernel_kwargs = {
+                "block": blk,
+                "grid": (int(updiv(self.n_angles, blk[0])), int(updiv(self.n_x, blk[1])), 1)
+            }
 
     def _normalize_chebyshev(self, sinos):
         if sinos.flags.c_contiguous:
@@ -224,24 +235,24 @@ class CudaSinoNormalization(SinoNormalization):
 
 
     #
-    # Array subtraction
+    # Array subtraction/division
     #
 
-    def _normalize_subtract_array(self, sino):
+    def _normalize_op(self, sino):
         if sino.ndim == 2:
             # Things can go wrong if "sino" is a non-contiguous 2D array
             # But this should be handled outside this function, as the processing is in-place
-            sino -= self._d_normalization_array
+            self._norm_kernel(sino, *self._norm_kernel_args, **self._norm_kernel_kwargs)
         else:
             if sino.flags.forc:
                 # Contiguous 3D array. But pycuda wants the same shape for both operands.
                 for i in range(sino.shape[0]):
-                    sino[i] -= self._d_normalization_array
+                    self._norm_kernel(sino[i], *self._norm_kernel_args, **self._norm_kernel_kwargs)
             else:
                 # Non-contiguous 2D array. Make a temp. copy
                 for i in range(sino.shape[0]):
                     self._d_tmp[:] = sino[i][:]
-                    self._d_tmp -= self._d_normalization_array
+                    self._norm_kernel(self._d_tmp, *self._norm_kernel_args, **self._norm_kernel_kwargs)
                     sino[i][:] = self._d_tmp[:]
         return sino