SPECTGPU: run pre-commit hooks on SPECTGPURotateAndGaussianInterpolate

Author: david.roddy

src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu

@@ -6,7 +6,7 @@
   \ingroup projection
   \ingroup SPECTGPU
 
-  \brief implementations for cuda kernel for rotating projector with gaussian interpolation 
+  \brief implementations for cuda kernel for rotating projector with gaussian interpolation
   a la Wallis et al 1997,TMI, doi: 10.1109/42.552061.
 
   \author Daniel Deidda
@@ -26,21 +26,19 @@
 #include <cuda_runtime.h>
 #include <numeric>
 
-//the following is a pull operation
-__global__ void rotateKernel_pull(const float* __restrict__ in_im,
-                             float* __restrict__ out_im,
-                             int3 dim,
-                             float3 spacing,
-                             float3 origin,
-                             float angle_rad)
-                             {
+// the following is a pull operation
+__global__ void
+rotateKernel_pull(
+    const float* __restrict__ in_im, float* __restrict__ out_im, int3 dim, float3 spacing, float3 origin, float angle_rad)
+{
   // parallelise the operation across all image voxels
   int i = threadIdx.x + blockDim.x * blockIdx.x;
   int j = threadIdx.y + blockDim.y * blockIdx.y;
   int k = threadIdx.z + blockDim.z * blockIdx.z;
 
   // check we are not outside the image
-  if (i >= dim.x || j >= dim.y || k >= dim.z) return;
+  if (i >= dim.x || j >= dim.y || k >= dim.z)
+    return;
 
   // get the 1-dimensional index
   int idx = i + (j * dim.x) + (k * dim.x * dim.y);
@@ -68,86 +66,88 @@ __global__ void rotateKernel_pull(const float* __restrict__ in_im,
   // first loop is for finding the value of a gaussian kernel at
   // each nearest neighbour position and summing them all. This will allow normalisation
   // of each NN contribution in the next loop.
-  float G = 0;  // accumulator variable
-  float sigma = 1;  // gaussian kernel sigma
-
-  for (int dr = -1; dr <= 1; dr++)  {
-    int r = (int)roundf(r_r) + dr;  // nearest neighbour z coordinate in rotated voxel space
-    float x_r = (r * spacing.z) - origin.z; // converted to image space
-
-    for (int dq = -1; dq <= 1; dq++)  {
-      int q = (int)roundf(q_r) + dq;
-      float x_q = (q * spacing.y) - origin.y;
-
-      for (int dp = -1; dp <= 1; dp++)  {
-        int p = (int)roundf(p_r) + dp;
-        float x_p = (p * spacing.x) - origin.x;
-
-        // get distance between nearest neighbour and central voxel
-        // both need to be in image space
-        float delta_x = X_rot_x - x_p;
-        float delta_y = X_rot_y - x_q;
-        float delta_z = X_rot_z - x_r;
-
-        // caulculate gaussian kernel 
-        float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma))
-        G += g;
-      };
+  float G = 0;     // accumulator variable
+  float sigma = 1; // gaussian kernel sigma
+
+  for (int dr = -1; dr <= 1; dr++)
+    {
+      int r = (int)roundf(r_r) + dr;          // nearest neighbour z coordinate in rotated voxel space
+      float x_r = (r * spacing.z) - origin.z; // converted to image space
+
+      for (int dq = -1; dq <= 1; dq++)
+        {
+          int q = (int)roundf(q_r) + dq;
+          float x_q = (q * spacing.y) - origin.y;
+
+          for (int dp = -1; dp <= 1; dp++)
+            {
+              int p = (int)roundf(p_r) + dp;
+              float x_p = (p * spacing.x) - origin.x;
+
+              // get distance between nearest neighbour and central voxel
+              // both need to be in image space
+              float delta_x = X_rot_x - x_p;
+              float delta_y = X_rot_y - x_q;
+              float delta_z = X_rot_z - x_r;
+
+              // caulculate gaussian kernel
+              float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)) G
+                  += g;
+            };
+        };
     };
-  };
 
   // loop again but this time actually fetch the image values
   float accumulation = 0;
 
