Merge pull request #184 from DiamondLightSource/vocentering_gpu_offload

Vocentering gpu offload

Author: dkazanc

.scripts/download_zenodo.py

@@ -19,15 +19,15 @@ def calculate_md5(filename):
 
 def download_zenodo_files(output_dir: Path):
     """
-    Download all files from Zenodo record 14627503 and verify their checksums.
+    Download all files from Zenodo record 14652312 and verify their checksums.
 
     Args:
         output_dir: Directory where files should be downloaded
     """
     try:
-        print("Fetching files from Zenodo record 14627503...")
+        print("Fetching files from Zenodo record 14652312...")
         with urllib.request.urlopen(
-            "https://zenodo.org/api/records/14627503"
+            "https://zenodo.org/api/records/14652312"
         ) as response:
             data = json.loads(response.read())
 

docs/source/examples/Cor_largesino.ipynb

@@ -30,5 +30,4 @@
 
     examples/pipeline1_FBP
     examples/pipeline2_iterative
-    examples/Cor_largesino
-    examples/DistortionCorr
+    examples/DistortionCorr

docs/source/index.rst

@@ -1,7 +1,7 @@
 #include <cupy/complex.cuh>
 
 extern "C" __global__ void
-shift_whole_shifts(const float *sino2, const float *sino3,
+shift_whole_shifts(const float *flip_sino, const float *comp_sino,
                    const float *__restrict__ list_shift, float *mat, int nx,
                    int nymat) {
   int xid = threadIdx.x + blockIdx.x * blockDim.x;
@@ -17,14 +17,14 @@ shift_whole_shifts(const float *sino2, const float *sino3,
   float frac_part = modf(shift_col, &int_part);
   if (abs(frac_part) > 1e-5f) {
     // we have a floating point shift, so we only roll in
-    // sino3, but we leave the rest for later using scipy
+    // comp_sino, but we leave the rest for later using scipy
     int shift_int =
         shift_col >= 0.0 ? int(ceil(shift_col)) : int(floor(shift_col));
     if (shift_int >= 0 && xid < shift_int) {
-      mat[zid * nymat * nx + yid * nx + xid] = sino3[yid * nx + xid];
+      mat[zid * nymat * nx + yid * nx + xid] = comp_sino[yid * nx + xid];
     }
     if (shift_int < 0 && xid >= nx + shift_int) {
-      mat[zid * nymat * nx + yid * nx + xid] = sino3[yid * nx + xid];
+      mat[zid * nymat * nx + yid * nx + xid] = comp_sino[yid * nx + xid];
     }
   } else {
     // we have an integer shift, so we can roll in directly
@@ -33,16 +33,16 @@ shift_whole_shifts(const float *sino2, const float *sino3,
     if (shift_int >= 0) {
       if (xid >= shift_int) {
         mat[zid * nymat * nx + yid * nx + xid] =
-            sino2[yid * nx + xid - shift_int];
+            flip_sino[yid * nx + xid - shift_int];
       } else {
-        mat[zid * nymat * nx + yid * nx + xid] = sino3[yid * nx + xid];
+        mat[zid * nymat * nx + yid * nx + xid] = comp_sino[yid * nx + xid];
       }
     } else {
       if (xid < nx + shift_int) {
         mat[zid * nymat * nx + yid * nx + xid] =
-            sino2[yid * nx + xid - shift_int];
+            flip_sino[yid * nx + xid - shift_int];
       } else {
-        mat[zid * nymat * nx + yid * nx + xid] = sino3[yid * nx + xid];
+        mat[zid * nymat * nx + yid * nx + xid] = comp_sino[yid * nx + xid];
       }
     }
   }

httomolibgpu/cuda_kernels/center_360_shifts.cu

@@ -48,4 +48,56 @@ extern "C" __global__ void generate_mask(const int ncol, const int nrow,
   }
 
   mask[j * (ncol/2+1) + outi] = outval;
+}
+
+extern "C" __global__ void generate_mask_full(const int ncol, const int nrow,
+                                              const int cen_col, const int cen_row,
+                                              const float du, const float dv,
+                                              const float radius, const float drop,
+                                              float *mask) {
+  int i = blockDim.x * blockIdx.x + threadIdx.x;
+  int j = blockIdx.y;
+
+  if (i >= ncol)
+    return;
+
+  // we only need to look at the right half as we're using a real2complex FFT
+  int outi = i;
+  //i += ncol-1;
+
+  int pos = __float2int_ru(((j - cen_row) * dv / radius) / du);
+  int pos1 = -pos + cen_col;
+  int pos2 = pos + cen_col;
+
+  if (pos1 > pos2) {
+    int temp = pos1;
+    pos1 = pos2;
+    pos2 = temp;
+    if (pos1 >= ncol) {
+      pos1 = ncol - 1;
+    }
+    if (pos2 < 0) {
+      pos2 = 0;
+    }
+  } else {
+    if (pos1 < 0) {
+      pos1 = 0;
+    }
+    if (pos2 >= ncol) {
+      pos2 = ncol - 1;
+    }
+  }
+
+  float outval = (pos1 <= i && i <= pos2) ? 1.0 : 0.0;
+
+  // mask[cen_row - drop: cen_row + drop + 1, :] = 0
+  if (j >= cen_row - drop && j <= cen_row + drop) {
+    outval = 0;
+  }
+  // mask[:, cen_col - 1: cen_col + 2] = 0
+  if (i >= cen_col - 1 && i <= cen_col + 1) {
+    outval = 0;
+  }
+
+  mask[j * ncol + outi] = outval;
 }