Refactorisations

KCL-BMEIS · Oct 5, 2023 · e1ec1f4 · e1ec1f4
1 parent 118e1da
commit e1ec1f4
Show file tree

Hide file tree

Showing 9 changed files with 142 additions and 133 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -71,7 +71,7 @@ option(BUILD_TESTING "To build the unit tests" OFF)
 option(USE_CUDA "To use the CUDA platform" OFF)
 option(USE_OPENCL "To use the OpenCL platform" OFF)
 option(USE_OPENMP "To use openMP for multi-CPU processing" ON)
-option(USE_SSE "To enable SEE computation in some case" ON)
+option(USE_SSE "To enable SSE computation in some case" ON)
 #-----------------------------------------------------------------------------
 option(USE_NRRD "To use the NRRD file format" OFF)
 mark_as_advanced(USE_NRRD)

diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt
@@ -1 +1 @@
-334
+335
diff --git a/reg-lib/_reg_f3d.cpp b/reg-lib/_reg_f3d.cpp
@@ -482,6 +482,7 @@ void reg_f3d<T>::SetOptimiser() {
 template<class T>
 void reg_f3d<T>::SmoothGradient() {
     // The gradient is smoothed using a Gaussian kernel if it is required
+    if (this->gradientSmoothingSigma == 0) return;
     this->compute->SmoothGradient(this->gradientSmoothingSigma);
     NR_FUNC_CALLED();
 }

diff --git a/reg-lib/_reg_f3d2.cpp b/reg-lib/_reg_f3d2.cpp
@@ -348,11 +348,10 @@ void reg_f3d2<T>::GetLandmarkDistanceGradient() {
 /* *************************************************************** */
 template <class T>
 void reg_f3d2<T>::SmoothGradient() {
-    reg_f3d<T>::SmoothGradient();
-
     // The gradient is smoothed using a Gaussian kernel if it is required
+    if (this->gradientSmoothingSigma == 0) return;
+    reg_f3d<T>::SmoothGradient();
     computeBw->SmoothGradient(this->gradientSmoothingSigma);
-
     NR_FUNC_CALLED();
 }
 /* *************************************************************** */

diff --git a/reg-lib/cpu/_reg_maths.h b/reg-lib/cpu/_reg_maths.h
@@ -56,6 +56,10 @@ DEVICE inline T Square(const T& x) {
     return x * x;
 }
 template<typename T>
+DEVICE inline T Cube(const T& x) {
+    return x * x * x;
+}
+template<typename T>
 DEVICE inline int Floor(const T& x) {
     const int i = static_cast<int>(x);
     return i - (x < i);

diff --git a/reg-lib/cpu/_reg_tools.cpp b/reg-lib/cpu/_reg_tools.cpp
diff --git a/reg-lib/cpu/_reg_tools.h b/reg-lib/cpu/_reg_tools.h
@@ -81,10 +81,14 @@ void reg_tools_removeSCLInfo(nifti_image *img);
 void reg_getRealImageSpacing(nifti_image *image,
                              float *spacingValues);
 /* *************************************************************** */
-/** @brief Smooth an image using a Gaussian kernel
+/** @brief Smooth an image using a specified kernel
  * @param image Image to be smoothed
- * @param sigma Standard deviation of the Gaussian kernel
- * to use. The kernel is bounded between +/- 3 sigma.
+ * @param sigma Standard deviation of the kernel to use.
+ * The kernel is bounded between +/- 3 sigma.
+ * @param kernelType Type of kernel to use.
+ * @param mask An integer mask over which the smoothing should occur.
+ * @param timePoints Boolean array to specify which time points have to be
+ * smoothed. The array follow the dim array of the nifti header.
  * @param axis Boolean array to specify which axis have to be
  * smoothed. The array follow the dim array of the nifti header.
  */
@@ -100,16 +104,16 @@ void reg_tools_kernelConvolution(nifti_image *image,
  * @param varianceX The variance of the Gaussian kernel in X
  * @param varianceY The variance of the Gaussian kernel in Y
  * @param varianceZ The variance of the Gaussian kernel in Z
- * @param mask An integer mask over which the Gaussian smoothing should occur
- * @param timePoint Boolean array to specify which timepoints have to be
+ * @param mask An integer mask over which the Gaussian smoothing should occur.
+ * @param timePoints Boolean array to specify which time points have to be
  * smoothed.
  */
 void reg_tools_labelKernelConvolution(nifti_image *image,
                                       float varianceX,
                                       float varianceY,
                                       float varianceZ,
                                       int *mask = nullptr,
-                                      bool *timePoint = nullptr);
+                                      bool *timePoints = nullptr);
 /* *************************************************************** */
 /** @brief Downsample an image by a ratio of two
  * @param image Image to be downsampled

diff --git a/reg-lib/cuda/CudaCommon.cu b/reg-lib/cuda/CudaCommon.cu
@@ -82,22 +82,23 @@ void TransferNiftiToDevice(cudaArray *arrayCuda, const nifti_image *img) {
 template <class DataType>
 void TransferNiftiToDevice(cudaArray *arrayCuda, const nifti_image *img) {
     if (sizeof(DataType) == sizeof(float4)) {
-        if (img->datatype != NIFTI_TYPE_FLOAT32 || img->dim[5] < 2 || img->dim[4] > 1)
-            NR_FATAL_ERROR("The specified image is not a single precision deformation field image");
+        if (img->datatype != NIFTI_TYPE_FLOAT32)
+            NR_FATAL_ERROR("The specified image is not a single precision image");
         const float *niftiImgValues = static_cast<float*>(img->data);
         const size_t voxelNumber = NiftiImage::calcVoxelNumber(img, 3);
+        const auto timePointCount = img->dim[4] * img->dim[5];
         unique_ptr<float4[]> array(new float4[voxelNumber]());
         for (size_t i = 0; i < voxelNumber; i++)
             array[i].x = *niftiImgValues++;
-        if (img->dim[5] >= 2) {
+        if (timePointCount >= 2) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array[i].y = *niftiImgValues++;
         }
-        if (img->dim[5] >= 3) {
+        if (timePointCount >= 3) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array[i].z = *niftiImgValues++;
         }
-        if (img->dim[5] >= 4) {
+        if (timePointCount >= 4) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array[i].w = *niftiImgValues++;
         }
@@ -153,29 +154,30 @@ void TransferNiftiToDevice(cudaArray *array1Cuda, cudaArray *array2Cuda, const n
 template <class DataType>
 void TransferNiftiToDevice(cudaArray *array1Cuda, cudaArray *array2Cuda, const nifti_image *img) {
     if (sizeof(DataType) == sizeof(float4)) {
-        if (img->datatype != NIFTI_TYPE_FLOAT32 || img->dim[5] < 2 || img->dim[4] > 1)
-            NR_FATAL_ERROR("The specified image is not a single precision deformation field image");
+        if (img->datatype != NIFTI_TYPE_FLOAT32)
+            NR_FATAL_ERROR("The specified image is not a single precision image");
         const float *niftiImgValues = static_cast<float*>(img->data);
         const size_t voxelNumber = NiftiImage::calcVoxelNumber(img, 3);
+        const auto timePointCount = img->dim[4] * img->dim[5];
         unique_ptr<float4[]> array1(new float4[voxelNumber]());
         unique_ptr<float4[]> array2(new float4[voxelNumber]());
         for (size_t i = 0; i < voxelNumber; i++)
             array1[i].x = *niftiImgValues++;
         for (size_t i = 0; i < voxelNumber; i++)
             array2[i].x = *niftiImgValues++;
-        if (img->dim[5] >= 2) {
+        if (timePointCount >= 2) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array1[i].y = *niftiImgValues++;
             for (size_t i = 0; i < voxelNumber; i++)
                 array2[i].y = *niftiImgValues++;
         }
-        if (img->dim[5] >= 3) {
+        if (timePointCount >= 3) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array1[i].z = *niftiImgValues++;
             for (size_t i = 0; i < voxelNumber; i++)
                 array2[i].z = *niftiImgValues++;
         }
-        if (img->dim[5] >= 4) {
+        if (timePointCount >= 4) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array1[i].w = *niftiImgValues++;
             for (size_t i = 0; i < voxelNumber; i++)
@@ -223,22 +225,23 @@ void TransferNiftiToDevice(DataType *arrayCuda, const nifti_image *img) {
 template <class DataType>
 void TransferNiftiToDevice(DataType *arrayCuda, const nifti_image *img) {
     if (sizeof(DataType) == sizeof(float4)) {
-        if (img->datatype != NIFTI_TYPE_FLOAT32 || img->dim[5] < 2 || img->dim[4] > 1)
-            NR_FATAL_ERROR("The specified image is not a single precision deformation field image");
+        if (img->datatype != NIFTI_TYPE_FLOAT32)
+            NR_FATAL_ERROR("The specified image is not a single precision image");
         const float *niftiImgValues = static_cast<float*>(img->data);
         const size_t voxelNumber = NiftiImage::calcVoxelNumber(img, 3);
+        const auto timePointCount = img->dim[4] * img->dim[5];
         unique_ptr<float4[]> array(new float4[voxelNumber]());
         for (size_t i = 0; i < voxelNumber; i++)
             array[i].x = *niftiImgValues++;
-        if (img->dim[5] >= 2) {
+        if (timePointCount >= 2) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array[i].y = *niftiImgValues++;
         }
-        if (img->dim[5] >= 3) {
+        if (timePointCount >= 3) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array[i].z = *niftiImgValues++;
         }
-        if (img->dim[5] >= 4) {
+        if (timePointCount >= 4) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array[i].w = *niftiImgValues++;
         }
@@ -273,29 +276,30 @@ void TransferNiftiToDevice(DataType *array1Cuda, DataType *array2Cuda, const nif
 template <class DataType>
 void TransferNiftiToDevice(DataType *array1Cuda, DataType *array2Cuda, const nifti_image *img) {
     if (sizeof(DataType) == sizeof(float4)) {
-        if (img->datatype != NIFTI_TYPE_FLOAT32 || img->dim[5] < 2 || img->dim[4] > 1)
-            NR_FATAL_ERROR("The specified image is not a single precision deformation field image");
+        if (img->datatype != NIFTI_TYPE_FLOAT32)
+            NR_FATAL_ERROR("The specified image is not a single precision image");
         const float *niftiImgValues = static_cast<float*>(img->data);
         const size_t voxelNumber = NiftiImage::calcVoxelNumber(img, 3);
+        const auto timePointCount = img->dim[4] * img->dim[5];
         unique_ptr<float4[]> array1(new float4[voxelNumber]());
         unique_ptr<float4[]> array2(new float4[voxelNumber]());
         for (size_t i = 0; i < voxelNumber; i++)
             array1[i].x = *niftiImgValues++;
         for (size_t i = 0; i < voxelNumber; i++)
             array2[i].x = *niftiImgValues++;
-        if (img->dim[5] >= 2) {
+        if (timePointCount >= 2) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array1[i].y = *niftiImgValues++;
             for (size_t i = 0; i < voxelNumber; i++)
                 array2[i].y = *niftiImgValues++;
         }
-        if (img->dim[5] >= 3) {
+        if (timePointCount >= 3) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array1[i].z = *niftiImgValues++;
             for (size_t i = 0; i < voxelNumber; i++)
                 array2[i].z = *niftiImgValues++;
         }
-        if (img->dim[5] >= 4) {
+        if (timePointCount >= 4) {
             for (size_t i = 0; i < voxelNumber; i++)
                 array1[i].w = *niftiImgValues++;
             for (size_t i = 0; i < voxelNumber; i++)
@@ -350,23 +354,24 @@ template <class DataType>
 void TransferFromDeviceToNifti(nifti_image *img, const DataType *arrayCuda) {
     if (sizeof(DataType) == sizeof(float4)) {
         // A nifti 5D volume is expected
-        if (img->dim[0] < 5 || img->dim[4]>1 || img->dim[5] < 2 || img->datatype != NIFTI_TYPE_FLOAT32)
-            NR_FATAL_ERROR("The nifti image is not a 5D volume");
+        if (img->datatype != NIFTI_TYPE_FLOAT32)
+            NR_FATAL_ERROR("The specified image is not a single precision image");
         const size_t voxelNumber = NiftiImage::calcVoxelNumber(img, 3);
+        const auto timePointCount = img->dim[4] * img->dim[5];
         thrust::device_ptr<const float4> arrayCudaPtr(reinterpret_cast<const float4*>(arrayCuda));
         const thrust::host_vector<float4> array(arrayCudaPtr, arrayCudaPtr + voxelNumber);
         float *niftiImgValues = static_cast<float*>(img->data);
         for (size_t i = 0; i < voxelNumber; i++)
             *niftiImgValues++ = array[i].x;
-        if (img->dim[5] >= 2) {
+        if (timePointCount >= 2) {
             for (size_t i = 0; i < voxelNumber; i++)
                 *niftiImgValues++ = array[i].y;
         }
-        if (img->dim[5] >= 3) {
+        if (timePointCount >= 3) {
             for (size_t i = 0; i < voxelNumber; i++)
                 *niftiImgValues++ = array[i].z;
         }
-        if (img->dim[5] >= 4) {
+        if (timePointCount >= 4) {
             for (size_t i = 0; i < voxelNumber; i++)
                 *niftiImgValues++ = array[i].w;
         }
@@ -399,9 +404,10 @@ template <class DataType>
 void TransferFromDeviceToNifti(nifti_image *img, const DataType *array1Cuda, const DataType *array2Cuda) {
     if (sizeof(DataType) == sizeof(float4)) {
         // A nifti 5D volume is expected
-        if (img->dim[0] < 5 || img->dim[4]>1 || img->dim[5] < 2 || img->datatype != NIFTI_TYPE_FLOAT32)
-            NR_FATAL_ERROR("The nifti image is not a 5D volume");
+        if (img->datatype != NIFTI_TYPE_FLOAT32)
+            NR_FATAL_ERROR("The specified image is not a single precision image");
         const size_t voxelNumber = NiftiImage::calcVoxelNumber(img, 3);
+        const auto timePointCount = img->dim[4] * img->dim[5];
         thrust::device_ptr<const float4> array1CudaPtr(reinterpret_cast<const float4*>(array1Cuda));
         thrust::device_ptr<const float4> array2CudaPtr(reinterpret_cast<const float4*>(array2Cuda));
         const thrust::host_vector<float4> array1(array1CudaPtr, array1CudaPtr + voxelNumber);
@@ -411,19 +417,19 @@ void TransferFromDeviceToNifti(nifti_image *img, const DataType *array1Cuda, con
             *niftiImgValues++ = array1[i].x;
         for (size_t i = 0; i < voxelNumber; i++)
             *niftiImgValues++ = array2[i].x;
-        if (img->dim[5] >= 2) {
+        if (timePointCount >= 2) {
             for (size_t i = 0; i < voxelNumber; i++)
                 *niftiImgValues++ = array1[i].y;
             for (size_t i = 0; i < voxelNumber; i++)
                 *niftiImgValues++ = array2[i].y;
         }
-        if (img->dim[5] >= 3) {
+        if (timePointCount >= 3) {
             for (size_t i = 0; i < voxelNumber; i++)
                 *niftiImgValues++ = array1[i].z;
             for (size_t i = 0; i < voxelNumber; i++)
                 *niftiImgValues++ = array2[i].z;
         }
-        if (img->dim[5] >= 4) {
+        if (timePointCount >= 4) {
             for (size_t i = 0; i < voxelNumber; i++)
                 *niftiImgValues++ = array1[i].w;
             for (size_t i = 0; i < voxelNumber; i++)

diff --git a/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp b/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp
@@ -7,15 +7,15 @@
  *  to ensure the CPU and CUDA versions yield the same output
 **/
 
-class ApproxLinearEnergyGradient {
+class ApproxLinearEnergyGradientTest {
 protected:
     using TestData = std::tuple<std::string, NiftiImage&, NiftiImage&, NiftiImage&, float>;
     using TestCase = std::tuple<std::string, double, double, NiftiImage, NiftiImage>;
 
     inline static vector<TestCase> testCases;
 
 public:
-    ApproxLinearEnergyGradient() {
+    ApproxLinearEnergyGradientTest() {
         if (!testCases.empty())
             return;
 
@@ -124,7 +124,7 @@ class ApproxLinearEnergyGradient {
     }
 };
 
-TEST_CASE_METHOD(ApproxLinearEnergyGradient, "Regression Approximate Linear Energy Gradient", "[regression]") {
+TEST_CASE_METHOD(ApproxLinearEnergyGradientTest, "Regression Approximate Linear Energy Gradient", "[regression]") {
     // Loop over all generated test cases
     for (auto&& testCase : testCases) {
         // Retrieve test information