Fix a bug causing wrong calculation of the affine transformation matrix

niftyreg_build_version.txt

@@ -1 +1 @@
-265
+266

reg-apps/reg_aladin.cpp

@@ -163,9 +163,9 @@ int main(int argc, char **argv) {
     float floatingSigma = 0;
     float referenceSigma = 0;
 
-    float referenceLowerThr = std::numeric_limits<PrecisionType>::min();
+    float referenceLowerThr = std::numeric_limits<PrecisionType>::lowest();
     float referenceUpperThr = std::numeric_limits<PrecisionType>::max();
-    float floatingLowerThr = std::numeric_limits<PrecisionType>::min();
+    float floatingLowerThr = std::numeric_limits<PrecisionType>::lowest();
     float floatingUpperThr = std::numeric_limits<PrecisionType>::max();
     float paddingValue = std::numeric_limits<PrecisionType>::quiet_NaN();
 

reg-lib/_reg_aladin.cpp

@@ -32,11 +32,11 @@ reg_aladin<T>::reg_aladin() {
     this->floatingSigma = 0;
     this->referenceSigma = 0;
 
+    this->referenceLowerThreshold = std::numeric_limits<T>::lowest();
     this->referenceUpperThreshold = std::numeric_limits<T>::max();
-    this->referenceLowerThreshold = std::numeric_limits<T>::min();
 
+    this->floatingLowerThreshold = std::numeric_limits<T>::lowest();
     this->floatingUpperThreshold = std::numeric_limits<T>::max();
-    this->floatingLowerThreshold = std::numeric_limits<T>::min();
 
     this->warpedPaddingValue = std::numeric_limits<T>::quiet_NaN();
 

reg-lib/_reg_aladin_sym.cpp

@@ -11,9 +11,6 @@ reg_aladin_sym<T>::reg_aladin_sym()
 
     this->backwardBlockMatchingParams = nullptr;
 
-    this->floatingUpperThreshold = std::numeric_limits<T>::max();
-    this->floatingLowerThreshold = std::numeric_limits<T>::min();
-
 #ifndef NDEBUG
     reg_print_msg_debug("reg_aladin_sym constructor called");
 #endif
@@ -63,7 +60,7 @@ void reg_aladin_sym<T>::InitialiseRegistration() {
             }
         }
     }
