Add cuda implementation

nabu/cuda/src/halftomo.cu

@@ -30,8 +30,6 @@ __global__ void halftomo_kernel(
         output[y * n_x2 + x] = sinogram[y * n_x + x];
     }
 
-
-
     // output[:, nx - d : nx] = (1 - weights) * sino[:n_a2, nx - d :]
     //                        + weights * sino[n_a2:, ::-1][:, d : 2 * d]
     else if (x < n_x) { // x in [n_x - d, n_x [
@@ -51,3 +49,49 @@ __global__ void halftomo_kernel(
     }
 
 }
+
+
+// Same as previous kernel, but now the CoR is outside the image support
+// i.e rotation_axis_position >= n_x
+// Weigting is different !
+__global__ void halftomo_kernel_cor_outside(
+    float* sinogram,
+    float* output,
+    float* weights,
+    int n_angles,
+    int n_x,
+    int rotation_axis_position
+) {
+    int x = blockDim.x * blockIdx.x + threadIdx.x;
+    int y = blockDim.y * blockIdx.y + threadIdx.y;
+
+    int n_a2 = n_angles / 2;
+    int d = rotation_axis_position - n_x;
+    int n_x2  = 2 * rotation_axis_position;
+
+    if ((x >= n_x2) || (y >= n_a2)) return;
+
+    // output[:, :nx] = sino[:n_a2, :nx]
+    if (x < n_x) {
+        output[y * n_x2 + x] = sinogram[y * n_x + x];
+    }
+
+    // output[:, nx : nx + d] = (1 - weights) * sino[:n_a2, -d :][:, ::-1]
+    else if (x < n_x + d) {
+        float w = weights[x - n_x];
+        // output[y * n_x2 + x] = (1.0f - w) * sinogram[y*n_x + n_x - 1 - (n_x - x)];
+        output[y * n_x2 + x] = (1.0f - w) * sinogram[y*n_x + 2 * n_x - 1 - x];
+    }
+
+    // output[:, nx + d : nx + 2*d] = weights * sino[n_a2:, ::-1][:, :d] = sino[n_a2:, -d-1:-n_x-1:-1]
+    else if (x < n_x + 2 * d) {
+        float w = weights[x - (n_x + d)];
+        output[y * n_x2 + x] = w * sinogram[(n_a2 + y)*n_x + x - 2*d];
+    }
+
+    // output[:, nx+2*d:] = sino[n_a2:, ::-1][:, 2 * d :] = sino[n_a2:, -2*d-1:-n_x-1:-1]
+    else {
+        output[y * n_x2 + x] = sinogram[(n_a2 + y)*n_x + (n_x2 - 1 - x)];
+    }
+
+}
\ No newline at end of file

nabu/preproc/sinogram_cuda.py

@@ -26,8 +26,12 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
     def _init_cuda_halftomo(self):
         if not(self.halftomo):
             return
+        if self._rot_center_int < self.n_x:
+            kernel_name = "halftomo_kernel"
+        else:
+            kernel_name = "halftomo_kernel_cor_outside"
         self.halftomo_kernel = CudaKernel(
-            "halftomo_kernel",
+            kernel_name,
             get_cuda_srcfile("halftomo.cu"),
             signature="PPPiii",
         )
@@ -40,7 +44,7 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
             1
         )
         d = self.n_x - rc # will have to be adapted for varying axis pos
-        self.halftomo_weights = np.linspace(0, 1, d, endpoint=True, dtype="f")
+        self.halftomo_weights = np.linspace(0, 1, abs(d), endpoint=True, dtype="f")
         self.d_halftomo_weights = garray.to_gpu(self.halftomo_weights)
         if self._halftomo_flip:
             self.xflip_kernel = CudaKernel(