-  for (int dr = -1; dr <= 1; dr++)  {
-    float r = dr + (int)roundf(r_r);
-    float x_r = origin.z + (r * spacing.z);
+  for (int dr = -1; dr <= 1; dr++)
+    {
+      float r = dr + (int)roundf(r_r);
+      float x_r = origin.z + (r * spacing.z);
 
-    for (int dq = -1; dq <= 1; dq++)  {
-      float q = dq + (int)roundf(q_r);
-      float x_q = origin.y + (q * spacing.y);
+      for (int dq = -1; dq <= 1; dq++)
+        {
+          float q = dq + (int)roundf(q_r);
+          float x_q = origin.y + (q * spacing.y);
 
-      for (int dp = -1; dp <= 1; dp++)  {
-        float p = dp + (int)roundf(p_r);
-        float x_p = origin.x + (p * spacing.x);
+          for (int dp = -1; dp <= 1; dp++)
+            {
+              float p = dp + (int)roundf(p_r);
+              float x_p = origin.x + (p * spacing.x);
 
-        float delta_x = Xr_x - x_p;
-        float delta_y = Xr_y - x_q;
-        float delta_z = Xr_z - x_r;
+              float delta_x = Xr_x - x_p;
+              float delta_y = Xr_y - x_q;
+              float delta_z = Xr_z - x_r;
 
-        // get input image voxel index that we are at
-        int i_idx = p + (q * dim.x) + (r * dim.x * dim.y);
+              // get input image voxel index that we are at
+              int i_idx = p + (q * dim.x) + (r * dim.x * dim.y);
 
-        // calculate gaussian weighting for this voxel
-        float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma));
+              // calculate gaussian weighting for this voxel
+              float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma));
 
-        // record the weighted value from this voxel
-        accumulation += in_im[i_idx] * g / G;
-
-      };
+              // record the weighted value from this voxel
+              accumulation += in_im[i_idx] * g / G;
+            };
+        };
     };
-  };
 
   // assign the pulled voxel values to a single voxel in the output image
   out_im[idx] = accumulation;
-
 };
 
 // the following is the adjoint operation (push)