-    if (this->floatingLowerThreshold != std::numeric_limits<T>::min()) {
+    if (this->floatingLowerThreshold != std::numeric_limits<T>::lowest()) {
         for (unsigned l = 0; l < this->levelsToPerform; ++l) {
             T *refPtr = static_cast<T *>(this->floatingPyramid[l]->data);
             int *mskPtr = this->floatingMaskPyramid[l].get();

reg-lib/_reg_base.cpp

@@ -34,14 +34,14 @@ reg_base<T>::reg_base(int refTimePoint, int floTimePoint) {
     referenceSmoothingSigma = 0;
     floatingSmoothingSigma = 0;
 
+    referenceThresholdLow.reset(new T[referenceTimePoint]);
+    std::fill(referenceThresholdLow.get(), referenceThresholdLow.get() + referenceTimePoint, std::numeric_limits<T>::lowest());
     referenceThresholdUp.reset(new T[referenceTimePoint]);
     std::fill(referenceThresholdUp.get(), referenceThresholdUp.get() + referenceTimePoint, std::numeric_limits<T>::max());
-    referenceThresholdLow.reset(new T[referenceTimePoint]);
-    std::fill(referenceThresholdLow.get(), referenceThresholdLow.get() + referenceTimePoint, std::numeric_limits<T>::min());
+    floatingThresholdLow.reset(new T[floatingTimePoint]);
+    std::fill(floatingThresholdLow.get(), floatingThresholdLow.get() + floatingTimePoint, std::numeric_limits<T>::lowest());
     floatingThresholdUp.reset(new T[floatingTimePoint]);
     std::fill(floatingThresholdUp.get(), floatingThresholdUp.get() + floatingTimePoint, std::numeric_limits<T>::max());
-    floatingThresholdLow.reset(new T[floatingTimePoint]);
-    std::fill(floatingThresholdLow.get(), floatingThresholdLow.get() + floatingTimePoint, std::numeric_limits<T>::min());
 
     robustRange = false;
     warpedPaddingValue = std::numeric_limits<T>::quiet_NaN();
@@ -504,7 +504,7 @@ void reg_base<T>::Initialise() {
         T *refDataPtr = static_cast<T *>(tmpReference->data);
         reg_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>::min())
+        if (referenceThresholdLow[0] == std::numeric_limits<T>::lowest())
             referenceThresholdLow[0] = refDataPtr[(int)reg_round((float)tmpReference->nvox * 0.02f)];
         if (referenceThresholdUp[0] == std::numeric_limits<T>::max())
             referenceThresholdUp[0] = refDataPtr[(int)reg_round((float)tmpReference->nvox * 0.98f)];
@@ -516,7 +516,7 @@ void reg_base<T>::Initialise() {
         T *floDataPtr = static_cast<T *>(tmpFloating->data);
         reg_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>::min())
+        if (floatingThresholdLow[0] == std::numeric_limits<T>::lowest())
             floatingThresholdLow[0] = floDataPtr[(int)reg_round((float)tmpFloating->nvox * 0.02f)];
         if (floatingThresholdUp[0] == std::numeric_limits<T>::max())
             floatingThresholdUp[0] = floDataPtr[(int)reg_round((float)tmpFloating->nvox * 0.98f)];

reg-lib/cpu/_reg_tools.cpp

@@ -104,35 +104,35 @@ void reg_intensityRescale_core(nifti_image *image,
     switch (image->datatype) {
     case NIFTI_TYPE_UINT8:
         currentMin = (DataType)std::numeric_limits<unsigned char>::max();
-        currentMax = 0;
+        currentMax = (DataType)std::numeric_limits<unsigned char>::lowest();
         break;
     case NIFTI_TYPE_INT8:
         currentMin = (DataType)std::numeric_limits<char>::max();
-        currentMax = (DataType)std::numeric_limits<char>::min();
+        currentMax = (DataType)std::numeric_limits<char>::lowest();
         break;
     case NIFTI_TYPE_UINT16:
         currentMin = (DataType)std::numeric_limits<unsigned short>::max();
-        currentMax = (DataType)std::numeric_limits<unsigned short>::min();
+        currentMax = (DataType)std::numeric_limits<unsigned short>::lowest();
         break;
     case NIFTI_TYPE_INT16:
         currentMin = (DataType)std::numeric_limits<short>::max();
-        currentMax = (DataType)std::numeric_limits<short>::min();
+        currentMax = (DataType)std::numeric_limits<short>::lowest();
         break;
     case NIFTI_TYPE_UINT32:
         currentMin = (DataType)std::numeric_limits<unsigned>::max();
-        currentMax = (DataType)std::numeric_limits<unsigned>::min();
+        currentMax = (DataType)std::numeric_limits<unsigned>::lowest();
         break;
     case NIFTI_TYPE_INT32:
         currentMin = (DataType)std::numeric_limits<int>::max();
-        currentMax = (DataType)std::numeric_limits<int>::min();
+        currentMax = (DataType)std::numeric_limits<int>::lowest();
         break;
     case NIFTI_TYPE_FLOAT32:
         currentMin = (DataType)std::numeric_limits<float>::max();
-        currentMax = (DataType)std::numeric_limits<float>::min();
+        currentMax = (DataType)std::numeric_limits<float>::lowest();
         break;
     case NIFTI_TYPE_FLOAT64:
         currentMin = (DataType)std::numeric_limits<double>::max();
-        currentMax = (DataType)std::numeric_limits<double>::min();
+        currentMax = (DataType)std::numeric_limits<double>::lowest();
         break;
     }
 
@@ -284,7 +284,7 @@ template<class T, class DataType>
 void reg_thresholdImage(nifti_image *image, T lowThr, T upThr) {
     DataType *imagePtr = static_cast<DataType*>(image->data);
     T currentMin = std::numeric_limits<T>::max();
-    T currentMax = std::numeric_limits<T>::min();
+    T currentMax = std::numeric_limits<T>::lowest();
 
     if (image->scl_slope == 0)image->scl_slope = 1.0;
 
@@ -1338,7 +1338,7 @@ void reg_tools_labelKernelConvolution_core(nifti_image *image,
                             }
                             currIterator = tmp_lab.begin();
                             maxindex = 0;
-                            maxval = std::numeric_limits<float>::min();
+                            maxval = std::numeric_limits<float>::lowest();
                             while (currIterator != tmp_lab.end()) {
                                 if (currIterator->second > maxval) {
                                     maxindex = currIterator->first;
@@ -2008,7 +2008,7 @@ DataType reg_tools_getMinMaxValue(const nifti_image *image, int timepoint, bool
         reg_print_msg_error("reg_tools_getMinMaxValue. The required time point does not exists");
 
     const DataType *imgPtr = static_cast<DataType*>(image->data);
-    DataType retValue = calcMin ? std::numeric_limits<DataType>::max() : std::numeric_limits<DataType>::min();
+    DataType retValue = calcMin ? std::numeric_limits<DataType>::max() : std::numeric_limits<DataType>::lowest();
     const size_t voxelNumber = CalcVoxelNumber(*image);
     const float sclSlope = image->scl_slope == 0 ? 1 : image->scl_slope;