Refactor _reg_maths* and _reg_common_cuda_kernels

niftyreg_build_version.txt

@@ -1 +1 @@
-418
+419

reg-apps/reg_average.cpp

@@ -18,7 +18,6 @@
 #include "_reg_resampling.h"
 #include "_reg_globalTrans.h"
 #include "_reg_localTrans.h"
-#include "_reg_maths_eigen.h"
 
 using PrecisionType = float;
 
@@ -112,7 +111,7 @@ mat44 compute_average_matrices(size_t matrixNumber,
    // Input matrices are logged in place
    for(size_t m=0; m<matrixNumber; ++m)
    {
-      matrices[m] = reg_mat44_logm(&matrices[m]);
+      matrices[m] = Mat44Logm(&matrices[m]);
       matrixWeight[m]=1.;
    }
    // The number of iteration to perform is defined based on the use of lts
@@ -196,7 +195,7 @@ mat44 compute_average_matrices(size_t matrixNumber,
             matrixIndexSorted[m]=m;
          }
          // Sort the computed distances
-         reg_heapSort(matrixWeight, matrixIndexSorted, matrixNumber);
+         HeapSort(matrixWeight, matrixIndexSorted, matrixNumber);
          // Re-assign the weights for the next iteration
          memset(matrixWeight, 0, matrixNumber*sizeof(float));
          for(size_t m=matrixNumber-1; m>lts_inlier * matrixNumber; --m)
@@ -205,7 +204,7 @@ mat44 compute_average_matrices(size_t matrixNumber,
          }
       }
       // The average matrix is exponentiated
-      average_matrix = reg_mat44_expm(&average_matrix);
+      average_matrix = Mat44Expm(&average_matrix);
    } // iteration number
    // Free the allocated array
    free(matrixWeight);
@@ -230,15 +229,15 @@ mat44 compute_affine_demean(size_t matrixNumber,
       tempMatrix=nifti_quatern_to_mat44(qb,qc,qd,qx,qy,qz,1.f,1.f,1.f,qfac);
       // remove the rigid componenent from the affine matrix
       tempMatrix=nifti_mat44_inverse(tempMatrix);
-      tempMatrix=reg_mat44_mul(&tempMatrix,&current_affine);
+      tempMatrix=tempMatrix * current_affine;
       // sum up all the affine matrices
-      tempMatrix = reg_mat44_logm(&tempMatrix);
+      tempMatrix = Mat44Logm(&tempMatrix);
       demeanMatrix = demeanMatrix + tempMatrix;
    }
    // The average matrix is normalised
-   demeanMatrix = reg_mat44_mul(&demeanMatrix,1.f/(float)matrixNumber);
+   demeanMatrix = demeanMatrix * (1.f / (float)matrixNumber);
    // The average matrix is exponentiated
-   demeanMatrix = reg_mat44_expm(&demeanMatrix);
+   demeanMatrix = Mat44Expm(&demeanMatrix);
    // The average matrix is inverted
    demeanMatrix = nifti_mat44_inverse(demeanMatrix);
    return demeanMatrix;
@@ -293,11 +292,7 @@ int compute_nrr_demean(nifti_image *demean_field,
             affineTransformation=*reinterpret_cast<mat44 *>(transformation->ext_list[0].edata);
             // Note that if the transformation is a flow field, only half-of the affine has be used
             if(transformation->num_ext>1 && deformationField->intent_p1!=DEF_VEL_FIELD)
-            {
-               affineTransformation=reg_mat44_mul(
-                                       reinterpret_cast<mat44 *>(transformation->ext_list[1].edata),
-                                       &affineTransformation);
-            }
+               affineTransformation=reinterpret_cast<mat44&>(*transformation->ext_list[1].edata) * affineTransformation;
          }
          else reg_tool_ReadAffineFile(&affineTransformation,inputAffName[t]);
          // The affine component is substracted

reg-apps/reg_benchmark.cpp

@@ -153,7 +153,7 @@ int main(int argc, char **argv)
    originIndex[0] = -1.0f;
    originIndex[1] = -1.0f;
    originIndex[2] = -1.0f;
-   reg_mat44_mul(&(controlPointImage->qto_xyz), originIndex, originReal);
+   Mat44Mul(controlPointImage->qto_xyz, originIndex, originReal);
    controlPointImage->qto_xyz.m[0][3] = controlPointImage->qoffset_x = originReal[0];
    controlPointImage->qto_xyz.m[1][3] = controlPointImage->qoffset_y = originReal[1];
    controlPointImage->qto_xyz.m[2][3] = controlPointImage->qoffset_z = originReal[2];