-__global__ void rotateKernel_push(const float* __restrict__ in_im,
-                             float* __restrict__ out_im,
-                             int3 dim,
-                             float3 spacing,
-                             float3 origin,
-                             float angle_rad)
-                             {
-                              // parallelise the operation across all image voxels
+__global__ void
+rotateKernel_push(
+    const float* __restrict__ in_im, float* __restrict__ out_im, int3 dim, float3 spacing, float3 origin, float angle_rad)
+{
+  // parallelise the operation across all image voxels
   int i = threadIdx.x + blockDim.x * blockIdx.x;
   int j = threadIdx.y + blockDim.y * blockIdx.y;
   int k = threadIdx.z + blockDim.z * blockIdx.z;
 
   // check we are not outside the image
-  if (i >= dim.x || j >= dim.y || k >= dim.z) return;
+  if (i >= dim.x || j >= dim.y || k >= dim.z)
+    return;
 
   // get the 1-dimensional index
   int idx = i + (j * dim.x) + (k * dim.x * dim.y);
@@ -175,61 +175,66 @@ __global__ void rotateKernel_push(const float* __restrict__ in_im,
   // first loop is for finding the value of a gaussian kernel at
   // each nearest neighbour position and summing them all. This will allow normalisation
   // of each NN contribution in the next loop.
-  float G = 0;  // accumulator variable
-  float sigma = 1;  // gaussian kernel sigma
-
-  for (int dr = -1; dr <= 1; dr++)  {
-    int r = (int)roundf(r_r) + dr;  // nearest neighbour z coordinate in rotated voxel space
-    float x_r = (r * spacing.z) - origin.z; // converted to image space
-
-    for (int dq = -1; dq <= 1; dq++)  {
-      int q = (int)roundf(q_r) + dq;
-      float x_q = (q * spacing.y) - origin.y;
-
-      for (int dp = -1; dp <= 1; dp++)  {
-        int p = (int)roundf(p_r) + dp;
-        float x_p = (p * spacing.x) - origin.x;
-
-        // get distance between nearest neighbour and central voxel
-        // both need to be in image space
-        float delta_x = X_rot_x - x_p;
-        float delta_y = X_rot_y - x_q;
-        float delta_z = X_rot_z - x_r;
-
-        // caulculate gaussian kernel 
-        float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma))
-        G += g;
-      };
+  float G = 0;     // accumulator variable
+  float sigma = 1; // gaussian kernel sigma
+
+  for (int dr = -1; dr <= 1; dr++)
+    {
+      int r = (int)roundf(r_r) + dr;          // nearest neighbour z coordinate in rotated voxel space
+      float x_r = (r * spacing.z) - origin.z; // converted to image space
+
+      for (int dq = -1; dq <= 1; dq++)
+        {
+          int q = (int)roundf(q_r) + dq;
+          float x_q = (q * spacing.y) - origin.y;
+
+          for (int dp = -1; dp <= 1; dp++)
+            {
+              int p = (int)roundf(p_r) + dp;
+              float x_p = (p * spacing.x) - origin.x;
+
+              // get distance between nearest neighbour and central voxel
+              // both need to be in image space
+              float delta_x = X_rot_x - x_p;
+              float delta_y = X_rot_y - x_q;
+              float delta_z = X_rot_z - x_r;
+
+              // caulculate gaussian kernel
+              float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma)) G
+                  += g;
+            };
+        };
     };
-  };
 
   // loop again but this time actually fetch the image values
   float accumulation = 0;
 
-  for (int dr = -1; dr <= 1; dr++)  {
-    float r = dr + (int)roundf(r_r);
-    float x_r = origin.z + (r * spacing.z);
-
-    for (int dq = -1; dq <= 1; dq++)  {
-      float q = dq + (int)roundf(q_r);
-      float x_q = origin.y + (q * spacing.y);
-
-      for (int dp = -1; dp <= 1; dp++)  {
-        float p = dp + (int)roundf(p_r);
-        float x_p = origin.x + (p * spacing.x);
-
-        float delta_x = Xr_x - x_p;
-        float delta_y = Xr_y - x_q;
-        float delta_z = Xr_z - x_r;
-
-        // get input image voxel index that we are at
-        int i_idx = p + (q * dim.x) + (r * dim.x * dim.y);
-
-        // Pushing the weighted counts to NN
-        float g = expf(-1. * ((delta_x*delta_x) + (delta_y*delta_y) + (delta_z*delta_z)) / (2*sigma*sigma));
-        atomicAdd(&out_im[i_idx], in_im[idx]*g/G); 
-                              
-         };
+  for (int dr = -1; dr <= 1; dr++)
+    {
+      float r = dr + (int)roundf(r_r);
+      float x_r = origin.z + (r * spacing.z);
+
+      for (int dq = -1; dq <= 1; dq++)
+        {
+          float q = dq + (int)roundf(q_r);
+          float x_q = origin.y + (q * spacing.y);
+
+          for (int dp = -1; dp <= 1; dp++)
+            {
+              float p = dp + (int)roundf(p_r);
+              float x_p = origin.x + (p * spacing.x);
+
+              float delta_x = Xr_x - x_p;
+              float delta_y = Xr_y - x_q;
+              float delta_z = Xr_z - x_r;
+
+              // get input image voxel index that we are at
+              int i_idx = p + (q * dim.x) + (r * dim.x * dim.y);
+
+              // Pushing the weighted counts to NN
+              float g = expf(-1. * ((delta_x * delta_x) + (delta_y * delta_y) + (delta_z * delta_z)) / (2 * sigma * sigma));
+              atomicAdd(&out_im[i_idx], in_im[idx] * g / G);
+            };
+        };
     };
-  };
 }