#92: fix issue in reg_resample - fix failing test - fix output display

.gitignore

@@ -42,6 +42,7 @@ CMakeSettings.json
 
 # Build
 build*
+out*
 
 # Doxygen
 html

niftyreg_build_version.txt

@@ -1 +1 @@
-315
+316

reg-apps/reg_average.cpp

@@ -382,23 +382,8 @@ int compute_average_image(nifti_image *averageImage,
    // Loop over all input images
    for(size_t i=0; i<imageNumber; ++i){
       // Generate a deformation field defined by the average final
-      nifti_image *deformationField=nifti_copy_nim_info(averageImage);
-      deformationField->ndim=deformationField->dim[0]=5;
-      deformationField->nt=deformationField->dim[4]=1;
-      deformationField->nu=deformationField->dim[5]=deformationField->nz>1?3:2;
-      deformationField->nvox=NiftiImage::calcVoxelNumber(deformationField, deformationField->ndim);
-      deformationField->nbyper=sizeof(float);
-      deformationField->datatype=NIFTI_TYPE_FLOAT32;
-      deformationField->intent_code=NIFTI_INTENT_VECTOR;
-      memset(deformationField->intent_name, 0, 16);
-      strcpy(deformationField->intent_name,"NREG_TRANS");
-      deformationField->scl_slope=1.f;
-      deformationField->scl_inter=0.f;
-      deformationField->intent_p1=DISP_FIELD;
-      deformationField->data=calloc(deformationField->nvox, deformationField->nbyper);
-      reg_tools_multiplyValueToImage(deformationField,deformationField,0.f);
-      // Set the transformation to identity
-      reg_getDeformationFromDisplacement(deformationField);
+      NiftiImage deformationField;
+      reg_createDeformationField<float>(deformationField, averageImage);
       // Compute the transformation if required
       if(inputNRRName!=nullptr){
          nifti_image *current_transformation = reg_io_ReadImageFile(inputNRRName[i]);
@@ -465,7 +450,6 @@ int compute_average_image(nifti_image *averageImage,
                         nullptr,
                         interpolation_order,
                         std::numeric_limits<float>::quiet_NaN());
-      nifti_image_free(deformationField);
       nifti_image_free(current_input_image);
       // Add the image to the average
       remove_nan_and_add(averageImage, warpedImage, definedValue);

reg-apps/reg_resample.cpp

@@ -277,7 +277,7 @@ int main(int argc, char **argv)
    NR_VERBOSE_APP("Floating image name: " << floatingImage->fname);
    NR_VERBOSE_APP("\t" << floatingImage->nx << "x" << floatingImage->ny << "x" << floatingImage->nz << " voxels, " << floatingImage->nt << " volumes");
    NR_VERBOSE_APP("\t" << floatingImage->dx << "x" << floatingImage->dy << "x" << floatingImage->dz << " mm");
-   NR_VERBOSE_APP("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n");
+   NR_VERBOSE_APP("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *");
 
    /* *********************** */
    /* READ THE TRANSFORMATION */
@@ -313,9 +313,6 @@ int main(int argc, char **argv)
    // Create a deformation field
    nifti_image *deformationFieldImage = nifti_copy_nim_info(referenceImage);
    deformationFieldImage->dim[0]=deformationFieldImage->ndim=5;
-   deformationFieldImage->dim[1]=deformationFieldImage->nx=referenceImage->nx;
-   deformationFieldImage->dim[2]=deformationFieldImage->ny=referenceImage->ny;
-   deformationFieldImage->dim[3]=deformationFieldImage->nz=referenceImage->nz;
    deformationFieldImage->dim[4]=deformationFieldImage->nt=1;
    deformationFieldImage->pixdim[4]=deformationFieldImage->dt=1.0;
    deformationFieldImage->dim[5]=deformationFieldImage->nu=referenceImage->nz>1?3:2;
@@ -336,10 +333,14 @@ int main(int argc, char **argv)
    }
    deformationFieldImage->data = calloc(deformationFieldImage->nvox, deformationFieldImage->nbyper);
 
-   // Initialise the deformation field with an identity transformation
+   // Initialise as a displacement field with an identity transformation
+   deformationFieldImage->intent_code = NIFTI_INTENT_VECTOR;
+   memset(deformationFieldImage->intent_name, 0, 16);
+   strcpy(deformationFieldImage->intent_name, "NREG_TRANS");
+   deformationFieldImage->intent_p1 = DISP_FIELD;
    reg_tools_multiplyValueToImage(deformationFieldImage,deformationFieldImage,0.f);
+   // Convert it then to an deformation field with identity
    reg_getDeformationFromDisplacement(deformationFieldImage);
-   deformationFieldImage->intent_p1=DEF_FIELD;
 
    // Compute the transformation to apply
    if(inputTransformationImage!=nullptr)
@@ -348,40 +349,51 @@ int main(int argc, char **argv)
       {
       case LIN_SPLINE_GRID:
       case CUB_SPLINE_GRID:
-         reg_spline_getDeformationField(inputTransformationImage,
-                                        deformationFieldImage,
-                                        nullptr,
-                                        false,
-                                        true);
-         break;
+          NR_VERBOSE_APP("Input transformation is a cubic spline grid");
+          reg_spline_getDeformationField(inputTransformationImage,
+              deformationFieldImage,
+              nullptr, // no mask
+              true, // composition is used,
+              true); // b-spline are used
+          NR_VERBOSE_APP("Input transformation is converted to a deformation field");
+          break;
       case DISP_VEL_FIELD:
-         reg_getDeformationFromDisplacement(inputTransformationImage);
+          NR_VERBOSE_APP("Input transformation is a displacement velocity field");
+          reg_getDeformationFromDisplacement(inputTransformationImage);
+          NR_VERBOSE_APP("Input transformation is converted to a deformation velocity field");
       case DEF_VEL_FIELD:
          {
-            nifti_image *tempFlowField = nifti_dup(*deformationFieldImage);
-            reg_defField_compose(inputTransformationImage,
-                                 tempFlowField,
-                                 nullptr);
-            tempFlowField->intent_p1=inputTransformationImage->intent_p1;
-            tempFlowField->intent_p2=inputTransformationImage->intent_p2;
-            reg_defField_getDeformationFieldFromFlowField(tempFlowField,
-                                                          deformationFieldImage,
-                                                          false);
-            nifti_image_free(tempFlowField);
-         }
-         break;
+          NR_VERBOSE_APP("Input transformation is a deformation velocity field");
+          nifti_image *tempFlowField = nifti_dup(*deformationFieldImage);
+          reg_defField_compose(inputTransformationImage,
+                               tempFlowField,
+                               nullptr);
+          tempFlowField->intent_p1=inputTransformationImage->intent_p1;
+          tempFlowField->intent_p2=inputTransformationImage->intent_p2;
+          reg_defField_getDeformationFieldFromFlowField(tempFlowField,
+                                                        deformationFieldImage,
+                                                        false);
+          nifti_image_free(tempFlowField);
+          NR_VERBOSE_APP("Input transformation is converted to a deformation field");
+          }
+          break;
       case SPLINE_VEL_GRID:
-         reg_spline_getDefFieldFromVelocityGrid(inputTransformationImage,
+          NR_VERBOSE_APP("Input transformation is a spine velocity grid");
+          reg_spline_getDefFieldFromVelocityGrid(inputTransformationImage,
                                                 deformationFieldImage,
                                                 false);
-         break;
+          NR_VERBOSE_APP("Input transformation is converted to a deformation field");
+          break;
       case DISP_FIELD:
-         reg_getDeformationFromDisplacement(inputTransformationImage);
+          NR_VERBOSE_APP("Input transformation is a displacement field");
+          reg_getDeformationFromDisplacement(inputTransformationImage);
+          NR_VERBOSE_APP("Input transformation is converted to a deformation field");
       default:
-         reg_defField_compose(inputTransformationImage,
-                              deformationFieldImage,
-                              nullptr);
-         break;
+          NR_VERBOSE_APP("Input transformation is a deformation field");
+          reg_defField_compose(inputTransformationImage,
+                               deformationFieldImage,
+                               nullptr);
+          break;
       }
       nifti_image_free(inputTransformationImage);
       inputTransformationImage=nullptr;

reg-lib/Content.cpp

@@ -69,10 +69,14 @@ void Content::AllocateDeformationField(size_t bytes) {
     deformationField->intent_code = NIFTI_INTENT_VECTOR;
     memset(deformationField->intent_name, 0, sizeof(deformationField->intent_name));
     strcpy(deformationField->intent_name, "NREG_TRANS");
-    deformationField->intent_p1 = DEF_FIELD;
+    // First create a displacement field filled with 0 to obtain an identity disp
+    deformationField->intent_p1 = DISP_FIELD;
     deformationField->scl_slope = 1;
     deformationField->scl_inter = 0;
     deformationField->data = calloc(deformationField->nvox, deformationField->nbyper);
+    reg_tools_multiplyValueToImage(deformationField, deformationField, 0.f);
+    // Convert to an identity deformation field
+    reg_getDeformationFromDisplacement(deformationField);
 }
 /* *************************************************************** */
 void Content::DeallocateDeformationField() {

reg-lib/Debug.hpp

@@ -58,8 +58,8 @@ inline std::string StripFunctionName(const std::string& funcName) {
 #else
 #define NR_FUNC_CALLED()
 #define NR_DEBUG(msg)
-#define NR_VERBOSE(msg)     if (this->verbose) NR_COUT << "[NiftyReg DEBUG] " << msg << std::endl
-#define NR_VERBOSE_APP(msg) if (verbose) NR_COUT << "[NiftyReg DEBUG] " << msg << std::endl
+#define NR_VERBOSE(msg)     if (this->verbose) NR_COUT << "[NiftyReg INFO] " << msg << std::endl
+#define NR_VERBOSE_APP(msg) if (verbose) NR_COUT << "[NiftyReg INFO] " << msg << std::endl
 #endif
 /* *************************************************************** */
 #define NR_WARN(msg)        NR_COUT << "[NiftyReg WARNING] " << msg << std::endl
@@ -68,14 +68,19 @@ inline std::string StripFunctionName(const std::string& funcName) {
 #define NR_INFO(msg)        NR_COUT << "[NiftyReg INFO] " << msg << std::endl
 /* *************************************************************** */
 #ifndef NDEBUG
-#define NR_MAT33(mat, title)            reg_mat33_disp(mat, "[NiftyReg DEBUG] "s + (title))
-#define NR_MAT33_VERBOSE(mat, title)    NR_MAT33(mat, title)
-#define NR_MAT44(mat, title)            reg_mat44_disp(mat, "[NiftyReg DEBUG] "s + (title))
-#define NR_MAT44_VERBOSE(mat, title)    NR_MAT44(mat, title)
+#define NR_MAT33(mat, title)          reg_mat33_disp(mat, "[NiftyReg DEBUG] "s + (title))
+#define NR_MAT44(mat, title)          reg_mat44_disp(mat, "[NiftyReg DEBUG] "s + (title))
+#define NR_MAT33_DEBUG(mat, title)    NR_MAT33(mat, title)
+#define NR_MAT44_DEBUG(mat, title)    NR_MAT44(mat, title)
+#define NR_MAT33_VERBOSE(mat, title)  NR_MAT33(mat, title)
+#define NR_MAT44_VERBOSE(mat, title)  NR_MAT44(mat, title)
 #else
-#define NR_MAT33(mat, title)
-#define NR_MAT33_VERBOSE(mat, title)    if (this->verbose) reg_mat33_disp(mat, "[NiftyReg DEBUG] "s + (title))
-#define NR_MAT44(mat, title)
-#define NR_MAT44_VERBOSE(mat, title)    if (this->verbose) reg_mat44_disp(mat, "[NiftyReg DEBUG] "s + (title))
+#define NR_MAT33(mat, title)          reg_mat33_disp(mat, title)
+#define NR_MAT44(mat, title)          reg_mat44_disp(mat, title)
+#define NR_MAT33_DEBUG(mat, title)    
+#define NR_MAT44_DEBUG(mat, title)    
+#define NR_MAT33_VERBOSE(mat, title)  if (this->verbose) NR_MAT33(mat, "[NiftyReg INFO] "s + (title))
+#define NR_MAT44_VERBOSE(mat, title)  if (this->verbose) NR_MAT44(mat, "[NiftyReg INFO] "s + (title))
+
 #endif
 /* *************************************************************** */

reg-lib/_reg_aladin_sym.cpp

@@ -135,8 +135,8 @@ void reg_aladin_sym<T>::UpdateTransformationMatrix(int type) {
     this->bBlockMatchingKernel->template castTo<BlockMatchingKernel>()->Calculate();
     this->bLtsKernel->template castTo<LtsKernel>()->Calculate(type);
 
-    NR_MAT44_VERBOSE(*this->affineTransformation, "The pre-updated forward transformation matrix");
-    NR_MAT44_VERBOSE(*this->affineTransformationBw, "The pre-updated backward transformation matrix");
+    NR_MAT44_DEBUG(*this->affineTransformation, "The pre-updated forward transformation matrix");
+    NR_MAT44_DEBUG(*this->affineTransformationBw, "The pre-updated backward transformation matrix");
 
     // Forward and backward matrix are inverted
     mat44 fInverted = nifti_mat44_inverse(*this->affineTransformation);
@@ -153,8 +153,8 @@ void reg_aladin_sym<T>::UpdateTransformationMatrix(int type) {
     this->affineTransformation->m[3][3] = 1.f;
     this->affineTransformationBw->m[3][3] = 1.f;
 
-    NR_MAT44_VERBOSE(*this->affineTransformation, "The updated forward transformation matrix");
-    NR_MAT44_VERBOSE(*this->affineTransformationBw, "The updated backward transformation matrix");
+    NR_MAT44_DEBUG(*this->affineTransformation, "The updated forward transformation matrix");
+    NR_MAT44_DEBUG(*this->affineTransformationBw, "The updated backward transformation matrix");
 }
 /* *************************************************************** */
 template <class T>

reg-lib/_reg_f3d.cpp

@@ -176,13 +176,14 @@ void reg_f3d<T>::Initialise() {
         if (this->referencePyramid[0]->nz > 1)
             gridSpacing[2] = spacingInMillimetre[2] * powf(2, this->levelNumber - 1);
 
-        // Create and allocate the control point image
+        // Create and allocate the control point image - by default the transformation is initialised
+        // to an identity transformation
         reg_createControlPointGrid<T>(controlPointGrid, this->referencePyramid[0], gridSpacing);
 
-        // The control point position image is initialised with the affine transformation
-        if (!this->affineTransformation) {
-            reg_getDeformationFromDisplacement(controlPointGrid);
-        } else reg_affine_getDeformationField(this->affineTransformation.get(), controlPointGrid);
+        // The control point grid is updated with an identity transformation
+        if (this->affineTransformation) {
+            reg_affine_getDeformationField(this->affineTransformation.get(), controlPointGrid);
+        }
     } else {
         // The control point grid image is initialised with the provided grid
         controlPointGrid = inputControlPointGrid;
@@ -419,15 +420,16 @@ void reg_f3d<T>::DisplayCurrentLevelParameters(int currentLevel) {
     NR_VERBOSE("\t* image dimension: " << controlPointGrid->nx << " x " << controlPointGrid->ny << " x " << controlPointGrid->nz);
     NR_VERBOSE("\t* image spacing: " << controlPointGrid->dx << " x " << controlPointGrid->dy << " x " << controlPointGrid->dz << " mm");
 
+    // Input matrices are only printed out in debug
     if (reference->sform_code > 0)
-        NR_MAT44_VERBOSE(reference->sto_xyz, "Reference sform");
-    else NR_MAT44_VERBOSE(reference->qto_xyz, "Reference qform");
+        NR_MAT44_DEBUG(reference->sto_xyz, "Reference sform");
+    else NR_MAT44_DEBUG(reference->qto_xyz, "Reference qform");
     if (floating->sform_code > 0)
-        NR_MAT44_VERBOSE(floating->sto_xyz, "Floating sform");
-    else NR_MAT44_VERBOSE(floating->qto_xyz, "Floating qform");
+        NR_MAT44_DEBUG(floating->sto_xyz, "Floating sform");
+    else NR_MAT44_DEBUG(floating->qto_xyz, "Floating qform");
     if (controlPointGrid->sform_code > 0)
-        NR_MAT44_VERBOSE(controlPointGrid->sto_xyz, "CPP sform");
-    else NR_MAT44_VERBOSE(controlPointGrid->qto_xyz, "CPP qform");
+        NR_MAT44_DEBUG(controlPointGrid->sto_xyz, "CPP sform");
+    else NR_MAT44_DEBUG(controlPointGrid->qto_xyz, "CPP qform");
 
     NR_FUNC_CALLED();
 }

reg-lib/_reg_f3d2.cpp

@@ -422,8 +422,8 @@ void reg_f3d2<T>::DisplayCurrentLevelParameters(int currentLevel) {
     NR_VERBOSE("\t* image spacing: " << controlPointGridBw->dx << " x " << controlPointGridBw->dy << " x " << controlPointGridBw->dz << " mm");
 
     if (controlPointGridBw->sform_code > 0)
-        NR_MAT44_VERBOSE(controlPointGridBw->sto_xyz, "Backward CPP sform");
-    else NR_MAT44_VERBOSE(controlPointGridBw->qto_xyz, "Backward CPP qform");
+        NR_MAT44_DEBUG(controlPointGridBw->sto_xyz, "Backward CPP sform");
+    else NR_MAT44_DEBUG(controlPointGridBw->qto_xyz, "Backward CPP qform");
 
     NR_FUNC_CALLED();
 }

reg-lib/cpu/_reg_globalTrans.cpp

@@ -38,8 +38,6 @@ void reg_affine_deformationField2D(mat44 *affineTransformation,
       transformationMatrix = *affineTransformation;
    else transformationMatrix = reg_mat44_mul(affineTransformation, referenceMatrix);
 
-   NR_MAT44(transformationMatrix, "Global affine transformation");
-
    double voxel[3]={0,0,0}, position[3]={0,0,0};
    int x=0, y=0;
    size_t index=0;
@@ -99,8 +97,6 @@ void reg_affine_deformationField3D(mat44 *affineTransformation,
       transformationMatrix = *affineTransformation;
    else transformationMatrix = reg_mat44_mul(affineTransformation, referenceMatrix);
 
-   NR_MAT44(transformationMatrix, "Global affine transformation");
-
    double voxel[3]={0,0,0}, position[3]={0,0,0};
    int x=0, y=0, z=0;
    size_t index=0;