reg-apps/reg_gpuinfo.cpp

@@ -1,4 +1,4 @@
-#include "_reg_maths.h"
+#include "Maths.hpp"
 #include "Platform.h"
 
 #ifdef USE_CUDA

reg-apps/reg_jacobian.cpp

@@ -233,7 +233,7 @@ int main(int argc, char **argv)
       if(!reg_isAnImageFileName(param->inputTransName)){
          mat44 *affineTransformation=(mat44 *)malloc(sizeof(mat44));
          reg_tool_ReadAffineFile(affineTransformation,param->inputTransName);
-         NR_COUT << reg_mat44_det<double>(affineTransformation) << std::endl;
+         NR_COUT << Mat44Det<double>(affineTransformation) << std::endl;
          return EXIT_SUCCESS;
       }
 

reg-apps/reg_resample.cpp

@@ -307,7 +307,7 @@ int main(int argc, char **argv)
    else
    {
       // No transformation is specified, an identity transformation is used
-      reg_mat44_eye(&inputAffineTransformation);
+      Mat44Eye(&inputAffineTransformation);
    }
 
    // Create a deformation field

reg-apps/reg_tools.cpp

@@ -498,7 +498,7 @@ int main(int argc, char **argv)
     {
         reg_tools_changeDatatype<float>(image);
         nifti_image *normImage = nifti_dup(*image);
-        reg_heapSort(static_cast<float *>(normImage->data), normImage->nvox);
+        HeapSort(static_cast<float *>(normImage->data), normImage->nvox);
         float minValue = static_cast<float *>(normImage->data)[Floor(03*(int)normImage->nvox/100)];
         float maxValue = static_cast<float *>(normImage->data)[Floor(97*(int)normImage->nvox/100)];
         reg_tools_subtractValueFromImage(image,normImage,minValue);
@@ -892,7 +892,7 @@ int main(int argc, char **argv)
         const size_t jacobianVoxelNumber = NiftiImage::calcVoxelNumber(def, 3);
         mat33 *jacobian = (mat33 *)malloc(jacobianVoxelNumber * sizeof(mat33));
         for (size_t i = 0; i < jacobianVoxelNumber; ++i)
-            reg_mat33_eye(&jacobian[i]);
+            Mat33Eye(&jacobian[i]);
         // resample the original image into the space of the new image
         if(flag->interpFlag == 0){
             param->interpOrder = 3;

reg-apps/reg_transform.cpp

@@ -16,7 +16,6 @@
 #include "_reg_globalTrans.h"
 #include "_reg_localTrans.h"
 #include "_reg_tools.h"
-#include "_reg_maths_eigen.h"
 
 #include "reg_transform.h"
 
@@ -532,7 +531,7 @@ int main(int argc, char **argv) {
         if (affine1Trans != nullptr && affine2Trans != nullptr) {
             NR_INFO("Transformation 2 is an affine parametrisation:");
             NR_INFO(param->input2TransName);
-            *affine1Trans = reg_mat44_mul(affine2Trans, affine1Trans);
+            *affine1Trans = *affine2Trans * *affine1Trans;
             reg_tool_WriteAffineFile(affine1Trans, param->outputTransName);
         } else {
             // Check if the reference image is required
@@ -955,10 +954,10 @@ int main(int argc, char **argv) {
 
         // Update the sform
         if (image->sform_code > 0) {
-            image->sto_xyz = reg_mat44_mul(affineTransformation, &(image->sto_xyz));
+            image->sto_xyz = *affineTransformation * image->sto_xyz;
         } else {
             image->sform_code = 1;
-            image->sto_xyz = reg_mat44_mul(affineTransformation, &(image->qto_xyz));
+            image->sto_xyz = *affineTransformation * image->qto_xyz;
         }
         image->sto_ijk = nifti_mat44_inverse(image->sto_xyz);
 
@@ -980,9 +979,9 @@ int main(int argc, char **argv) {
             affineTrans = (mat44 *)malloc(sizeof(mat44));
             reg_tool_ReadAffineFile(affineTrans, param->inputTransName);
             // The affine transformation is halfed
-            *affineTrans = reg_mat44_logm(affineTrans);
-            *affineTrans = reg_mat44_mul(affineTrans, 0.5);
-            *affineTrans = reg_mat44_expm(affineTrans);
+            *affineTrans = Mat44Logm(affineTrans);
+            *affineTrans = *affineTrans * 0.5;
+            *affineTrans = Mat44Expm(affineTrans);
             // The affine transformation is saved
             reg_tool_WriteAffineFile(affineTrans, param->outputTransName);
         } else {
@@ -1183,17 +1182,17 @@ int main(int argc, char **argv) {
     if (flag->makeAffFlag) {
         // Create all the required matrices
         mat44 rotationX;
-        reg_mat44_eye(&rotationX);
+        Mat44Eye(&rotationX);
         mat44 translation;
-        reg_mat44_eye(&translation);
+        Mat44Eye(&translation);
         mat44 rotationY;
-        reg_mat44_eye(&rotationY);
+        Mat44Eye(&rotationY);
         mat44 rotationZ;
-        reg_mat44_eye(&rotationZ);
+        Mat44Eye(&rotationZ);
         mat44 scaling;
-        reg_mat44_eye(&scaling);
+        Mat44Eye(&scaling);
         mat44 shearing;
-        reg_mat44_eye(&shearing);
+        Mat44Eye(&shearing);
         // Set up the rotation matrix along the YZ plane
         rotationX.m[1][1] = cosf(param->affTransParam[0]);
         rotationX.m[1][2] = -sinf(param->affTransParam[0]);
@@ -1222,11 +1221,11 @@ int main(int argc, char **argv) {
         shearing.m[2][0] = param->affTransParam[10];
         shearing.m[2][1] = param->affTransParam[11];
         // Combine all the transformations
-        mat44 affine = reg_mat44_mul(&rotationY, &rotationZ);
-        affine = reg_mat44_mul(&rotationX, &affine);
-        affine = reg_mat44_mul(&scaling, &affine);
-        affine = reg_mat44_mul(&shearing, &affine);
-        affine = reg_mat44_mul(&translation, &affine);
+        mat44 affine = rotationY * rotationZ;
+        affine = rotationX * affine;
+        affine = scaling * affine;
+        affine = shearing * affine;
+        affine = translation * affine;
         // Save the new matrix
         reg_tool_WriteAffineFile(&affine, param->outputTransName);
     }

reg-io/_reg_ReadWriteMatrix.cpp

@@ -70,11 +70,11 @@ void reg_tool_ReadAffineFile(mat44 *mat,
         absoluteFloating = nifti_mat44_inverse(absoluteFloating);
         *mat = nifti_mat44_inverse(*mat);
 
-        *mat = reg_mat44_mul(&absoluteFloating, mat);
-        *mat = reg_mat44_mul(mat, &absoluteReference);
-        *mat = reg_mat44_mul(floatingMatrix, mat);
+        *mat = absoluteFloating * *mat;
+        *mat = *mat * absoluteReference;
+        *mat = *floatingMatrix * *mat;
         mat44 tmp = nifti_mat44_inverse(*referenceMatrix);
-        *mat = reg_mat44_mul(mat, &tmp);
+        *mat = *mat * tmp;
     }
 
     NR_MAT44_DEBUG(*mat, "Affine matrix");
@@ -168,7 +168,7 @@ template<class T>
 T** reg_tool_ReadMatrixFile(char *filename, size_t nbLine, size_t nbColumn) {
     //THEN CONSTRUCT THE MATRIX
     // Allocate the matrices
-    T** mat = reg_matrix2DAllocate<T>(nbLine, nbColumn);
+    T** mat = Matrix2dAlloc<T>(nbLine, nbColumn);
     //STORE THE VALUES
     std::string line;
     std::ifstream matrixFile(filename);

reg-io/nrrd/reg_nrrd.cpp

@@ -167,7 +167,7 @@ nifti_image *reg_io_nrdd2nifti(Nrrd *nrrdImage)
 
    // The space orientation is extracted and converted into a matrix
    mat44 qform_orientation_matrix;
-   reg_mat44_eye(&qform_orientation_matrix);
+   Mat44Eye(&qform_orientation_matrix);
    if(nrrdImage->space==nrrdSpaceRightAnteriorSuperior ||
          nrrdImage->space==nrrdSpaceRightAnteriorSuperiorTime ||
          nrrdImage->space==nrrdSpace3DRightHanded ||
@@ -251,7 +251,7 @@ nifti_image *reg_io_nrdd2nifti(Nrrd *nrrdImage)
    if(nrrdImage->axis[1].spaceDirection[0]!=std::numeric_limits<double>::quiet_NaN())
    {
       niiImage->sform_code=1;
-      reg_mat44_eye(&niiImage->sto_xyz);
+      Mat44Eye(&niiImage->sto_xyz);
       for(int i=0; i<(niiImage->ndim<3?niiImage->ndim:3); ++i)
       {
          for(int j=0; j<(niiImage->ndim<3?niiImage->ndim:3); ++j)

reg-lib/CMakeLists.txt

@@ -12,8 +12,7 @@ endif(USE_OPENCL)
 ##BUILD THE CPU LIBRARIES
 #-----------------------------------------------------------------------------
 add_library(_reg_maths ${NIFTYREG_LIBRARY_TYPE}
-  cpu/_reg_maths.cpp
-  cpu/_reg_maths_eigen.cpp
+  cpu/Maths.cpp
 )
 install(TARGETS _reg_maths
   RUNTIME DESTINATION bin

reg-lib/Debug.hpp

@@ -68,11 +68,11 @@ inline std::string StripFunctionName(const std::string& funcName) {
 #define NR_INFO(msg)        NR_COUT << "[NiftyReg INFO] " << msg << std::endl
 /* *************************************************************** */
 #ifndef NDEBUG
-#define NR_MAT44(mat, title)          reg_mat44_disp(mat, "[NiftyReg DEBUG] "s + (title))
+#define NR_MAT44(mat, title)          Mat44Disp(mat, "[NiftyReg DEBUG] "s + (title))
 #define NR_MAT44_DEBUG(mat, title)    NR_MAT44(mat, title)
 #define NR_MAT44_VERBOSE(mat, title)  NR_MAT44(mat, title)
 #else
-#define NR_MAT44(mat, title)          reg_mat44_disp(mat, title)
+#define NR_MAT44(mat, title)          Mat44Disp(mat, title)
 #define NR_MAT44_DEBUG(mat, title)
 #define NR_MAT44_VERBOSE(mat, title)  if (this->verbose) NR_MAT44(mat, "[NiftyReg INFO] "s + (title))
 #endif

reg-lib/_reg_aladin.cpp

@@ -177,9 +177,9 @@ void reg_aladin<T>::InitialiseRegistration() {
             referenceCenter[2] = (float)(this->inputReference->nz) / 2.0f;
             //From pixel coordinates to real coordinates
             float floatingRealPosition[3];
-            reg_mat44_mul(floatingMatrix, floatingCenter, floatingRealPosition);
+            Mat44Mul(*floatingMatrix, floatingCenter, floatingRealPosition);
             float referenceRealPosition[3];
-            reg_mat44_mul(referenceMatrix, referenceCenter, referenceRealPosition);
+            Mat44Mul(*referenceMatrix, referenceCenter, referenceRealPosition);
             //Set translation to the transformation matrix
             this->affineTransformation->m[0][3] = floatingRealPosition[0] - referenceRealPosition[0];
             this->affineTransformation->m[1][3] = floatingRealPosition[1] - referenceRealPosition[1];
@@ -207,7 +207,7 @@ void reg_aladin<T>::InitialiseRegistration() {
             referenceCentre[2] /= referenceCount;
             float refCOM[3]{};
             if (this->inputReference->sform_code > 0)
-                reg_mat44_mul(&this->inputReference->sto_xyz, referenceCentre, refCOM);
+                Mat44Mul(this->inputReference->sto_xyz, referenceCentre, refCOM);
 
             float floatingCentre[3] = { 0, 0, 0 };
             float floatingCount = 0;
@@ -231,8 +231,8 @@ void reg_aladin<T>::InitialiseRegistration() {
             floatingCentre[2] /= floatingCount;
             float floCOM[3]{};
             if (this->inputFloating->sform_code > 0)
-                reg_mat44_mul(&this->inputFloating->sto_xyz, floatingCentre, floCOM);
-            reg_mat44_eye(this->affineTransformation.get());
+                Mat44Mul(this->inputFloating->sto_xyz, floatingCentre, floCOM);
+            Mat44Eye(this->affineTransformation.get());
             this->affineTransformation->m[0][3] = floCOM[0] - refCOM[0];
             this->affineTransformation->m[1][3] = floCOM[1] - refCOM[1];
             this->affineTransformation->m[2][3] = floCOM[2] - refCOM[2];

reg-lib/_reg_aladin_sym.cpp

@@ -1,5 +1,4 @@
 #include "_reg_aladin_sym.h"
-#include "_reg_maths_eigen.h"
 
 /* *************************************************************** */
 template <class T>
@@ -81,7 +80,7 @@ void reg_aladin_sym<T>::InitialiseRegistration() {
         referenceCentre[2] /= referenceCount;
         float refCOG[3]{};
         if (this->inputReference->sform_code > 0)
-            reg_mat44_mul(&(this->inputReference->sto_xyz), referenceCentre, refCOG);
+            Mat44Mul(this->inputReference->sto_xyz, referenceCentre, refCOG);
 
         float floatingCentre[3] = { 0, 0, 0 };
         float floatingCount = 0;
@@ -106,8 +105,8 @@ void reg_aladin_sym<T>::InitialiseRegistration() {
         floatingCentre[2] /= floatingCount;
         float floCOG[3]{};
         if (this->inputFloating->sform_code > 0)
-            reg_mat44_mul(&(this->inputFloating->sto_xyz), floatingCentre, floCOG);
-        reg_mat44_eye(this->affineTransformation.get());
+            Mat44Mul(this->inputFloating->sto_xyz, floatingCentre, floCOG);
+        Mat44Eye(this->affineTransformation.get());
         this->affineTransformation->m[0][3] = floCOG[0] - refCOG[0];
         this->affineTransformation->m[1][3] = floCOG[1] - refCOG[1];
         this->affineTransformation->m[2][3] = floCOG[2] - refCOG[2];
@@ -143,9 +142,9 @@ void reg_aladin_sym<T>::UpdateTransformationMatrix(int type) {
     mat44 bInverted = nifti_mat44_inverse(*this->affineTransformationBw);
 
     // We average the forward and inverted backward matrix
-    *this->affineTransformation = reg_mat44_avg2(this->affineTransformation.get(), &bInverted);
+    *this->affineTransformation = Mat44Avg2(this->affineTransformation.get(), &bInverted);
     // We average the inverted forward and backward matrix
-    *this->affineTransformationBw = reg_mat44_avg2(&fInverted, this->affineTransformationBw.get());
+    *this->affineTransformationBw = Mat44Avg2(&fInverted, this->affineTransformationBw.get());
     for (int i = 0; i < 3; ++i) {
         this->affineTransformation->m[3][i] = 0.f;
         this->affineTransformationBw->m[3][i] = 0.f;

reg-lib/_reg_base.cpp

@@ -402,7 +402,7 @@ void reg_base<T>::Initialise() {
         reg_tools_changeDatatype<T>(tmpReference);
         // Extract the robust range of the reference image
         T *refDataPtr = static_cast<T *>(tmpReference->data);
-        reg_heapSort(refDataPtr, tmpReference->nvox);
+        HeapSort(refDataPtr, tmpReference->nvox);
         // Update the reference threshold values if no value has been setup by the user
         if (referenceThresholdLow[0] == std::numeric_limits<T>::lowest())
             referenceThresholdLow[0] = refDataPtr[Round((float)tmpReference->nvox * 0.02f)];
@@ -414,7 +414,7 @@ void reg_base<T>::Initialise() {
         reg_tools_changeDatatype<T>(tmpFloating);
         // Extract the robust range of the floating image
         T *floDataPtr = static_cast<T *>(tmpFloating->data);
-        reg_heapSort(floDataPtr, tmpFloating->nvox);
+        HeapSort(floDataPtr, tmpFloating->nvox);
         // Update the floating threshold values if no value has been setup by the user
         if (floatingThresholdLow[0] == std::numeric_limits<T>::lowest())
             floatingThresholdLow[0] = floDataPtr[Round((float)tmpFloating->nvox * 0.02f)];

reg-lib/cl/ClAffineDeformationFieldKernel.cpp

@@ -89,7 +89,7 @@ void ClAffineDeformationFieldKernel::Calculate(bool compose) {
     const size_t globalWorkSize[dims] = {xBlocks * xThreads, yBlocks * yThreads, zBlocks * zThreads};
     const size_t localWorkSize[dims] = {xThreads, yThreads, zThreads};
 
-    mat44 transformationMatrix = compose ? *affineTransformation : reg_mat44_mul(affineTransformation, referenceMatrix);
+    mat44 transformationMatrix = compose ? *affineTransformation : *affineTransformation * *referenceMatrix;
 
     float* trans = (float *)malloc(16 * sizeof(float));
     mat44ToCptr(transformationMatrix, trans);