From a447e474dfda8c3428deebe8e0c4a2fbe63785a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Onur=20=C3=9Clgen?= Date: Mon, 31 Jul 2023 14:06:30 +0100 Subject: [PATCH] Rearchitect reg_measure to handle forward and backward similarity measure values #92 --- niftyreg_build_version.txt | 2 +- reg-apps/reg_tools.cpp | 4 +- reg-lib/cpu/_reg_dti.cpp | 582 ++++++++++++++------------------ reg-lib/cpu/_reg_dti.h | 14 +- reg-lib/cpu/_reg_kld.cpp | 154 +++------ reg-lib/cpu/_reg_kld.h | 18 +- reg-lib/cpu/_reg_lncc.cpp | 391 ++++++++++----------- reg-lib/cpu/_reg_lncc.h | 70 +--- reg-lib/cpu/_reg_measure.h | 41 ++- reg-lib/cpu/_reg_mind.cpp | 487 ++++++++++++-------------- reg-lib/cpu/_reg_mind.h | 50 +-- reg-lib/cpu/_reg_nmi.cpp | 120 ++++--- reg-lib/cpu/_reg_nmi.h | 14 +- reg-lib/cpu/_reg_ssd.cpp | 226 ++++++------- reg-lib/cpu/_reg_ssd.h | 49 +-- reg-lib/cuda/_reg_measure_gpu.h | 21 +- reg-lib/cuda/_reg_nmi_gpu.cu | 101 +++--- reg-lib/cuda/_reg_nmi_gpu.h | 12 +- reg-lib/cuda/_reg_ssd_gpu.cu | 6 +- reg-lib/cuda/_reg_ssd_gpu.h | 6 +- 20 files changed, 1067 insertions(+), 1301 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 594cd09d..9530e048 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -295 +296 diff --git a/reg-apps/reg_tools.cpp b/reg-apps/reg_tools.cpp index 4f2ea7b8..5c1d5eeb 100755 --- a/reg-apps/reg_tools.cpp +++ b/reg-apps/reg_tools.cpp @@ -1037,7 +1037,7 @@ int main(int argc, char **argv) outputImage->data = malloc(outputImage->nvox * outputImage->nbyper); // Compute the MIND descriptor int *mask = (int *)calloc(image->nvox, sizeof(int)); - GetMINDImageDescriptor(image, outputImage, mask, 1, 0); + GetMindImageDescriptor(image, outputImage, mask, 1, 0); free(mask); // Save the MIND descriptor image if(flag->outputImageFlag) @@ -1064,7 +1064,7 @@ int main(int argc, char **argv) outputImage->data = malloc(outputImage->nvox * outputImage->nbyper); // Compute the MIND-SSC descriptor int *mask = (int *)calloc(image->nvox, sizeof(int)); - GetMINDSSCImageDescriptor(image, outputImage, mask, 1, 0); + GetMindSscImageDescriptor(image, outputImage, mask, 1, 0); free(mask); // Save the MIND descriptor image if(flag->outputImageFlag) diff --git a/reg-lib/cpu/_reg_dti.cpp b/reg-lib/cpu/_reg_dti.cpp index e9c99a2f..d4fa63be 100755 --- a/reg-lib/cpu/_reg_dti.cpp +++ b/reg-lib/cpu/_reg_dti.cpp @@ -13,11 +13,9 @@ #include "_reg_dti.h" /* *************************************************************** */ -reg_dti::reg_dti() - : reg_measure() -{ +reg_dti::reg_dti(): reg_measure() { #ifndef NDEBUG - reg_print_msg_debug("reg_dti constructor called"); + reg_print_msg_debug("reg_dti constructor called"); #endif } /* *************************************************************** */ @@ -32,89 +30,82 @@ void reg_dti::InitialiseMeasure(nifti_image *refImg, int *floMask, nifti_image *warpedImgBw, nifti_image *warpedGradBw, - nifti_image *voxelBasedGradBw) -{ - // Set the pointers using the parent class function - reg_measure::InitialiseMeasure(refImg, - floImg, - refMask, - warpedImg, - warpedGrad, - voxelBasedGrad, - localWeightSim, - floMask, - warpedImgBw, - warpedGradBw, - voxelBasedGradBw); + nifti_image *voxelBasedGradBw) { + // Set the pointers using the parent class function + reg_measure::InitialiseMeasure(refImg, + floImg, + refMask, + warpedImg, + warpedGrad, + voxelBasedGrad, + localWeightSim, + floMask, + warpedImgBw, + warpedGradBw, + voxelBasedGradBw); - // Check that the input images have the same number of time point - if(this->referenceImage->nt != this->floatingImage->nt) - { - reg_print_fct_error("reg_dti::InitialiseMeasure"); - reg_print_msg_error("This number of time point should be the same for both input images"); - reg_exit(); - } + // Check that the input images have the same number of time point + if (this->referenceImage->nt != this->floatingImage->nt) { + reg_print_fct_error("reg_dti::InitialiseMeasure"); + reg_print_msg_error("This number of time point should be the same for both input images"); + reg_exit(); + } - int j=0; - for(int i=0; int; ++i) - { - //JM - note, the specific value of timePointWeight is not used for DTI images - //any value > 0 indicates the 'time point' is active - if(this->timePointWeight[i]>0) - { - this->dtIndicies[j++]=i; + int j = 0; + for (int i = 0; i < refImg->nt; ++i) { + //JM - note, the specific value of timePointWeight is not used for DTI images + //any value > 0 indicates the 'time point' is active + if (this->timePointWeight[i] > 0) { + this->dtIndicies[j++] = i; #ifndef NDEBUG - reg_print_msg_debug("reg_dti::InitialiseMeasure()."); - char text[255]; - sprintf(text, "Active time point: %i", i); - reg_print_msg_debug(text); + reg_print_msg_debug("reg_dti::InitialiseMeasure()"); + char text[255]; + sprintf(text, "Active time point: %i", i); + reg_print_msg_debug(text); #endif - } - } - if((refImg->nz>1 && j!=6) && (refImg->nz==1 && j!=3)) - { - reg_print_fct_error("reg_dti::InitialiseMeasure"); - reg_print_msg_error("Unexpected number of DTI components"); - reg_exit(); - } + } + } + if ((refImg->nz > 1 && j != 6) && (refImg->nz == 1 && j != 3)) { + reg_print_fct_error("reg_dti::InitialiseMeasure"); + reg_print_msg_error("Unexpected number of DTI components"); + reg_exit(); + } } /* *************************************************************** */ template -double reg_getDTIMeasureValue(nifti_image *referenceImage, - nifti_image *warpedImage, - int *mask, - unsigned *dtIndicies - ) -{ +double reg_getDTIMeasureValue(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const int *mask, + const unsigned *dtIndicies) { #ifdef _WIN32 - long voxel; - const long voxelNumber = (long)NiftiImage::calcVoxelNumber(referenceImage, 3); + long voxel; + const long voxelNumber = (long)NiftiImage::calcVoxelNumber(referenceImage, 3); #else - size_t voxel; - const size_t voxelNumber = NiftiImage::calcVoxelNumber(referenceImage, 3); + size_t voxel; + const size_t voxelNumber = NiftiImage::calcVoxelNumber(referenceImage, 3); #endif - /* As the tensor has 6 unique components that we need to worry about, read them out - for the floating and reference images. */ - DataType *firstWarpedVox = static_cast(warpedImage->data); - DataType *warpedIntensityXX = &firstWarpedVox[voxelNumber*dtIndicies[0]]; - DataType *warpedIntensityXY = &firstWarpedVox[voxelNumber*dtIndicies[1]]; - DataType *warpedIntensityYY = &firstWarpedVox[voxelNumber*dtIndicies[2]]; - DataType *warpedIntensityXZ = &firstWarpedVox[voxelNumber*dtIndicies[3]]; - DataType *warpedIntensityYZ = &firstWarpedVox[voxelNumber*dtIndicies[4]]; - DataType *warpedIntensityZZ = &firstWarpedVox[voxelNumber*dtIndicies[5]]; + /* As the tensor has 6 unique components that we need to worry about, read them out + for the floating and reference images. */ + const DataType *firstWarpedVox = static_cast(warpedImage->data); + const DataType *warpedIntensityXX = &firstWarpedVox[voxelNumber * dtIndicies[0]]; + const DataType *warpedIntensityXY = &firstWarpedVox[voxelNumber * dtIndicies[1]]; + const DataType *warpedIntensityYY = &firstWarpedVox[voxelNumber * dtIndicies[2]]; + const DataType *warpedIntensityXZ = &firstWarpedVox[voxelNumber * dtIndicies[3]]; + const DataType *warpedIntensityYZ = &firstWarpedVox[voxelNumber * dtIndicies[4]]; + const DataType *warpedIntensityZZ = &firstWarpedVox[voxelNumber * dtIndicies[5]]; - DataType *firstRefVox = static_cast(referenceImage->data); - DataType *referenceIntensityXX = &firstRefVox[voxelNumber*dtIndicies[0]]; - DataType *referenceIntensityXY = &firstRefVox[voxelNumber*dtIndicies[1]]; - DataType *referenceIntensityYY = &firstRefVox[voxelNumber*dtIndicies[2]]; - DataType *referenceIntensityXZ = &firstRefVox[voxelNumber*dtIndicies[3]]; - DataType *referenceIntensityYZ = &firstRefVox[voxelNumber*dtIndicies[4]]; - DataType *referenceIntensityZZ = &firstRefVox[voxelNumber*dtIndicies[5]]; + const DataType *firstRefVox = static_cast(referenceImage->data); + const DataType *referenceIntensityXX = &firstRefVox[voxelNumber * dtIndicies[0]]; + const DataType *referenceIntensityXY = &firstRefVox[voxelNumber * dtIndicies[1]]; + const DataType *referenceIntensityYY = &firstRefVox[voxelNumber * dtIndicies[2]]; + const DataType *referenceIntensityXZ = &firstRefVox[voxelNumber * dtIndicies[3]]; + const DataType *referenceIntensityYZ = &firstRefVox[voxelNumber * dtIndicies[4]]; + const DataType *referenceIntensityZZ = &firstRefVox[voxelNumber * dtIndicies[5]]; - double DTI_cost=0, n=0; - const double twoThirds = (2.0/3.0); - DataType rXX, rXY, rYY, rXZ, rYZ, rZZ; + double dtiCost = 0, n = 0; + constexpr double twoThirds = 2.0 / 3.0; + DataType rXX, rXY, rYY, rXZ, rYZ, rZZ; #ifdef _OPENMP #pragma omp parallel for default(none) \ shared(referenceImage, referenceIntensityXX, referenceIntensityXY, referenceIntensityXZ, \ @@ -122,115 +113,65 @@ double reg_getDTIMeasureValue(nifti_image *referenceImage, warpedIntensityXX,warpedIntensityXY,warpedIntensityXZ, \ warpedIntensityYY,warpedIntensityYZ, warpedIntensityZZ, mask,voxelNumber) \ private(rXX, rXY, rYY, rXZ, rYZ, rZZ) \ - reduction(+:DTI_cost, n) + reduction(+:dtiCost, n) #endif - for(voxel=0; voxel-1 ) - { - if(referenceIntensityXX[voxel]==referenceIntensityXX[voxel] && - warpedIntensityXX[voxel]==warpedIntensityXX[voxel]) - { - // Calculate the elementwise residual of the diffusion tensor components - rXX = referenceIntensityXX[voxel] - warpedIntensityXX[voxel]; - rXY = referenceIntensityXY[voxel] - warpedIntensityXY[voxel]; - rYY = referenceIntensityYY[voxel] - warpedIntensityYY[voxel]; - rXZ = referenceIntensityXZ[voxel] - warpedIntensityXZ[voxel]; - rYZ = referenceIntensityYZ[voxel] - warpedIntensityYZ[voxel]; - rZZ = referenceIntensityZZ[voxel] - warpedIntensityZZ[voxel]; - DTI_cost -= twoThirds * (reg_pow2(rXX) + reg_pow2(rYY) + reg_pow2(rZZ)) - + 2.0 * (reg_pow2(rXY) + reg_pow2(rXZ) + reg_pow2(rYZ)) - - twoThirds * (rXX*rYY+rXX*rZZ+rYY*rZZ); - n++; - } // check if values are defined - } // check if voxel belongs mask - } // loop over voxels - return DTI_cost/n; + for (voxel = 0; voxel < voxelNumber; ++voxel) { + // Check if the current voxel belongs to the mask and the intensities are not nans + if (mask[voxel] > -1) { + if (referenceIntensityXX[voxel] == referenceIntensityXX[voxel] && + warpedIntensityXX[voxel] == warpedIntensityXX[voxel]) { + // Calculate the elementwise residual of the diffusion tensor components + rXX = referenceIntensityXX[voxel] - warpedIntensityXX[voxel]; + rXY = referenceIntensityXY[voxel] - warpedIntensityXY[voxel]; + rYY = referenceIntensityYY[voxel] - warpedIntensityYY[voxel]; + rXZ = referenceIntensityXZ[voxel] - warpedIntensityXZ[voxel]; + rYZ = referenceIntensityYZ[voxel] - warpedIntensityYZ[voxel]; + rZZ = referenceIntensityZZ[voxel] - warpedIntensityZZ[voxel]; + dtiCost -= twoThirds * (reg_pow2(rXX) + reg_pow2(rYY) + reg_pow2(rZZ)) + + 2.0 * (reg_pow2(rXY) + reg_pow2(rXZ) + reg_pow2(rYZ)) + - twoThirds * (rXX * rYY + rXX * rZZ + rYY * rZZ); + n++; + } // check if values are defined + } // check if voxel belongs mask + } // loop over voxels + return dtiCost / n; } -template double reg_getDTIMeasureValue(nifti_image *,nifti_image *,int *, unsigned *); -template double reg_getDTIMeasureValue(nifti_image *,nifti_image *,int *, unsigned *); /* *************************************************************** */ -double reg_dti::GetSimilarityMeasureValue() -{ - // Check that all the specified image are of the same datatype - if(this->warpedImage->datatype != this->referenceImage->datatype) - { - reg_print_fct_error("reg_dti::GetSimilarityMeasureValue"); - reg_print_msg_error("Both input images are expected to have the same type"); - reg_exit(); - } - double DTIMeasureValue; - switch(this->referenceImage->datatype) - { - case NIFTI_TYPE_FLOAT32: - DTIMeasureValue = reg_getDTIMeasureValue - (this->referenceImage, - this->warpedImage, - this->referenceMask, - this->dtIndicies - ); - break; - case NIFTI_TYPE_FLOAT64: - DTIMeasureValue = reg_getDTIMeasureValue - (this->referenceImage, - this->warpedImage, - this->referenceMask, - this->dtIndicies - ); - break; - default: - reg_print_fct_error("reg_dti::GetSimilarityMeasureValue"); - reg_print_msg_error("Result pixel type unsupported in the DTI computation function"); - reg_exit(); - } - - // Backward computation - if(this->isSymmetric) - { - // Check that all the specified image are of the same datatype - if(this->warpedImageBw->datatype != this->floatingImage->datatype) - { - reg_print_fct_error("reg_dti::GetSimilarityMeasureValue"); - reg_print_msg_error("Both input images are expected to have the same type"); - reg_exit(); - } - switch(this->floatingImage->datatype) - { - case NIFTI_TYPE_FLOAT32: - DTIMeasureValue += reg_getDTIMeasureValue - (this->floatingImage, - this->warpedImageBw, - this->floatingMask, - this->dtIndicies - ); - break; - case NIFTI_TYPE_FLOAT64: - DTIMeasureValue += reg_getDTIMeasureValue - (this->floatingImage, - this->warpedImageBw, - this->floatingMask, - this->dtIndicies - ); - break; - default: - reg_print_fct_error("reg_dti::GetSimilarityMeasureValue"); - reg_print_msg_error("Warped pixel type unsupported in the DTI computation function"); - reg_exit(); - } - } - return DTIMeasureValue; +double GetSimilarityMeasureValue(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const int *mask, + const unsigned *dtIndicies) { + return std::visit([&](auto&& refImgDataType) { + using RefImgDataType = std::decay_t; + return reg_getDTIMeasureValue(referenceImage, + warpedImage, + mask, + dtIndicies); + }, NiftiImage::getFloatingDataType(referenceImage)); +} +/* *************************************************************** */ +double reg_dti::GetSimilarityMeasureValueFw() { + return ::GetSimilarityMeasureValue(this->referenceImage, + this->warpedImage, + this->referenceMask, + this->dtIndicies); +} +/* *************************************************************** */ +double reg_dti::GetSimilarityMeasureValueBw() { + return ::GetSimilarityMeasureValue(this->floatingImage, + this->warpedImageBw, + this->floatingMask, + this->dtIndicies); } /* *************************************************************** */ template void reg_getVoxelBasedDTIMeasureGradient(nifti_image *referenceImage, - nifti_image *warpedImage, - nifti_image *warpedGradient, - nifti_image *dtiMeasureGradientImage, - int *mask, - unsigned *dtIndicies) -{ - // Create pointers to the reference and warped images + nifti_image *warpedImage, + nifti_image *warpedGradient, + nifti_image *dtiMeasureGradientImage, + int *mask, + unsigned *dtIndicies) { + // Create pointers to the reference and warped images #ifdef _WIN32 long voxel; const long voxelNumber = (long)NiftiImage::calcVoxelNumber(referenceImage, 3); @@ -239,45 +180,45 @@ void reg_getVoxelBasedDTIMeasureGradient(nifti_image *referenceImage, const size_t voxelNumber = NiftiImage::calcVoxelNumber(referenceImage, 3); #endif - /* As the tensor has 6 unique components that we need to worry about, read them out - for the floating and reference images. */ - DataType *firstWarpedVox = static_cast(warpedImage->data); - DataType *warpedIntensityXX = &firstWarpedVox[voxelNumber*dtIndicies[0]]; - DataType *warpedIntensityXY = &firstWarpedVox[voxelNumber*dtIndicies[1]]; - DataType *warpedIntensityYY = &firstWarpedVox[voxelNumber*dtIndicies[2]]; - DataType *warpedIntensityXZ = &firstWarpedVox[voxelNumber*dtIndicies[3]]; - DataType *warpedIntensityYZ = &firstWarpedVox[voxelNumber*dtIndicies[4]]; - DataType *warpedIntensityZZ = &firstWarpedVox[voxelNumber*dtIndicies[5]]; + /* As the tensor has 6 unique components that we need to worry about, read them out + for the floating and reference images. */ + DataType *firstWarpedVox = static_cast(warpedImage->data); + DataType *warpedIntensityXX = &firstWarpedVox[voxelNumber * dtIndicies[0]]; + DataType *warpedIntensityXY = &firstWarpedVox[voxelNumber * dtIndicies[1]]; + DataType *warpedIntensityYY = &firstWarpedVox[voxelNumber * dtIndicies[2]]; + DataType *warpedIntensityXZ = &firstWarpedVox[voxelNumber * dtIndicies[3]]; + DataType *warpedIntensityYZ = &firstWarpedVox[voxelNumber * dtIndicies[4]]; + DataType *warpedIntensityZZ = &firstWarpedVox[voxelNumber * dtIndicies[5]]; - DataType *firstRefVox = static_cast(referenceImage->data); - DataType *referenceIntensityXX = &firstRefVox[voxelNumber*dtIndicies[0]]; - DataType *referenceIntensityXY = &firstRefVox[voxelNumber*dtIndicies[1]]; - DataType *referenceIntensityYY = &firstRefVox[voxelNumber*dtIndicies[2]]; - DataType *referenceIntensityXZ = &firstRefVox[voxelNumber*dtIndicies[3]]; - DataType *referenceIntensityYZ = &firstRefVox[voxelNumber*dtIndicies[4]]; - DataType *referenceIntensityZZ = &firstRefVox[voxelNumber*dtIndicies[5]]; + DataType *firstRefVox = static_cast(referenceImage->data); + DataType *referenceIntensityXX = &firstRefVox[voxelNumber * dtIndicies[0]]; + DataType *referenceIntensityXY = &firstRefVox[voxelNumber * dtIndicies[1]]; + DataType *referenceIntensityYY = &firstRefVox[voxelNumber * dtIndicies[2]]; + DataType *referenceIntensityXZ = &firstRefVox[voxelNumber * dtIndicies[3]]; + DataType *referenceIntensityYZ = &firstRefVox[voxelNumber * dtIndicies[4]]; + DataType *referenceIntensityZZ = &firstRefVox[voxelNumber * dtIndicies[5]]; - // THE FOLLOWING IS WRONG - reg_print_msg_error("ERROR IN THE DTI GRADIENT COMPUTATION - TO FIX"); - reg_exit(); - unsigned gradientVoxels = warpedGradient->nu*voxelNumber; - DataType *firstGradVox = static_cast(warpedGradient->data); - DataType *spatialGradXX = &firstGradVox[gradientVoxels*dtIndicies[0]]; - DataType *spatialGradXY = &firstGradVox[gradientVoxels*dtIndicies[1]]; - DataType *spatialGradYY = &firstGradVox[gradientVoxels*dtIndicies[2]]; - DataType *spatialGradXZ = &firstGradVox[gradientVoxels*dtIndicies[3]]; - DataType *spatialGradYZ = &firstGradVox[gradientVoxels*dtIndicies[4]]; - DataType *spatialGradZZ = &firstGradVox[gradientVoxels*dtIndicies[5]]; + // THE FOLLOWING IS WRONG + reg_print_msg_error("ERROR IN THE DTI GRADIENT COMPUTATION - TO FIX"); + reg_exit(); + unsigned gradientVoxels = warpedGradient->nu * voxelNumber; + DataType *firstGradVox = static_cast(warpedGradient->data); + DataType *spatialGradXX = &firstGradVox[gradientVoxels * dtIndicies[0]]; + DataType *spatialGradXY = &firstGradVox[gradientVoxels * dtIndicies[1]]; + DataType *spatialGradYY = &firstGradVox[gradientVoxels * dtIndicies[2]]; + DataType *spatialGradXZ = &firstGradVox[gradientVoxels * dtIndicies[3]]; + DataType *spatialGradYZ = &firstGradVox[gradientVoxels * dtIndicies[4]]; + DataType *spatialGradZZ = &firstGradVox[gradientVoxels * dtIndicies[5]]; - // Create an array to store the computed gradient per time point - DataType *dtiMeasureGradPtrX=static_cast(dtiMeasureGradientImage->data); - DataType *dtiMeasureGradPtrY = &dtiMeasureGradPtrX[voxelNumber]; - DataType *dtiMeasureGradPtrZ = &dtiMeasureGradPtrY[voxelNumber]; + // Create an array to store the computed gradient per time point + DataType *dtiMeasureGradPtrX = static_cast(dtiMeasureGradientImage->data); + DataType *dtiMeasureGradPtrY = &dtiMeasureGradPtrX[voxelNumber]; + DataType *dtiMeasureGradPtrZ = &dtiMeasureGradPtrY[voxelNumber]; - const double twoThirds = 2.0/3.0; - const double fourThirds = 4.0/3.0; + const double twoThirds = 2.0 / 3.0; + const double fourThirds = 4.0 / 3.0; - DataType rXX, rXY, rYY, rXZ, rYZ, rZZ, xxGrad, yyGrad, zzGrad, xyGrad, xzGrad, yzGrad; + DataType rXX, rXY, rYY, rXZ, rYZ, rZZ, xxGrad, yyGrad, zzGrad, xyGrad, xzGrad, yzGrad; #ifdef _OPENMP #pragma omp parallel for default(none) \ shared(referenceIntensityXX, referenceIntensityXY, referenceIntensityXZ, \ @@ -287,133 +228,114 @@ void reg_getVoxelBasedDTIMeasureGradient(nifti_image *referenceImage, dtiMeasureGradPtrX, dtiMeasureGradPtrY, dtiMeasureGradPtrZ, voxelNumber) \ private(rXX, rXY, rYY, rXZ, rYZ, rZZ, xxGrad, yyGrad, zzGrad, xyGrad, xzGrad, yzGrad) #endif - for(voxel=0; voxel-1 ) - { - if(referenceIntensityXX[voxel]==referenceIntensityXX[voxel] && - warpedIntensityXX[voxel]==warpedIntensityXX[voxel]) - { - rXX = referenceIntensityXX[voxel] - warpedIntensityXX[voxel]; - rXY = referenceIntensityXY[voxel] - warpedIntensityXY[voxel]; - rYY = referenceIntensityYY[voxel] - warpedIntensityYY[voxel]; - rXZ = referenceIntensityXZ[voxel] - warpedIntensityXZ[voxel]; - rYZ = referenceIntensityYZ[voxel] - warpedIntensityYZ[voxel]; - rZZ = referenceIntensityZZ[voxel] - warpedIntensityZZ[voxel]; + for (voxel = 0; voxel < voxelNumber; voxel++) { + if (mask[voxel] > -1) { + if (referenceIntensityXX[voxel] == referenceIntensityXX[voxel] && + warpedIntensityXX[voxel] == warpedIntensityXX[voxel]) { + rXX = referenceIntensityXX[voxel] - warpedIntensityXX[voxel]; + rXY = referenceIntensityXY[voxel] - warpedIntensityXY[voxel]; + rYY = referenceIntensityYY[voxel] - warpedIntensityYY[voxel]; + rXZ = referenceIntensityXZ[voxel] - warpedIntensityXZ[voxel]; + rYZ = referenceIntensityYZ[voxel] - warpedIntensityYZ[voxel]; + rZZ = referenceIntensityZZ[voxel] - warpedIntensityZZ[voxel]; - xxGrad = fourThirds*rXX-twoThirds*(rYY+rZZ); - yyGrad = fourThirds*rYY-twoThirds*(rXX+rZZ); - zzGrad = fourThirds*rZZ-twoThirds*(rYY+rXX); - xyGrad = 4.0*rXY; - xzGrad = 4.0*rXZ; - yzGrad = 4.0*rYZ; + xxGrad = static_cast(fourThirds * rXX - twoThirds * (rYY + rZZ)); + yyGrad = static_cast(fourThirds * rYY - twoThirds * (rXX + rZZ)); + zzGrad = static_cast(fourThirds * rZZ - twoThirds * (rYY + rXX)); + xyGrad = 4.f * rXY; + xzGrad = 4.f * rXZ; + yzGrad = 4.f * rYZ; - dtiMeasureGradPtrX[voxel] -= (spatialGradXX[voxel]*xxGrad+spatialGradYY[voxel]*yyGrad+spatialGradZZ[voxel]*zzGrad \ - + spatialGradXY[voxel]*xyGrad + spatialGradXZ[voxel]*xzGrad + spatialGradYZ[voxel]*yzGrad); + dtiMeasureGradPtrX[voxel] -= (spatialGradXX[voxel] * xxGrad + spatialGradYY[voxel] * yyGrad + spatialGradZZ[voxel] * zzGrad + + spatialGradXY[voxel] * xyGrad + spatialGradXZ[voxel] * xzGrad + spatialGradYZ[voxel] * yzGrad); - dtiMeasureGradPtrY[voxel] -= (spatialGradXX[voxel+voxelNumber]*xxGrad+spatialGradYY[voxel+voxelNumber]*yyGrad+spatialGradZZ[voxel+voxelNumber]*zzGrad \ - + spatialGradXY[voxel+voxelNumber]*xyGrad + spatialGradXZ[voxel+voxelNumber]*xzGrad + spatialGradYZ[voxel+voxelNumber]*yzGrad); + dtiMeasureGradPtrY[voxel] -= (spatialGradXX[voxel + voxelNumber] * xxGrad + spatialGradYY[voxel + voxelNumber] * yyGrad + spatialGradZZ[voxel + voxelNumber] * zzGrad + + spatialGradXY[voxel + voxelNumber] * xyGrad + spatialGradXZ[voxel + voxelNumber] * xzGrad + spatialGradYZ[voxel + voxelNumber] * yzGrad); - dtiMeasureGradPtrZ[voxel] -= (spatialGradXX[voxel+2*voxelNumber]*xxGrad+spatialGradYY[voxel+2*voxelNumber]*yyGrad \ - + spatialGradZZ[voxel+2*voxelNumber]*zzGrad + spatialGradXY[voxel+2*voxelNumber]*xyGrad \ - + spatialGradXZ[voxel+2*voxelNumber]*xzGrad + spatialGradYZ[voxel+2*voxelNumber]*yzGrad); - } - } - } + dtiMeasureGradPtrZ[voxel] -= (spatialGradXX[voxel + 2 * voxelNumber] * xxGrad + spatialGradYY[voxel + 2 * voxelNumber] * yyGrad + + spatialGradZZ[voxel + 2 * voxelNumber] * zzGrad + spatialGradXY[voxel + 2 * voxelNumber] * xyGrad + + spatialGradXZ[voxel + 2 * voxelNumber] * xzGrad + spatialGradYZ[voxel + 2 * voxelNumber] * yzGrad); + } + } + } } /* *************************************************************** */ -template void reg_getVoxelBasedDTIMeasureGradient -(nifti_image *,nifti_image *,nifti_image *,nifti_image *, int *, unsigned *); -template void reg_getVoxelBasedDTIMeasureGradient -(nifti_image *,nifti_image *,nifti_image *,nifti_image *, int *, unsigned *); -/* *************************************************************** */ -void reg_dti::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) -{ - // Check if the specified time point exists and is active - reg_measure::GetVoxelBasedSimilarityMeasureGradient(currentTimepoint); - if(this->timePointWeight[currentTimepoint]==0) - return; +void reg_dti::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { + // Check if the specified time point exists and is active + reg_measure::GetVoxelBasedSimilarityMeasureGradient(currentTimepoint); + if (this->timePointWeight[currentTimepoint] == 0) + return; - // Check if all required input images are of the same data type - int dtype = this->referenceImage->datatype; - if(this->warpedImage->datatype != dtype || - this->warpedGradient->datatype != dtype || - this->voxelBasedGradient->datatype != dtype - ) - { - reg_print_fct_error("reg_dti::GetVoxelBasedSimilarityMeasureGradient"); - reg_print_msg_error("Input images are expected to be of the same type"); - reg_exit(); - } - // Compute the gradient of the ssd for the forward transformation - switch(dtype) - { - case NIFTI_TYPE_FLOAT32: - reg_getVoxelBasedDTIMeasureGradient - (this->referenceImage, - this->warpedImage, - this->warpedGradient, - this->voxelBasedGradient, - this->referenceMask, - this->dtIndicies - ); - break; - case NIFTI_TYPE_FLOAT64: - reg_getVoxelBasedDTIMeasureGradient - (this->referenceImage, - this->warpedImage, - this->warpedGradient, - this->voxelBasedGradient, - this->referenceMask, - this->dtIndicies - ); - break; - default: - reg_print_fct_error("reg_dti::GetVoxelBasedSimilarityMeasureGradient"); - reg_print_msg_error("The input image data type is not supported"); - reg_exit(); - } - // Compute the gradient of the ssd for the backward transformation - if(this->isSymmetric) - { - dtype = this->floatingImage->datatype; - if(this->warpedImageBw->datatype != dtype || + // Check if all required input images are of the same data type + int dtype = this->referenceImage->datatype; + if (this->warpedImage->datatype != dtype || + this->warpedGradient->datatype != dtype || + this->voxelBasedGradient->datatype != dtype + ) { + reg_print_fct_error("reg_dti::GetVoxelBasedSimilarityMeasureGradient"); + reg_print_msg_error("Input images are expected to be of the same type"); + reg_exit(); + } + // Compute the gradient of the ssd for the forward transformation + switch (dtype) { + case NIFTI_TYPE_FLOAT32: + reg_getVoxelBasedDTIMeasureGradient + (this->referenceImage, + this->warpedImage, + this->warpedGradient, + this->voxelBasedGradient, + this->referenceMask, + this->dtIndicies); + break; + case NIFTI_TYPE_FLOAT64: + reg_getVoxelBasedDTIMeasureGradient + (this->referenceImage, + this->warpedImage, + this->warpedGradient, + this->voxelBasedGradient, + this->referenceMask, + this->dtIndicies); + break; + default: + reg_print_fct_error("reg_dti::GetVoxelBasedSimilarityMeasureGradient"); + reg_print_msg_error("The input image data type is not supported"); + reg_exit(); + } + // Compute the gradient of the ssd for the backward transformation + if (this->isSymmetric) { + dtype = this->floatingImage->datatype; + if (this->warpedImageBw->datatype != dtype || this->warpedGradientBw->datatype != dtype || - this->voxelBasedGradientBw->datatype != dtype - ) - { - reg_print_fct_error("reg_dti::GetVoxelBasedSimilarityMeasureGradient"); - reg_print_msg_error("Input images are expected to be of the same type"); - reg_exit(); - } - // Compute the gradient of the nmi for the backward transformation - switch(dtype) - { - case NIFTI_TYPE_FLOAT32: - reg_getVoxelBasedDTIMeasureGradient - (this->floatingImage, - this->warpedImageBw, - this->warpedGradientBw, - this->voxelBasedGradientBw, - this->floatingMask, - this->dtIndicies - ); - break; - case NIFTI_TYPE_FLOAT64: - reg_getVoxelBasedDTIMeasureGradient - (this->floatingImage, - this->warpedImageBw, - this->warpedGradientBw, - this->voxelBasedGradientBw, - this->floatingMask, - this->dtIndicies - ); - break; - default: - reg_print_fct_error("reg_dti::GetVoxelBasedSimilarityMeasureGradient"); - reg_print_msg_error("The input image data type is not supported"); - reg_exit(); - } - } + this->voxelBasedGradientBw->datatype != dtype) { + reg_print_fct_error("reg_dti::GetVoxelBasedSimilarityMeasureGradient"); + reg_print_msg_error("Input images are expected to be of the same type"); + reg_exit(); + } + // Compute the gradient of the nmi for the backward transformation + switch (dtype) { + case NIFTI_TYPE_FLOAT32: + reg_getVoxelBasedDTIMeasureGradient + (this->floatingImage, + this->warpedImageBw, + this->warpedGradientBw, + this->voxelBasedGradientBw, + this->floatingMask, + this->dtIndicies); + break; + case NIFTI_TYPE_FLOAT64: + reg_getVoxelBasedDTIMeasureGradient + (this->floatingImage, + this->warpedImageBw, + this->warpedGradientBw, + this->voxelBasedGradientBw, + this->floatingMask, + this->dtIndicies); + break; + default: + reg_print_fct_error("reg_dti::GetVoxelBasedSimilarityMeasureGradient"); + reg_print_msg_error("The input image data type is not supported"); + reg_exit(); + } + } } /* *************************************************************** */ diff --git a/reg-lib/cpu/_reg_dti.h b/reg-lib/cpu/_reg_dti.h index 580382af..0e6dc21c 100755 --- a/reg-lib/cpu/_reg_dti.h +++ b/reg-lib/cpu/_reg_dti.h @@ -37,8 +37,10 @@ class reg_dti: public reg_measure { nifti_image *warpedImgBw = nullptr, nifti_image *warpedGradBw = nullptr, nifti_image *voxelBasedGradBw = nullptr) override; - /// @brief Returns the value - virtual double GetSimilarityMeasureValue() override; + /// @brief Returns the dti value forwards + virtual double GetSimilarityMeasureValueFw() override; + /// @brief Returns the dti value backwards + virtual double GetSimilarityMeasureValueBw() override; /// @brief Compute the voxel based gradient for DTI images virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override; @@ -57,10 +59,10 @@ class reg_dti: public reg_measure { * @return Returns an L2 measure of the distance between the anisotropic components of the diffusion tensors */ extern "C++" template -double reg_getDTIMeasureValue(nifti_image *referenceImage, - nifti_image *warpedImage, - int *mask, - unsigned *dtIndicies); +double reg_getDTIMeasureValue(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const int *mask, + const unsigned *dtIndicies); /* *************************************************************** */ /** * @brief Compute a voxel based gradient of the sum squared difference. diff --git a/reg-lib/cpu/_reg_kld.cpp b/reg-lib/cpu/_reg_kld.cpp index 39a8b84b..01302e80 100755 --- a/reg-lib/cpu/_reg_kld.cpp +++ b/reg-lib/cpu/_reg_kld.cpp @@ -12,7 +12,6 @@ #include "_reg_kld.h" -/* *************************************************************** */ /* *************************************************************** */ reg_kld::reg_kld(): reg_measure() { #ifndef NDEBUG @@ -20,7 +19,6 @@ reg_kld::reg_kld(): reg_measure() { #endif } /* *************************************************************** */ -/* *************************************************************** */ void reg_kld::InitialiseMeasure(nifti_image *refImg, nifti_image *floImg, int *refMask, @@ -55,11 +53,12 @@ void reg_kld::InitialiseMeasure(nifti_image *refImg, // are meant to be probabilities for (int t = 0; t < this->referenceImage->nt; ++t) { if (this->timePointWeight[t] > 0) { - float min_ref = reg_tools_getMinValue(this->referenceImage, t); - float max_ref = reg_tools_getMaxValue(this->referenceImage, t); - float min_flo = reg_tools_getMinValue(this->floatingImage, t); - float max_flo = reg_tools_getMaxValue(this->floatingImage, t); - if (min_ref < 0.f || min_flo < 0.f || max_ref>1.f || max_flo>1.f) { + const float minRef = reg_tools_getMinValue(this->referenceImage, t); + const float maxRef = reg_tools_getMaxValue(this->referenceImage, t); + const float minFlo = reg_tools_getMinValue(this->floatingImage, t); + const float maxFlo = reg_tools_getMaxValue(this->floatingImage, t); + if (minRef < 0.f || minFlo < 0.f || maxRef > 1.f || maxFlo > 1.f) { + reg_print_fct_error("reg_kld::InitialiseMeasure"); reg_print_msg_error("The input images are expected to be probabilities to use the kld measure"); reg_exit(); } @@ -67,7 +66,7 @@ void reg_kld::InitialiseMeasure(nifti_image *refImg, } #ifndef NDEBUG char text[255]; - reg_print_msg_debug("reg_kld::InitialiseMeasure()."); + reg_print_msg_debug("reg_kld::InitialiseMeasure()"); for (int i = 0; i < this->referenceImage->nt; ++i) { sprintf(text, "Weight for timepoint %i: %f", i, this->timePointWeight[i]); reg_print_msg_debug(text); @@ -75,13 +74,12 @@ void reg_kld::InitialiseMeasure(nifti_image *refImg, #endif } /* *************************************************************** */ -/* *************************************************************** */ template -double reg_getKLDivergence(nifti_image *referenceImage, - nifti_image *warpedImage, - double *timePointWeight, - nifti_image *jacobianDetImg, - int *mask) { +double reg_getKLDivergence(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const double *timePointWeight, + const nifti_image *jacobianDetImg, + const int *mask) { #ifdef _WIN32 long voxel; const long voxelNumber = (long)NiftiImage::calcVoxelNumber(referenceImage, 3); @@ -90,119 +88,77 @@ double reg_getKLDivergence(nifti_image *referenceImage, const size_t voxelNumber = NiftiImage::calcVoxelNumber(referenceImage, 3); #endif - DataType *refPtr = static_cast(referenceImage->data); - DataType *warPtr = static_cast(warpedImage->data); - int *maskPtr = nullptr; - bool MrClean = false; - if (mask == nullptr) { - maskPtr = (int*)calloc(voxelNumber, sizeof(int)); - MrClean = true; - } else maskPtr = &mask[0]; - - DataType *jacPtr = nullptr; + const DataType *refPtr = static_cast(referenceImage->data); + const DataType *warPtr = static_cast(warpedImage->data); + const DataType *jacPtr = nullptr; if (jacobianDetImg != nullptr) jacPtr = static_cast(jacobianDetImg->data); - double measure = 0, measure_tp = 0, num = 0, tempRefValue, tempWarValue, tempValue; + + double measure = 0, measureTp = 0, num = 0, tempRefValue, tempWarValue, tempValue; for (int time = 0; time < referenceImage->nt; ++time) { if (timePointWeight[time] > 0) { - DataType *currentRefPtr = &refPtr[time * voxelNumber]; - DataType *currentWarPtr = &warPtr[time * voxelNumber]; + const DataType *currentRefPtr = &refPtr[time * voxelNumber]; + const DataType *currentWarPtr = &warPtr[time * voxelNumber]; #ifdef _OPENMP #pragma omp parallel for default(none) \ - shared(voxelNumber,currentRefPtr, currentWarPtr, \ - maskPtr, jacobianDetImg, jacPtr) \ + shared(voxelNumber,currentRefPtr, currentWarPtr, mask, jacobianDetImg, jacPtr) \ private(tempRefValue, tempWarValue, tempValue) \ - reduction(+:measure_tp, num) + reduction(+:measureTp, num) #endif for (voxel = 0; voxel < voxelNumber; ++voxel) { - if (maskPtr[voxel] > -1) { + if (mask[voxel] > -1) { tempRefValue = currentRefPtr[voxel] + 1e-16; tempWarValue = currentWarPtr[voxel] + 1e-16; tempValue = tempRefValue * log(tempRefValue / tempWarValue); if (tempValue == tempValue && tempValue != std::numeric_limits::infinity()) { if (jacobianDetImg == nullptr) { - measure_tp -= tempValue; + measureTp -= tempValue; num++; } else { - measure_tp -= tempValue * jacPtr[voxel]; + measureTp -= tempValue * jacPtr[voxel]; num += jacPtr[voxel]; } } } } - measure += measure_tp * timePointWeight[time] / num; + measure += measureTp * timePointWeight[time] / num; } } - if (MrClean) free(maskPtr); return measure; } -template double reg_getKLDivergence(nifti_image*, nifti_image*, double*, nifti_image*, int*); -template double reg_getKLDivergence(nifti_image*, nifti_image*, double*, nifti_image*, int*); /* *************************************************************** */ +double GetSimilarityMeasureValue(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const double *timePointWeight, + const nifti_image *jacobianDetImg, + const int *mask) { + return std::visit([&](auto&& refImgDataType) { + using RefImgDataType = std::decay_t; + return reg_getKLDivergence(referenceImage, + warpedImage, + timePointWeight, + jacobianDetImg, + mask); + }, NiftiImage::getFloatingDataType(referenceImage)); +} /* *************************************************************** */ -double reg_kld::GetSimilarityMeasureValue() { - // Check that all the specified image are of the same datatype - if (this->warpedImage->datatype != this->referenceImage->datatype) { - reg_print_fct_error("reg_kld::GetSimilarityMeasureValue"); - reg_print_msg_error("Both input images are expected to have the same type"); - reg_exit(); - } - double KLDValue; - switch (this->referenceImage->datatype) { - case NIFTI_TYPE_FLOAT32: - KLDValue = reg_getKLDivergence(this->referenceImage, - this->warpedImage, - this->timePointWeight, - nullptr, // TODO this->forwardJacDetImagePointer, - this->referenceMask); - break; - case NIFTI_TYPE_FLOAT64: - KLDValue = reg_getKLDivergence(this->referenceImage, - this->warpedImage, - this->timePointWeight, - nullptr, // TODO this->forwardJacDetImagePointer, - this->referenceMask); - break; - default: - reg_print_fct_error("reg_kld::GetSimilarityMeasureValue"); - reg_print_msg_error("Warped pixel type unsupported"); - reg_exit(); - } - - // Backward computation - if (this->isSymmetric) { - // Check that all the specified image are of the same datatype - if (this->warpedImageBw->datatype != this->floatingImage->datatype) { - reg_print_fct_error("reg_kld::GetSimilarityMeasureValue"); - reg_print_msg_error("Both input images are expected to have the same type"); - reg_exit(); - } - switch (this->floatingImage->datatype) { - case NIFTI_TYPE_FLOAT32: - KLDValue += reg_getKLDivergence(this->floatingImage, - this->warpedImageBw, - this->timePointWeight, - nullptr, // TODO this->backwardJacDetImagePointer, - this->floatingMask); - break; - case NIFTI_TYPE_FLOAT64: - KLDValue += reg_getKLDivergence(this->floatingImage, - this->warpedImageBw, - this->timePointWeight, - nullptr, // TODO this->backwardJacDetImagePointer, - this->floatingMask); - break; - default: - reg_print_fct_error("reg_kld::GetSimilarityMeasureValue"); - reg_print_msg_error("Warped pixel type unsupported"); - reg_exit(); - } - } - return KLDValue; +double reg_kld::GetSimilarityMeasureValueFw() { + return ::GetSimilarityMeasureValue(this->referenceImage, + this->warpedImage, + this->timePointWeight, + nullptr, // TODO this->forwardJacDetImagePointer, + this->referenceMask); } /* *************************************************************** */ +double reg_kld::GetSimilarityMeasureValueBw() { + return ::GetSimilarityMeasureValue(this->floatingImage, + this->warpedImageBw, + this->timePointWeight, + nullptr, // TODO this->backwardJacDetImagePointer, + this->floatingMask); +} /* *************************************************************** */ template void reg_getKLDivergenceVoxelBasedGradient(nifti_image *referenceImage, @@ -313,11 +269,6 @@ void reg_getKLDivergenceVoxelBasedGradient(nifti_image *referenceImage, } if (MrClean) free(maskPtr); } -template void reg_getKLDivergenceVoxelBasedGradient -(nifti_image*, nifti_image*, nifti_image*, nifti_image*, nifti_image*, int*, int, double); -template void reg_getKLDivergenceVoxelBasedGradient -(nifti_image*, nifti_image*, nifti_image*, nifti_image*, nifti_image*, int*, int, double); -/* *************************************************************** */ /* *************************************************************** */ void reg_kld::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { // Check if the specified time point exists and is active @@ -401,4 +352,3 @@ void reg_kld::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { } } /* *************************************************************** */ -/* *************************************************************** */ diff --git a/reg-lib/cpu/_reg_kld.h b/reg-lib/cpu/_reg_kld.h index aaf70556..ae5f4cb2 100755 --- a/reg-lib/cpu/_reg_kld.h +++ b/reg-lib/cpu/_reg_kld.h @@ -34,8 +34,10 @@ class reg_kld: public reg_measure { nifti_image *warpedImgBw = nullptr, nifti_image *warpedGradBw = nullptr, nifti_image *voxelBasedGradBw = nullptr) override; - /// @brief Returns the kld value - virtual double GetSimilarityMeasureValue() override; + /// @brief Returns the kld value forwards + virtual double GetSimilarityMeasureValueFw() override; + /// @brief Returns the kld value backwards + virtual double GetSimilarityMeasureValueBw() override; /// @brief Compute the voxel based kld gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override; }; @@ -50,15 +52,15 @@ class reg_kld: public reg_measure { * image is used to modulate the KLD. The argument is ignored if the * pointer is set to nullptr * @param mask Array that contains a mask to specify which voxel - * should be considered. If set to nullptr, all voxels are considered + * should be considered * @return Returns the computed sum squared difference */ extern "C++" template -double reg_getKLDivergence(nifti_image *reference, - nifti_image *warped, - double *timePointWeight, - nifti_image *jacobianDeterminantImage, - int *mask); +double reg_getKLDivergence(const nifti_image *reference, + const nifti_image *warped, + const double *timePointWeight, + const nifti_image *jacobianDeterminantImage, + const int *mask); /* *************************************************************** */ /** @brief Compute a voxel based gradient of the sum squared difference. diff --git a/reg-lib/cpu/_reg_lncc.cpp b/reg-lib/cpu/_reg_lncc.cpp index fca452e3..2d1c3848 100644 --- a/reg-lib/cpu/_reg_lncc.cpp +++ b/reg-lib/cpu/_reg_lncc.cpp @@ -79,16 +79,18 @@ reg_lncc::~reg_lncc() { } /* *************************************************************** */ template -void reg_lncc::UpdateLocalStatImages(nifti_image *refImage, - nifti_image *warImage, - nifti_image *meanImage, - nifti_image *warpedMeanImage, - nifti_image *stdDevImage, - nifti_image *warpedSdevImage, - int *refMask, - int *combinedMask, - int currentTimepoint) { - // Generate the forward mask to ignore all NaN values +void UpdateLocalStatImages(const nifti_image *refImage, + const nifti_image *warImage, + nifti_image *meanImage, + nifti_image *warpedMeanImage, + nifti_image *sdevImage, + nifti_image *warpedSdevImage, + const int *refMask, + int *combinedMask, + const float *kernelStandardDeviation, + const int& kernelType, + const int& currentTimepoint) { + // Generate the combined mask to ignore all NaN values #ifdef _WIN32 long voxel; const long voxelNumber = (long)NiftiImage::calcVoxelNumber(refImage, 3); @@ -100,25 +102,25 @@ void reg_lncc::UpdateLocalStatImages(nifti_image *refImage, reg_tools_removeNanFromMask(refImage, combinedMask); reg_tools_removeNanFromMask(warImage, combinedMask); - DataType *origRefPtr = static_cast(refImage->data); + const DataType *origRefPtr = static_cast(refImage->data); DataType *meanImgPtr = static_cast(meanImage->data); - DataType *sdevImgPtr = static_cast(stdDevImage->data); + DataType *sdevImgPtr = static_cast(sdevImage->data); memcpy(meanImgPtr, &origRefPtr[currentTimepoint * voxelNumber], voxelNumber * refImage->nbyper); memcpy(sdevImgPtr, &origRefPtr[currentTimepoint * voxelNumber], voxelNumber * refImage->nbyper); - reg_tools_multiplyImageToImage(stdDevImage, stdDevImage, stdDevImage); - reg_tools_kernelConvolution(meanImage, this->kernelStandardDeviation, this->kernelType, combinedMask); - reg_tools_kernelConvolution(stdDevImage, this->kernelStandardDeviation, this->kernelType, combinedMask); + reg_tools_multiplyImageToImage(sdevImage, sdevImage, sdevImage); + reg_tools_kernelConvolution(meanImage, kernelStandardDeviation, kernelType, combinedMask); + reg_tools_kernelConvolution(sdevImage, kernelStandardDeviation, kernelType, combinedMask); - DataType *origWarPtr = static_cast(warImage->data); + const DataType *origWarPtr = static_cast(warImage->data); DataType *warMeanPtr = static_cast(warpedMeanImage->data); DataType *warSdevPtr = static_cast(warpedSdevImage->data); memcpy(warMeanPtr, &origWarPtr[currentTimepoint * voxelNumber], voxelNumber * warImage->nbyper); memcpy(warSdevPtr, &origWarPtr[currentTimepoint * voxelNumber], voxelNumber * warImage->nbyper); reg_tools_multiplyImageToImage(warpedSdevImage, warpedSdevImage, warpedSdevImage); - reg_tools_kernelConvolution(warpedMeanImage, this->kernelStandardDeviation, this->kernelType, combinedMask); - reg_tools_kernelConvolution(warpedSdevImage, this->kernelStandardDeviation, this->kernelType, combinedMask); + reg_tools_kernelConvolution(warpedMeanImage, kernelStandardDeviation, kernelType, combinedMask); + reg_tools_kernelConvolution(warpedSdevImage, kernelStandardDeviation, kernelType, combinedMask); #ifdef _OPENMP #pragma omp parallel for default(none) \ shared(voxelNumber, sdevImgPtr, meanImgPtr, warSdevPtr, warMeanPtr) @@ -243,7 +245,7 @@ void reg_lncc::InitialiseMeasure(nifti_image *refImg, } #ifndef NDEBUG char text[255]; - reg_print_msg_debug("reg_lncc::InitialiseMeasure()."); + reg_print_msg_debug("reg_lncc::InitialiseMeasure()"); for (int i = 0; i < this->referenceImage->nt; ++i) { sprintf(text, "Weight for timepoint %i: %f", i, this->timePointWeight[i]); reg_print_msg_debug(text); @@ -252,17 +254,17 @@ void reg_lncc::InitialiseMeasure(nifti_image *refImg, } /* *************************************************************** */ template -double reg_getLNCCValue(nifti_image *referenceImage, - nifti_image *meanImage, - nifti_image *sdevImage, - nifti_image *warpedImage, - nifti_image *warpedMeanImage, - nifti_image *warpedSdevImage, - int *combinedMask, - float *kernelStandardDeviation, +double reg_getLnccValue(const nifti_image *referenceImage, + const nifti_image *meanImage, + const nifti_image *sdevImage, + const nifti_image *warpedImage, + const nifti_image *warpedMeanImage, + const nifti_image *warpedSdevImage, + const int *combinedMask, + const float *kernelStandardDeviation, nifti_image *correlationImage, - int kernelType, - int currentTimepoint) { + const int& kernelType, + const int& currentTimepoint) { #ifdef _WIN32 long voxel; const long voxelNumber = (long)NiftiImage::calcVoxelNumber(referenceImage, 3); @@ -272,16 +274,16 @@ double reg_getLNCCValue(nifti_image *referenceImage, #endif // Compute the local correlation - DataType *refImagePtr = static_cast(referenceImage->data); - DataType *currentRefPtr = &refImagePtr[currentTimepoint * voxelNumber]; + const DataType *refImagePtr = static_cast(referenceImage->data); + const DataType *currentRefPtr = &refImagePtr[currentTimepoint * voxelNumber]; - DataType *warImagePtr = static_cast(warpedImage->data); - DataType *currentWarPtr = &warImagePtr[currentTimepoint * voxelNumber]; + const DataType *warImagePtr = static_cast(warpedImage->data); + const DataType *currentWarPtr = &warImagePtr[currentTimepoint * voxelNumber]; - DataType *meanImgPtr = static_cast(meanImage->data); - DataType *warMeanPtr = static_cast(warpedMeanImage->data); - DataType *sdevImgPtr = static_cast(sdevImage->data); - DataType *warSdevPtr = static_cast(warpedSdevImage->data); + const DataType *meanImgPtr = static_cast(meanImage->data); + const DataType *warMeanPtr = static_cast(warpedMeanImage->data); + const DataType *sdevImgPtr = static_cast(sdevImage->data); + const DataType *warSdevPtr = static_cast(warpedSdevImage->data); DataType *correlationPtr = static_cast(correlationImage->data); for (size_t i = 0; i < voxelNumber; ++i) @@ -289,156 +291,113 @@ double reg_getLNCCValue(nifti_image *referenceImage, reg_tools_kernelConvolution(correlationImage, kernelStandardDeviation, kernelType, combinedMask); - double lncc_value_sum = 0., lncc_value; - double activeVoxel_num = 0.; + double lnccSum = 0, lncc; + size_t activeVoxelNumber = 0; // Iteration over all voxels #ifdef _OPENMP #pragma omp parallel for default(none) \ shared(voxelNumber,combinedMask,meanImgPtr,warMeanPtr, \ sdevImgPtr,warSdevPtr,correlationPtr) \ - private(lncc_value) \ - reduction(+:lncc_value_sum) \ - reduction(+:activeVoxel_num) + private(lncc) \ + reduction(+:lnccSum, activeVoxelNumber) #endif for (voxel = 0; voxel < voxelNumber; ++voxel) { // Check if the current voxel belongs to the mask if (combinedMask[voxel] > -1) { - lncc_value = (correlationPtr[voxel] - (meanImgPtr[voxel] * warMeanPtr[voxel])) / (sdevImgPtr[voxel] * warSdevPtr[voxel]); - if (lncc_value == lncc_value && isinf(lncc_value) == 0) { - lncc_value_sum += fabs(lncc_value); - ++activeVoxel_num; + lncc = (correlationPtr[voxel] - (meanImgPtr[voxel] * warMeanPtr[voxel])) / (sdevImgPtr[voxel] * warSdevPtr[voxel]); + if (lncc == lncc && !isinf(lncc)) { + lnccSum += fabs(lncc); + ++activeVoxelNumber; } } } - return lncc_value_sum / activeVoxel_num; + return lnccSum / activeVoxelNumber; } /* *************************************************************** */ -double reg_lncc::GetSimilarityMeasureValue() { - double lncc_value = 0; - - for (int currentTimepoint = 0; currentTimepoint < this->referenceImage->nt; ++currentTimepoint) { - if (this->timePointWeight[currentTimepoint] > 0) { - double tp_value = 0; - // Compute the mean and variance of the reference and warped floating - switch (this->referenceImage->datatype) { - case NIFTI_TYPE_FLOAT32: - this->UpdateLocalStatImages(this->referenceImage, - this->warpedImage, - this->meanImage, - this->warpedMeanImage, - this->sdevImage, - this->warpedSdevImage, - this->referenceMask, - this->forwardMask, - currentTimepoint); - break; - case NIFTI_TYPE_FLOAT64: - this->UpdateLocalStatImages(this->referenceImage, - this->warpedImage, - this->meanImage, - this->warpedMeanImage, - this->sdevImage, - this->warpedSdevImage, - this->referenceMask, - this->forwardMask, - currentTimepoint); - break; - } - - // Compute the LNCC - Forward - switch (this->referenceImage->datatype) { - case NIFTI_TYPE_FLOAT32: - tp_value += reg_getLNCCValue(this->referenceImage, - this->meanImage, - this->sdevImage, - this->warpedImage, - this->warpedMeanImage, - this->warpedSdevImage, - this->forwardMask, - this->kernelStandardDeviation, - this->correlationImage, - this->kernelType, - currentTimepoint); - break; - case NIFTI_TYPE_FLOAT64: - tp_value += reg_getLNCCValue(this->referenceImage, - this->meanImage, - this->sdevImage, - this->warpedImage, - this->warpedMeanImage, - this->warpedSdevImage, - this->forwardMask, - this->kernelStandardDeviation, - this->correlationImage, - this->kernelType, - currentTimepoint); - break; - } - if (this->isSymmetric) { - // Compute the mean and variance of the floating and warped reference - switch (this->floatingImage->datatype) { - case NIFTI_TYPE_FLOAT32: - this->UpdateLocalStatImages(this->floatingImage, - this->warpedImageBw, - this->meanImageBw, - this->warpedMeanImageBw, - this->sdevImageBw, - this->warpedSdevImageBw, - this->floatingMask, - this->backwardMask, - currentTimepoint); - break; - case NIFTI_TYPE_FLOAT64: - this->UpdateLocalStatImages(this->floatingImage, - this->warpedImageBw, - this->meanImageBw, - this->warpedMeanImageBw, - this->sdevImageBw, - this->warpedSdevImageBw, - this->floatingMask, - this->backwardMask, - currentTimepoint); - break; - } - // Compute the LNCC - Backward - switch (this->floatingImage->datatype) { - case NIFTI_TYPE_FLOAT32: - tp_value += reg_getLNCCValue(this->floatingImage, - this->meanImageBw, - this->sdevImageBw, - this->warpedImageBw, - this->warpedMeanImageBw, - this->warpedSdevImageBw, - this->backwardMask, - this->kernelStandardDeviation, - this->correlationImageBw, - this->kernelType, +double GetSimilarityMeasureValue(const nifti_image *referenceImage, + nifti_image *meanImage, + nifti_image *sdevImage, + const nifti_image *warpedImage, + nifti_image *warpedMeanImage, + nifti_image *warpedSdevImage, + const int *referenceMask, + int *combinedMask, + const float *kernelStandardDeviation, + nifti_image *correlationImage, + const int& kernelType, + const int& referenceTimePoint, + const double *timePointWeight) { + double lncc = 0; + for (int currentTimepoint = 0; currentTimepoint < referenceTimePoint; ++currentTimepoint) { + if (timePointWeight[currentTimepoint] > 0) { + const double tp = std::visit([&](auto&& refImgDataType) { + using RefImgDataType = std::decay_t; + // Compute the mean and variance of the reference and warped floating + UpdateLocalStatImages(referenceImage, + warpedImage, + meanImage, + warpedMeanImage, + sdevImage, + warpedSdevImage, + referenceMask, + combinedMask, + kernelStandardDeviation, + kernelType, + currentTimepoint); + // Compute the LNCC value + return reg_getLnccValue(referenceImage, + meanImage, + sdevImage, + warpedImage, + warpedMeanImage, + warpedSdevImage, + combinedMask, + kernelStandardDeviation, + correlationImage, + kernelType, currentTimepoint); - break; - case NIFTI_TYPE_FLOAT64: - tp_value += reg_getLNCCValue(this->floatingImage, - this->meanImageBw, - this->sdevImageBw, - this->warpedImageBw, - this->warpedMeanImageBw, - this->warpedSdevImageBw, - this->backwardMask, - this->kernelStandardDeviation, - this->correlationImageBw, - this->kernelType, - currentTimepoint); - break; - } - } - lncc_value += tp_value * this->timePointWeight[currentTimepoint]; + }, NiftiImage::getFloatingDataType(referenceImage)); + lncc += tp * timePointWeight[currentTimepoint]; } } - return lncc_value; + return lncc; +} +/* *************************************************************** */ +double reg_lncc::GetSimilarityMeasureValueFw() { + return ::GetSimilarityMeasureValue(this->referenceImage, + this->meanImage, + this->sdevImage, + this->warpedImage, + this->warpedMeanImage, + this->warpedSdevImage, + this->referenceMask, + this->forwardMask, + this->kernelStandardDeviation, + this->correlationImage, + this->kernelType, + this->referenceTimePoint, + this->timePointWeight); +} +/* *************************************************************** */ +double reg_lncc::GetSimilarityMeasureValueBw() { + return ::GetSimilarityMeasureValue(this->floatingImage, + this->meanImageBw, + this->sdevImageBw, + this->warpedImageBw, + this->warpedMeanImageBw, + this->warpedSdevImageBw, + this->floatingMask, + this->backwardMask, + this->kernelStandardDeviation, + this->correlationImageBw, + this->kernelType, + this->referenceTimePoint, + this->timePointWeight); } /* *************************************************************** */ template -void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage, +void reg_getVoxelBasedLnccGradient(nifti_image *referenceImage, nifti_image *meanImage, nifti_image *sdevImage, nifti_image *warpedImage, @@ -480,7 +439,7 @@ void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage, double refMeanValue, warMeanValue, refSdevValue, warSdevValue, correlaValue; double temp1, temp2, temp3; - double activeVoxel_num = 0; + size_t activeVoxelNumber = 0; // Iteration over all voxels #ifdef _OPENMP @@ -489,12 +448,11 @@ void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage, sdevImgPtr,warSdevPtr,correlationPtr) \ private(refMeanValue,warMeanValue,refSdevValue, \ warSdevValue, correlaValue, temp1, temp2, temp3) \ - reduction(+:activeVoxel_num) + reduction(+:activeVoxelNumber) #endif for (voxel = 0; voxel < voxelNumber; ++voxel) { // Check if the current voxel belongs to the mask if (combinedMask[voxel] > -1) { - refMeanValue = meanImgPtr[voxel]; warMeanValue = warMeanPtr[voxel]; refSdevValue = sdevImgPtr[voxel]; @@ -502,8 +460,7 @@ void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage, correlaValue = correlationPtr[voxel] - (refMeanValue * warMeanValue); temp1 = 1.0 / (refSdevValue * warSdevValue); - temp2 = correlaValue / - (refSdevValue * warSdevValue * warSdevValue * warSdevValue); + temp2 = correlaValue / (refSdevValue * warSdevValue * warSdevValue * warSdevValue); temp3 = (correlaValue * warMeanValue) / (refSdevValue * warSdevValue * warSdevValue * warSdevValue) - @@ -520,13 +477,13 @@ void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage, warMeanPtr[voxel] = static_cast(temp1); warSdevPtr[voxel] = static_cast(temp2); correlationPtr[voxel] = static_cast(temp3); - activeVoxel_num++; + activeVoxelNumber++; } else warMeanPtr[voxel] = warSdevPtr[voxel] = correlationPtr[voxel] = 0; } else warMeanPtr[voxel] = warSdevPtr[voxel] = correlationPtr[voxel] = 0; } //adjust weight for number of voxels - double adjusted_weight = timepointWeight / activeVoxel_num; + double adjusted_weight = timepointWeight / activeVoxelNumber; // Smooth the newly computed values reg_tools_kernelConvolution(warpedMeanImage, kernelStandardDeviation, kernelType, combinedMask); @@ -593,33 +550,37 @@ void reg_lncc::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { // Compute the mean and variance of the reference and warped floating switch (this->referenceImage->datatype) { case NIFTI_TYPE_FLOAT32: - this->UpdateLocalStatImages(this->referenceImage, - this->warpedImage, - this->meanImage, - this->warpedMeanImage, - this->sdevImage, - this->warpedSdevImage, - this->referenceMask, - this->forwardMask, - currentTimepoint); + UpdateLocalStatImages(this->referenceImage, + this->warpedImage, + this->meanImage, + this->warpedMeanImage, + this->sdevImage, + this->warpedSdevImage, + this->referenceMask, + this->forwardMask, + this->kernelStandardDeviation, + this->kernelType, + currentTimepoint); break; case NIFTI_TYPE_FLOAT64: - this->UpdateLocalStatImages(this->referenceImage, - this->warpedImage, - this->meanImage, - this->warpedMeanImage, - this->sdevImage, - this->warpedSdevImage, - this->referenceMask, - this->forwardMask, - currentTimepoint); + UpdateLocalStatImages(this->referenceImage, + this->warpedImage, + this->meanImage, + this->warpedMeanImage, + this->sdevImage, + this->warpedSdevImage, + this->referenceMask, + this->forwardMask, + this->kernelStandardDeviation, + this->kernelType, + currentTimepoint); break; } // Compute the LNCC gradient - Forward switch (this->referenceImage->datatype) { case NIFTI_TYPE_FLOAT32: - reg_getVoxelBasedLNCCGradient(this->referenceImage, + reg_getVoxelBasedLnccGradient(this->referenceImage, this->meanImage, this->sdevImage, this->warpedImage, @@ -635,7 +596,7 @@ void reg_lncc::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { this->timePointWeight[currentTimepoint]); break; case NIFTI_TYPE_FLOAT64: - reg_getVoxelBasedLNCCGradient(this->referenceImage, + reg_getVoxelBasedLnccGradient(this->referenceImage, this->meanImage, this->sdevImage, this->warpedImage, @@ -655,32 +616,36 @@ void reg_lncc::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { // Compute the mean and variance of the floating and warped reference switch (this->floatingImage->datatype) { case NIFTI_TYPE_FLOAT32: - this->UpdateLocalStatImages(this->floatingImage, - this->warpedImageBw, - this->meanImageBw, - this->warpedMeanImageBw, - this->sdevImageBw, - this->warpedSdevImageBw, - this->floatingMask, - this->backwardMask, - currentTimepoint); + UpdateLocalStatImages(this->floatingImage, + this->warpedImageBw, + this->meanImageBw, + this->warpedMeanImageBw, + this->sdevImageBw, + this->warpedSdevImageBw, + this->floatingMask, + this->backwardMask, + this->kernelStandardDeviation, + this->kernelType, + currentTimepoint); break; case NIFTI_TYPE_FLOAT64: - this->UpdateLocalStatImages(this->floatingImage, - this->warpedImageBw, - this->meanImageBw, - this->warpedMeanImageBw, - this->sdevImageBw, - this->warpedSdevImageBw, - this->floatingMask, - this->backwardMask, - currentTimepoint); + UpdateLocalStatImages(this->floatingImage, + this->warpedImageBw, + this->meanImageBw, + this->warpedMeanImageBw, + this->sdevImageBw, + this->warpedSdevImageBw, + this->floatingMask, + this->backwardMask, + this->kernelStandardDeviation, + this->kernelType, + currentTimepoint); break; } // Compute the LNCC gradient - Backward switch (this->floatingImage->datatype) { case NIFTI_TYPE_FLOAT32: - reg_getVoxelBasedLNCCGradient(this->floatingImage, + reg_getVoxelBasedLnccGradient(this->floatingImage, this->meanImageBw, this->sdevImageBw, this->warpedImageBw, @@ -696,7 +661,7 @@ void reg_lncc::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { this->timePointWeight[currentTimepoint]); break; case NIFTI_TYPE_FLOAT64: - reg_getVoxelBasedLNCCGradient(this->floatingImage, + reg_getVoxelBasedLnccGradient(this->floatingImage, this->meanImageBw, this->sdevImageBw, this->warpedImageBw, diff --git a/reg-lib/cpu/_reg_lncc.h b/reg-lib/cpu/_reg_lncc.h index 5a7b5ef0..6c7dda5a 100644 --- a/reg-lib/cpu/_reg_lncc.h +++ b/reg-lib/cpu/_reg_lncc.h @@ -34,15 +34,17 @@ class reg_lncc: public reg_measure { nifti_image *warpedImgBw = nullptr, nifti_image *warpedGradBw = nullptr, nifti_image *voxelBasedGradBw = nullptr) override; - /// @brief Returns the lncc value - virtual double GetSimilarityMeasureValue() override; + /// @brief Returns the lncc value forwards + virtual double GetSimilarityMeasureValueFw() override; + /// @brief Returns the lncc value backwards + virtual double GetSimilarityMeasureValueBw() override; /// @brief Compute the voxel based lncc gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override; - /// @brief Stuff + /// @brief Set the kernel standard deviation virtual void SetKernelStandardDeviation(int t, float stddev) { this->kernelStandardDeviation[t] = stddev; } - /// @brief Stuff + /// @brief Set the kernel type virtual void SetKernelType(int t) { this->kernelType = t; } @@ -64,65 +66,5 @@ class reg_lncc: public reg_measure { int *backwardMask; int kernelType; - - template - void UpdateLocalStatImages(nifti_image *refImage, - nifti_image *warImage, - nifti_image *meanImage, - nifti_image *warpedMeanImage, - nifti_image *stdDevImage, - nifti_image *warpedSdevImage, - int *refMask, - int *mask, - int currentTimepoint); }; /* *************************************************************** */ -/** @brief Compute and return the LNCC between two input image - * @param referenceImage First input image to use to compute the metric - * @param warpedImage Second input image to use to compute the metric - * @param gaussianStandardDeviation Standard deviation of the Gaussian kernel - * to use. - * @param mask Array that contains a mask to specify which voxel - * should be considered. If set to nullptr, all voxels are considered - * @return Returns the computed LNCC - */ -extern "C++" template -double reg_getLNCCValue(nifti_image *referenceImage, - nifti_image *meanImage, - nifti_image *sdevImage, - nifti_image *warpedImage, - nifti_image *warpedMeanImage, - nifti_image *warpedSdevImage, - int *combinedMask, - float *kernelStandardDeviation, - nifti_image *correlationImage, - int kernelType, - int currentTimepoint); -/* *************************************************************** */ -/** @brief Compute a voxel based gradient of the LNCC. - * @param referenceImage First input image to use to compute the metric - * @param warpedImage Second input image to use to compute the metric - * @param warpedImageGradient Spatial gradient of the input warped image - * @param lnccGradientImage Output image that will be updated with the - * value of the LNCC gradient - * @param gaussianStandardDeviation Standard deviation of the Gaussian kernel - * to use. - * @param mask Array that contains a mask to specify which voxel - * should be considered. If set to nullptr, all voxels are considered - */ -extern "C++" template -void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage, - nifti_image *meanImage, - nifti_image *sdevImage, - nifti_image *warpedImage, - nifti_image *warpedMeanImage, - nifti_image *warpedStdDevImage, - int *combinedMask, - float *kernelStdDev, - nifti_image *correlationImage, - nifti_image *warpedGradient, - nifti_image *lnccGradientImage, - int kernelType, - int currentTimepoint, - double timepointWeight); -/* *************************************************************** */ diff --git a/reg-lib/cpu/_reg_measure.h b/reg-lib/cpu/_reg_measure.h index ee2a2625..56c42d50 100755 --- a/reg-lib/cpu/_reg_measure.h +++ b/reg-lib/cpu/_reg_measure.h @@ -16,7 +16,7 @@ class reg_measure { /// @brief Measure class constructor reg_measure() { #ifndef NDEBUG - printf("[NiftyReg DEBUG] reg_measure constructor called\n"); + reg_print_msg_debug("reg_measure constructor called"); #endif } /// @brief Measure class destructor @@ -56,12 +56,47 @@ class reg_measure { this->voxelBasedGradientBw = nullptr; } #ifndef NDEBUG - printf("[NiftyReg DEBUG] reg_measure::InitialiseMeasure()\n"); + reg_print_msg_debug("reg_measure::InitialiseMeasure()"); #endif } + /// @brief Returns the forward registration measure of similarity value + virtual double GetSimilarityMeasureValueFw() = 0; + /// @brief Returns the backward registration measure of similarity value + virtual double GetSimilarityMeasureValueBw() = 0; /// @brief Returns the registration measure of similarity value - virtual double GetSimilarityMeasureValue() = 0; + double GetSimilarityMeasureValue() { // Do not override + // Check that all the specified image are of the same datatype + if (this->referenceImage->datatype != NIFTI_TYPE_FLOAT32 && this->referenceImage->datatype != NIFTI_TYPE_FLOAT64) { + reg_print_fct_error("reg_measure::GetSimilarityMeasureValue()"); + reg_print_msg_error("Input images are expected to be of floating precision type"); + reg_exit(); + } + if (this->warpedImage->datatype != this->referenceImage->datatype) { + reg_print_fct_error("reg_measure::GetSimilarityMeasureValue()"); + reg_print_msg_error("Both input images are expected to have the same type"); + reg_exit(); + } + double sim = GetSimilarityMeasureValueFw(); + if (this->isSymmetric) { + // Check that all the specified image are of the same datatype + if (this->floatingImage->datatype != NIFTI_TYPE_FLOAT32 && this->floatingImage->datatype != NIFTI_TYPE_FLOAT64) { + reg_print_fct_error("reg_measure::GetSimilarityMeasureValue()"); + reg_print_msg_error("Input images are expected to be of floating precision type"); + reg_exit(); + } + if (this->floatingImage->datatype != this->warpedImageBw->datatype) { + reg_print_fct_error("reg_measure::GetSimilarityMeasureValue()"); + reg_print_msg_error("Both input images are expected to have the same type"); + reg_exit(); + } + sim += GetSimilarityMeasureValueBw(); + } +#ifndef NDEBUG + reg_print_msg_debug("reg_measure::GetSimilarityMeasureValue called"); +#endif + return sim; + } /// @brief Compute the voxel based measure of similarity gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { diff --git a/reg-lib/cpu/_reg_mind.cpp b/reg-lib/cpu/_reg_mind.cpp index 29aa32c9..7b289c27 100644 --- a/reg-lib/cpu/_reg_mind.cpp +++ b/reg-lib/cpu/_reg_mind.cpp @@ -14,14 +14,14 @@ /* *************************************************************** */ template -void ShiftImage(nifti_image* inputImgPtr, - nifti_image* shiftedImgPtr, - int *maskPtr, - int tx, - int ty, - int tz) { - DataType* inputData = static_cast(inputImgPtr->data); - DataType* shiftImageData = static_cast(shiftedImgPtr->data); +void ShiftImage(const nifti_image *inputImage, + nifti_image *shiftedImage, + const int *mask, + const int& tx, + const int& ty, + const int& tz) { + const DataType* inputData = static_cast(inputImage->data); + DataType* shiftImageData = static_cast(shiftedImage->data); int currentIndex; int shiftedIndex; @@ -30,23 +30,21 @@ void ShiftImage(nifti_image* inputImgPtr, #ifdef _OPENMP #pragma omp parallel for default(none) \ - shared(inputData, shiftImageData, shiftedImgPtr, inputImgPtr, \ - maskPtr, tx, ty, tz) \ - private(x, y, old_x, old_y, old_z, shiftedIndex, \ - currentIndex) + shared(inputData, shiftImageData, shiftedImage, inputImage, mask, tx, ty, tz) \ + private(x, y, old_x, old_y, old_z, shiftedIndex, currentIndex) #endif - for (z = 0; z < shiftedImgPtr->nz; z++) { - currentIndex = z * shiftedImgPtr->nx * shiftedImgPtr->ny; + for (z = 0; z < shiftedImage->nz; z++) { + currentIndex = z * shiftedImage->nx * shiftedImage->ny; old_z = z - tz; - for (y = 0; y < shiftedImgPtr->ny; y++) { + for (y = 0; y < shiftedImage->ny; y++) { old_y = y - ty; - for (x = 0; x < shiftedImgPtr->nx; x++) { + for (x = 0; x < shiftedImage->nx; x++) { old_x = x - tx; - if (old_x > -1 && old_xnx && - old_y>-1 && old_yny && - old_z>-1 && old_z < inputImgPtr->nz) { - shiftedIndex = (old_z * inputImgPtr->ny + old_y) * inputImgPtr->nx + old_x; - if (maskPtr[shiftedIndex] > -1) { + if (old_x > -1 && old_x < inputImage->nx && + old_y > -1 && old_y < inputImage->ny && + old_z > -1 && old_z < inputImage->nz) { + shiftedIndex = (old_z * inputImage->ny + old_y) * inputImage->nx + old_x; + if (mask[shiftedIndex] > -1) { shiftImageData[currentIndex] = inputData[shiftedIndex]; } // mask is not defined else { @@ -65,11 +63,11 @@ void ShiftImage(nifti_image* inputImgPtr, } /* *************************************************************** */ template -void GetMINDImageDescriptor_core(nifti_image* inputImage, - nifti_image* MINDImage, - int *maskPtr, - int descriptorOffset, - int currentTimepoint) { +void GetMindImageDescriptorCore(const nifti_image *inputImage, + nifti_image *mindImage, + const int *mask, + const int& descriptorOffset, + const int& currentTimepoint) { #ifdef WIN32 long voxelIndex; const long voxelNumber = (long)NiftiImage::calcVoxelNumber(inputImage, 3); @@ -79,7 +77,7 @@ void GetMINDImageDescriptor_core(nifti_image* inputImage, #endif // Create a pointer to the descriptor image - DataType* MINDImgDataPtr = static_cast(MINDImage->data); + DataType* mindImgDataPtr = static_cast(mindImage->data); // Allocate an image to store the current timepoint reference image nifti_image *currentInputImage = nifti_copy_nim_info(inputImage); @@ -87,7 +85,7 @@ void GetMINDImageDescriptor_core(nifti_image* inputImage, currentInputImage->nt = currentInputImage->dim[4] = 1; currentInputImage->nvox = voxelNumber; DataType *inputImagePtr = static_cast(inputImage->data); - currentInputImage->data = static_cast(&inputImagePtr[currentTimepoint * voxelNumber]); + currentInputImage->data = &inputImagePtr[currentTimepoint * voxelNumber]; // Allocate an image to store the mean image nifti_image *meanImage = nifti_dup(*currentInputImage, false); @@ -97,96 +95,95 @@ void GetMINDImageDescriptor_core(nifti_image* inputImage, nifti_image *shiftedImage = nifti_dup(*currentInputImage, false); // Allocation of the difference image - nifti_image *diff_image = nifti_dup(*currentInputImage, false); + nifti_image *diffImage = nifti_dup(*currentInputImage, false); // Define the sigma for the convolution float sigma = -0.5;// negative value denotes voxel width //2D version int samplingNbr = (currentInputImage->nz > 1) ? 6 : 4; - int RSampling3D_x[6] = { -descriptorOffset, descriptorOffset, 0, 0, 0, 0 }; - int RSampling3D_y[6] = { 0, 0, -descriptorOffset, descriptorOffset, 0, 0 }; - int RSampling3D_z[6] = { 0, 0, 0, 0, -descriptorOffset, descriptorOffset }; + int rSamplingX[6] = { -descriptorOffset, descriptorOffset, 0, 0, 0, 0 }; + int rSamplingY[6] = { 0, 0, -descriptorOffset, descriptorOffset, 0, 0 }; + int rSamplingZ[6] = { 0, 0, 0, 0, -descriptorOffset, descriptorOffset }; for (int i = 0; i < samplingNbr; i++) { - ShiftImage(currentInputImage, shiftedImage, maskPtr, - RSampling3D_x[i], RSampling3D_y[i], RSampling3D_z[i]); - reg_tools_subtractImageFromImage(currentInputImage, shiftedImage, diff_image); - reg_tools_multiplyImageToImage(diff_image, diff_image, diff_image); - reg_tools_kernelConvolution(diff_image, &sigma, GAUSSIAN_KERNEL, maskPtr); - reg_tools_addImageToImage(meanImage, diff_image, meanImage); + ShiftImage(currentInputImage, shiftedImage, mask, rSamplingX[i], rSamplingY[i], rSamplingZ[i]); + reg_tools_subtractImageFromImage(currentInputImage, shiftedImage, diffImage); + reg_tools_multiplyImageToImage(diffImage, diffImage, diffImage); + reg_tools_kernelConvolution(diffImage, &sigma, GAUSSIAN_KERNEL, mask); + reg_tools_addImageToImage(meanImage, diffImage, meanImage); // Store the current descriptor - const size_t index = i * diff_image->nvox; - memcpy(&MINDImgDataPtr[index], diff_image->data, diff_image->nbyper * diff_image->nvox); + const size_t index = i * diffImage->nvox; + memcpy(&mindImgDataPtr[index], diffImage->data, diffImage->nbyper * diffImage->nvox); } // Compute the mean over the number of sample reg_tools_divideValueToImage(meanImage, meanImage, samplingNbr); // Compute the MIND descriptor int mindIndex; - DataType meanValue, max_desc, descValue; + DataType meanValue, maxDesc, descValue; #ifdef _OPENMP #pragma omp parallel for default(none) \ - shared(voxelNumber, samplingNbr, maskPtr, meanImgDataPtr, \ - MINDImgDataPtr) \ - private(meanValue, max_desc, descValue, mindIndex) + shared(voxelNumber, samplingNbr, mask, meanImgDataPtr, \ + mindImgDataPtr) \ + private(meanValue, maxDesc, descValue, mindIndex) #endif for (voxelIndex = 0; voxelIndex < voxelNumber; voxelIndex++) { - if (maskPtr[voxelIndex] > -1) { + if (mask[voxelIndex] > -1) { // Get the mean value for the current voxel meanValue = meanImgDataPtr[voxelIndex]; if (meanValue == 0) { meanValue = std::numeric_limits::epsilon(); } - max_desc = 0; + maxDesc = 0; mindIndex = voxelIndex; for (int t = 0; t < samplingNbr; t++) { - descValue = (DataType)exp(-MINDImgDataPtr[mindIndex] / meanValue); - MINDImgDataPtr[mindIndex] = descValue; - max_desc = (std::max)(max_desc, descValue); + descValue = (DataType)exp(-mindImgDataPtr[mindIndex] / meanValue); + mindImgDataPtr[mindIndex] = descValue; + maxDesc = std::max(maxDesc, descValue); mindIndex += voxelNumber; } mindIndex = voxelIndex; for (int t = 0; t < samplingNbr; t++) { - descValue = MINDImgDataPtr[mindIndex]; - MINDImgDataPtr[mindIndex] = descValue / max_desc; + descValue = mindImgDataPtr[mindIndex]; + mindImgDataPtr[mindIndex] = descValue / maxDesc; mindIndex += voxelNumber; } } // mask } // voxIndex // Mr Propre - nifti_image_free(diff_image); + nifti_image_free(diffImage); nifti_image_free(shiftedImage); nifti_image_free(meanImage); currentInputImage->data = nullptr; nifti_image_free(currentInputImage); } /* *************************************************************** */ -void GetMINDImageDescriptor(nifti_image* inputImgPtr, - nifti_image* MINDImgPtr, - int *maskPtr, - int descriptorOffset, - int currentTimepoint) { +void GetMindImageDescriptor(const nifti_image *inputImage, + nifti_image *mindImage, + const int *mask, + const int& descriptorOffset, + const int& currentTimepoint) { #ifndef NDEBUG - reg_print_fct_debug("GetMINDImageDescriptor()"); + reg_print_fct_debug("GetMindImageDescriptor()"); #endif - if (inputImgPtr->datatype != MINDImgPtr->datatype) { - reg_print_fct_error("reg_mind -- GetMINDImageDescriptor"); + if (inputImage->datatype != mindImage->datatype) { + reg_print_fct_error("reg_mind::GetMindImageDescriptor"); reg_print_msg_error("The input image and the MIND image must have the same datatype !"); reg_exit(); } - switch (inputImgPtr->datatype) { + switch (inputImage->datatype) { case NIFTI_TYPE_FLOAT32: - GetMINDImageDescriptor_core(inputImgPtr, MINDImgPtr, maskPtr, descriptorOffset, currentTimepoint); + GetMindImageDescriptorCore(inputImage, mindImage, mask, descriptorOffset, currentTimepoint); break; case NIFTI_TYPE_FLOAT64: - GetMINDImageDescriptor_core(inputImgPtr, MINDImgPtr, maskPtr, descriptorOffset, currentTimepoint); + GetMindImageDescriptorCore(inputImage, mindImage, mask, descriptorOffset, currentTimepoint); break; default: - reg_print_fct_error("GetMINDImageDescriptor"); + reg_print_fct_error("GetMindImageDescriptor"); reg_print_msg_error("Input image datatype not supported"); reg_exit(); break; @@ -194,11 +191,11 @@ void GetMINDImageDescriptor(nifti_image* inputImgPtr, } /* *************************************************************** */ template -void GetMINDSSCImageDescriptor_core(nifti_image* inputImage, - nifti_image* MINDSSCImage, - int *maskPtr, - int descriptorOffset, - int currentTimepoint) { +void GetMindSscImageDescriptorCore(const nifti_image *inputImage, + nifti_image *mindSscImage, + const int *mask, + const int& descriptorOffset, + const int& currentTimepoint) { #ifdef WIN32 long voxelIndex; const long voxelNumber = (long)NiftiImage::calcVoxelNumber(inputImage, 3); @@ -208,7 +205,7 @@ void GetMINDSSCImageDescriptor_core(nifti_image* inputImage, #endif // Create a pointer to the descriptor image - DataType* MINDSSCImgDataPtr = static_cast(MINDSSCImage->data); + DataType* mindSscImgDataPtr = static_cast(mindSscImage->data); // Allocate an image to store the current timepoint reference image nifti_image *currentInputImage = nifti_copy_nim_info(inputImage); @@ -216,18 +213,17 @@ void GetMINDSSCImageDescriptor_core(nifti_image* inputImage, currentInputImage->nt = currentInputImage->dim[4] = 1; currentInputImage->nvox = voxelNumber; DataType *inputImagePtr = static_cast(inputImage->data); - currentInputImage->data = static_cast(&inputImagePtr[currentTimepoint * voxelNumber]); + currentInputImage->data = &inputImagePtr[currentTimepoint * voxelNumber]; // Allocate an image to store the mean image - nifti_image *mean_img = nifti_dup(*currentInputImage, false); - DataType* meanImgDataPtr = static_cast(mean_img->data); + nifti_image *meanImg = nifti_dup(*currentInputImage, false); + DataType* meanImgDataPtr = static_cast(meanImg->data); // Allocate an image to store the warped image nifti_image *shiftedImage = nifti_dup(*currentInputImage, false); // Define the sigma for the convolution - float sigma = -0.5;// negative value denotes voxel width - //float sigma = -1.0;// negative value denotes voxel width + float sigma = -0.5; // negative value denotes voxel width //2D version int samplingNbr = (currentInputImage->nz > 1) ? 6 : 2; @@ -236,14 +232,14 @@ void GetMINDSSCImageDescriptor_core(nifti_image* inputImage, // Allocation of the difference image //std::vector vectNiftiImage; //for(int i=0;invox, sizeof(int)); + nifti_image *diffImage = nifti_dup(*currentInputImage, false); + int *maskDiffImage = (int*)calloc(diffImage->nvox, sizeof(int)); - nifti_image *diff_imageShifted = nifti_dup(*currentInputImage, false); + nifti_image *diffImageShifted = nifti_dup(*currentInputImage, false); - int RSampling3D_x[6] = { +descriptorOffset, +descriptorOffset, -descriptorOffset, +0, +descriptorOffset, +0 }; - int RSampling3D_y[6] = { +descriptorOffset, -descriptorOffset, +0, -descriptorOffset, +0, +descriptorOffset }; - int RSampling3D_z[6] = { +0, +0, +descriptorOffset, +descriptorOffset, +descriptorOffset, +descriptorOffset }; + int rSamplingX[6] = { +descriptorOffset, +descriptorOffset, -descriptorOffset, +0, +descriptorOffset, +0 }; + int rSamplingY[6] = { +descriptorOffset, -descriptorOffset, +0, -descriptorOffset, +0, +descriptorOffset }; + int rSamplingZ[6] = { +0, +0, +descriptorOffset, +descriptorOffset, +descriptorOffset, +descriptorOffset }; int tx[12] = { -descriptorOffset, +0, -descriptorOffset, +0, +0, +descriptorOffset, +0, +0, +0, -descriptorOffset, +0, +0 }; int ty[12] = { +0, -descriptorOffset, +0, +descriptorOffset, +0, +0, +0, +descriptorOffset, +0, +0, +0, -descriptorOffset }; @@ -251,94 +247,91 @@ void GetMINDSSCImageDescriptor_core(nifti_image* inputImage, int compteurId = 0; for (int i = 0; i < samplingNbr; i++) { - ShiftImage(currentInputImage, shiftedImage, maskPtr, - RSampling3D_x[i], RSampling3D_y[i], RSampling3D_z[i]); - reg_tools_subtractImageFromImage(currentInputImage, shiftedImage, diff_image); - reg_tools_multiplyImageToImage(diff_image, diff_image, diff_image); - reg_tools_kernelConvolution(diff_image, &sigma, GAUSSIAN_KERNEL, maskPtr); + ShiftImage(currentInputImage, shiftedImage, mask, rSamplingX[i], rSamplingY[i], rSamplingZ[i]); + reg_tools_subtractImageFromImage(currentInputImage, shiftedImage, diffImage); + reg_tools_multiplyImageToImage(diffImage, diffImage, diffImage); + reg_tools_kernelConvolution(diffImage, &sigma, GAUSSIAN_KERNEL, mask); for (int j = 0; j < 2; j++) { - - ShiftImage(diff_image, diff_imageShifted, mask_diff_image, + ShiftImage(diffImage, diffImageShifted, maskDiffImage, tx[compteurId], ty[compteurId], tz[compteurId]); - reg_tools_addImageToImage(mean_img, diff_imageShifted, mean_img); + reg_tools_addImageToImage(meanImg, diffImageShifted, meanImg); // Store the current descriptor - const size_t index = compteurId * diff_imageShifted->nvox; - memcpy(&MINDSSCImgDataPtr[index], diff_imageShifted->data, - diff_imageShifted->nbyper * diff_imageShifted->nvox); + const size_t index = compteurId * diffImageShifted->nvox; + memcpy(&mindSscImgDataPtr[index], diffImageShifted->data, + diffImageShifted->nbyper * diffImageShifted->nvox); compteurId++; } } // Compute the mean over the number of sample - reg_tools_divideValueToImage(mean_img, mean_img, lengthDescriptor); + reg_tools_divideValueToImage(meanImg, meanImg, lengthDescriptor); - // Compute the MINDSSC descriptor + // Compute the MIND-SSC descriptor int mindIndex; - DataType meanValue, max_desc, descValue; + DataType meanValue, maxDesc, descValue; #ifdef _OPENMP #pragma omp parallel for default(none) \ - shared(voxelNumber, lengthDescriptor, samplingNbr, maskPtr, meanImgDataPtr, \ - MINDSSCImgDataPtr) \ - private(meanValue, max_desc, descValue, mindIndex) + shared(voxelNumber, lengthDescriptor, samplingNbr, mask, meanImgDataPtr, mindSscImgDataPtr) \ + private(meanValue, maxDesc, descValue, mindIndex) #endif for (voxelIndex = 0; voxelIndex < voxelNumber; voxelIndex++) { - if (maskPtr[voxelIndex] > -1) { + if (mask[voxelIndex] > -1) { // Get the mean value for the current voxel meanValue = meanImgDataPtr[voxelIndex]; if (meanValue == 0) { meanValue = std::numeric_limits::epsilon(); } - max_desc = 0; + maxDesc = 0; mindIndex = voxelIndex; for (int t = 0; t < lengthDescriptor; t++) { - descValue = (DataType)exp(-MINDSSCImgDataPtr[mindIndex] / meanValue); - MINDSSCImgDataPtr[mindIndex] = descValue; - max_desc = std::max(max_desc, descValue); + descValue = (DataType)exp(-mindSscImgDataPtr[mindIndex] / meanValue); + mindSscImgDataPtr[mindIndex] = descValue; + maxDesc = std::max(maxDesc, descValue); mindIndex += voxelNumber; } mindIndex = voxelIndex; for (int t = 0; t < lengthDescriptor; t++) { - descValue = MINDSSCImgDataPtr[mindIndex]; - MINDSSCImgDataPtr[mindIndex] = descValue / max_desc; + descValue = mindSscImgDataPtr[mindIndex]; + mindSscImgDataPtr[mindIndex] = descValue / maxDesc; mindIndex += voxelNumber; } } // mask } // voxIndex // Mr Propre - nifti_image_free(diff_imageShifted); - free(mask_diff_image); - nifti_image_free(diff_image); + nifti_image_free(diffImageShifted); + free(maskDiffImage); + nifti_image_free(diffImage); nifti_image_free(shiftedImage); - nifti_image_free(mean_img); + nifti_image_free(meanImg); currentInputImage->data = nullptr; nifti_image_free(currentInputImage); } /* *************************************************************** */ -void GetMINDSSCImageDescriptor(nifti_image* inputImgPtr, - nifti_image* MINDSSCImgPtr, - int *maskPtr, - int descriptorOffset, - int currentTimepoint) { +void GetMindSscImageDescriptor(const nifti_image *inputImage, + nifti_image *mindSscImage, + const int *mask, + const int& descriptorOffset, + const int& currentTimepoint) { #ifndef NDEBUG - reg_print_fct_debug("GetMINDSSCImageDescriptor()"); + reg_print_fct_debug("GetMindSscImageDescriptor()"); #endif - if (inputImgPtr->datatype != MINDSSCImgPtr->datatype) { - reg_print_fct_error("reg_mindssc -- GetMINDSSCImageDescriptor"); - reg_print_msg_error("The input image and the MINDSSC image must have the same datatype !"); + if (inputImage->datatype != mindSscImage->datatype) { + reg_print_fct_error("reg_mindssc::GetMindSscImageDescriptor"); + reg_print_msg_error("The input image and the MINDSSC image must have the same datatype!"); reg_exit(); } - switch (inputImgPtr->datatype) { + switch (inputImage->datatype) { case NIFTI_TYPE_FLOAT32: - GetMINDSSCImageDescriptor_core(inputImgPtr, MINDSSCImgPtr, maskPtr, descriptorOffset, currentTimepoint); + GetMindSscImageDescriptorCore(inputImage, mindSscImage, mask, descriptorOffset, currentTimepoint); break; case NIFTI_TYPE_FLOAT64: - GetMINDSSCImageDescriptor_core(inputImgPtr, MINDSSCImgPtr, maskPtr, descriptorOffset, currentTimepoint); + GetMindSscImageDescriptorCore(inputImage, mindSscImage, mask, descriptorOffset, currentTimepoint); break; default: - reg_print_fct_error("GetMINDSSCImageDescriptor"); + reg_print_fct_error("GetMindSscImageDescriptor"); reg_print_msg_error("Input image datatype not supported"); reg_exit(); break; @@ -350,7 +343,7 @@ reg_mind::reg_mind(): reg_ssd() { this->floatingImageDescriptor = nullptr; this->warpedFloatingImageDescriptor = nullptr; this->warpedReferenceImageDescriptor = nullptr; - this->mind_type = MIND_TYPE; + this->mindType = MIND_TYPE; this->descriptorOffset = 1; #ifndef NDEBUG reg_print_msg_debug("reg_mind constructor called"); @@ -408,23 +401,22 @@ void reg_mind::InitialiseMeasure(nifti_image *refImg, warpedGradBw, voxelBasedGradBw); - this->descriptor_number = 0; - if (this->mind_type == MIND_TYPE) { - descriptor_number = this->referenceImage->nz > 1 ? 6 : 4; - } else if (this->mind_type == MINDSSC_TYPE) { - descriptor_number = this->referenceImage->nz > 1 ? 12 : 4; - + this->descriptorNumber = 0; + if (this->mindType == MIND_TYPE) { + this->descriptorNumber = this->referenceImage->nz > 1 ? 6 : 4; + } else if (this->mindType == MINDSSC_TYPE) { + this->descriptorNumber = this->referenceImage->nz > 1 ? 12 : 4; } // Initialise the reference descriptor this->referenceImageDescriptor = nifti_copy_nim_info(this->referenceImage); this->referenceImageDescriptor->dim[0] = this->referenceImageDescriptor->ndim = 4; - this->referenceImageDescriptor->dim[4] = this->referenceImageDescriptor->nt = this->descriptor_number; + this->referenceImageDescriptor->dim[4] = this->referenceImageDescriptor->nt = this->descriptorNumber; this->referenceImageDescriptor->nvox = NiftiImage::calcVoxelNumber(this->referenceImageDescriptor, this->referenceImageDescriptor->ndim); this->referenceImageDescriptor->data = malloc(this->referenceImageDescriptor->nvox * this->referenceImageDescriptor->nbyper); // Initialise the warped floating descriptor this->warpedFloatingImageDescriptor = nifti_copy_nim_info(this->referenceImage); this->warpedFloatingImageDescriptor->dim[0] = this->warpedFloatingImageDescriptor->ndim = 4; - this->warpedFloatingImageDescriptor->dim[4] = this->warpedFloatingImageDescriptor->nt = this->descriptor_number; + this->warpedFloatingImageDescriptor->dim[4] = this->warpedFloatingImageDescriptor->nt = this->descriptorNumber; this->warpedFloatingImageDescriptor->nvox = NiftiImage::calcVoxelNumber(this->warpedFloatingImageDescriptor, this->warpedFloatingImageDescriptor->ndim); this->warpedFloatingImageDescriptor->data = malloc(this->warpedFloatingImageDescriptor->nvox * @@ -438,7 +430,7 @@ void reg_mind::InitialiseMeasure(nifti_image *refImg, // Initialise the floating descriptor this->floatingImageDescriptor = nifti_copy_nim_info(this->floatingImage); this->floatingImageDescriptor->dim[0] = this->floatingImageDescriptor->ndim = 4; - this->floatingImageDescriptor->dim[4] = this->floatingImageDescriptor->nt = this->descriptor_number; + this->floatingImageDescriptor->dim[4] = this->floatingImageDescriptor->nt = this->descriptorNumber; this->floatingImageDescriptor->nvox = NiftiImage::calcVoxelNumber(this->floatingImageDescriptor, this->floatingImageDescriptor->ndim); this->floatingImageDescriptor->data = malloc(this->floatingImageDescriptor->nvox * @@ -446,7 +438,7 @@ void reg_mind::InitialiseMeasure(nifti_image *refImg, // Initialise the warped floating descriptor this->warpedReferenceImageDescriptor = nifti_copy_nim_info(this->floatingImage); this->warpedReferenceImageDescriptor->dim[0] = this->warpedReferenceImageDescriptor->ndim = 4; - this->warpedReferenceImageDescriptor->dim[4] = this->warpedReferenceImageDescriptor->nt = this->descriptor_number; + this->warpedReferenceImageDescriptor->dim[4] = this->warpedReferenceImageDescriptor->nt = this->descriptorNumber; this->warpedReferenceImageDescriptor->nvox = NiftiImage::calcVoxelNumber(this->warpedReferenceImageDescriptor, this->warpedReferenceImageDescriptor->ndim); this->warpedReferenceImageDescriptor->data = malloc(this->warpedReferenceImageDescriptor->nvox * @@ -459,7 +451,7 @@ void reg_mind::InitialiseMeasure(nifti_image *refImg, #ifndef NDEBUG char text[255]; - reg_print_msg_debug("reg_mind::InitialiseMeasure()."); + reg_print_msg_debug("reg_mind::InitialiseMeasure()"); sprintf(text, "Active time point:"); for (int i = 0; i < this->referenceImageDescriptor->nt; ++i) if (this->timePointWeightDescriptor[i] > 0) @@ -468,127 +460,82 @@ void reg_mind::InitialiseMeasure(nifti_image *refImg, #endif } /* *************************************************************** */ -double reg_mind::GetSimilarityMeasureValue() { - double MINDValue = 0.; - for (int t = 0; t < this->referenceImage->nt; ++t) { - if (this->timePointWeight[t] > 0) { - size_t voxelNumber = NiftiImage::calcVoxelNumber(referenceImage, 3); - int *combinedMask = (int*)malloc(voxelNumber * sizeof(int)); - memcpy(combinedMask, this->referenceMask, voxelNumber * sizeof(int)); - reg_tools_removeNanFromMask(this->referenceImage, combinedMask); - reg_tools_removeNanFromMask(this->warpedImage, combinedMask); - - if (this->mind_type == MIND_TYPE) { - GetMINDImageDescriptor(this->referenceImage, - this->referenceImageDescriptor, - combinedMask, - this->descriptorOffset, - t); - GetMINDImageDescriptor(this->warpedImage, - this->warpedFloatingImageDescriptor, - combinedMask, - this->descriptorOffset, - t); - } else if (this->mind_type == MINDSSC_TYPE) { - GetMINDSSCImageDescriptor(this->referenceImage, - this->referenceImageDescriptor, - combinedMask, - this->descriptorOffset, - t); - GetMINDSSCImageDescriptor(this->warpedImage, - this->warpedFloatingImageDescriptor, - combinedMask, - this->descriptorOffset, - t); - } - - switch (this->referenceImageDescriptor->datatype) { - case NIFTI_TYPE_FLOAT32: - MINDValue += reg_getSSDValue(this->referenceImageDescriptor, - this->warpedFloatingImageDescriptor, - this->timePointWeightDescriptor, - nullptr, // TODO this->forwardJacDetImagePointer, - combinedMask, - this->currentValue, - nullptr); - break; - case NIFTI_TYPE_FLOAT64: - MINDValue += reg_getSSDValue(this->referenceImageDescriptor, - this->warpedFloatingImageDescriptor, - this->timePointWeightDescriptor, - nullptr, // TODO this->forwardJacDetImagePointer, - combinedMask, - this->currentValue, - nullptr); - break; - default: - reg_print_fct_error("reg_mind::GetSimilarityMeasureValue"); - reg_print_msg_error("Warped pixel type unsupported"); - reg_exit(); - } - free(combinedMask); - - // Backward computation - if (this->isSymmetric) { - voxelNumber = NiftiImage::calcVoxelNumber(floatingImage, 3); - combinedMask = (int*)malloc(voxelNumber * sizeof(int)); - memcpy(combinedMask, this->floatingMask, voxelNumber * sizeof(int)); - reg_tools_removeNanFromMask(this->floatingImage, combinedMask); - reg_tools_removeNanFromMask(this->warpedImageBw, combinedMask); - - if (this->mind_type == MIND_TYPE) { - GetMINDImageDescriptor(this->floatingImage, - this->floatingImageDescriptor, - combinedMask, - this->descriptorOffset, - t); - GetMINDImageDescriptor(this->warpedImageBw, - this->warpedReferenceImageDescriptor, - combinedMask, - this->descriptorOffset, - t); - } else if (this->mind_type == MINDSSC_TYPE) { - GetMINDSSCImageDescriptor(this->floatingImage, - this->floatingImageDescriptor, - combinedMask, - this->descriptorOffset, - t); - GetMINDSSCImageDescriptor(this->warpedImageBw, - this->warpedReferenceImageDescriptor, - combinedMask, - this->descriptorOffset, - t); - } +double GetSimilarityMeasureValue(nifti_image *referenceImage, + nifti_image *referenceImageDescriptor, + const int *referenceMask, + nifti_image *warpedImage, + nifti_image *warpedFloatingImageDescriptor, + const double *timePointWeight, + double *timePointWeightDescriptor, + nifti_image *jacobianDetImage, + float *currentValue, + int descriptorOffset, + const int& referenceTimePoint, + const int& mindType) { + if (referenceImageDescriptor->datatype != NIFTI_TYPE_FLOAT32 && + referenceImageDescriptor->datatype != NIFTI_TYPE_FLOAT64) { + reg_print_fct_error("reg_mind::GetSimilarityMeasureValue"); + reg_print_msg_error("The reference image descriptor is expected to be of floating precision type"); + reg_exit(); + } - switch (this->floatingImageDescriptor->datatype) { - case NIFTI_TYPE_FLOAT32: - MINDValue += reg_getSSDValue(this->floatingImageDescriptor, - this->warpedReferenceImageDescriptor, - this->timePointWeightDescriptor, - nullptr, // TODO this->backwardJacDetImagePointer, - combinedMask, - this->currentValue, + double mind = 0; + const size_t voxelNumber = NiftiImage::calcVoxelNumber(referenceImage, 3); + unique_ptr combinedMask(new int[voxelNumber]); + auto GetMindImgDesc = mindType == MIND_TYPE ? GetMindImageDescriptor : GetMindSscImageDescriptor; + + for (int currentTimepoint = 0; currentTimepoint < referenceTimePoint; ++currentTimepoint) { + if (timePointWeight[currentTimepoint] > 0) { + memcpy(combinedMask.get(), referenceMask, voxelNumber * sizeof(int)); + reg_tools_removeNanFromMask(referenceImage, combinedMask.get()); + reg_tools_removeNanFromMask(warpedImage, combinedMask.get()); + + GetMindImgDesc(referenceImage, referenceImageDescriptor, combinedMask.get(), descriptorOffset, currentTimepoint); + GetMindImgDesc(warpedImage, warpedFloatingImageDescriptor, combinedMask.get(), descriptorOffset, currentTimepoint); + + std::visit([&](auto&& refImgDataType) { + using RefImgDataType = std::decay_t; + mind += reg_getSsdValue(referenceImageDescriptor, + warpedFloatingImageDescriptor, + timePointWeightDescriptor, + jacobianDetImage, + combinedMask.get(), + currentValue, nullptr); - break; - case NIFTI_TYPE_FLOAT64: - MINDValue += reg_getSSDValue(this->floatingImageDescriptor, - this->warpedReferenceImageDescriptor, - this->timePointWeightDescriptor, - nullptr, // TODO this->backwardJacDetImagePointer, - combinedMask, - this->currentValue, - nullptr); - break; - default: - reg_print_fct_error("reg_mind::GetSimilarityMeasureValue"); - reg_print_msg_error("Warped pixel type unsupported"); - reg_exit(); - } - free(combinedMask); - } + }, NiftiImage::getFloatingDataType(referenceImageDescriptor)); } } - return MINDValue; // (double) this->referenceImageDescriptor->nt; + return mind; +} +/* *************************************************************** */ +double reg_mind::GetSimilarityMeasureValueFw() { + return ::GetSimilarityMeasureValue(this->referenceImage, + this->referenceImageDescriptor, + this->referenceMask, + this->warpedImage, + this->warpedFloatingImageDescriptor, + this->timePointWeight, + this->timePointWeightDescriptor, + nullptr, // TODO this->forwardJacDetImagePointer, + this->currentValue, + this->descriptorOffset, + this->referenceTimePoint, + this->mindType); +} +/* *************************************************************** */ +double reg_mind::GetSimilarityMeasureValueBw() { + return ::GetSimilarityMeasureValue(this->floatingImage, + this->floatingImageDescriptor, + this->floatingMask, + this->warpedImageBw, + this->warpedReferenceImageDescriptor, + this->timePointWeight, + this->timePointWeightDescriptor, + nullptr, // TODO this->backwardJacDetImagePointer, + this->currentValue, + this->descriptorOffset, + this->referenceTimePoint, + this->mindType); } /* *************************************************************** */ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { @@ -604,28 +551,28 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { reg_tools_removeNanFromMask(this->referenceImage, combinedMask); reg_tools_removeNanFromMask(this->warpedImage, combinedMask); - if (this->mind_type == MIND_TYPE) { + if (this->mindType == MIND_TYPE) { // Compute the reference image descriptors - GetMINDImageDescriptor(this->referenceImage, + GetMindImageDescriptor(this->referenceImage, this->referenceImageDescriptor, combinedMask, this->descriptorOffset, currentTimepoint); // Compute the warped floating image descriptors - GetMINDImageDescriptor(this->warpedImage, + GetMindImageDescriptor(this->warpedImage, this->warpedFloatingImageDescriptor, combinedMask, this->descriptorOffset, currentTimepoint); - } else if (this->mind_type == MINDSSC_TYPE) { + } else if (this->mindType == MINDSSC_TYPE) { // Compute the reference image descriptors - GetMINDSSCImageDescriptor(this->referenceImage, + GetMindSscImageDescriptor(this->referenceImage, this->referenceImageDescriptor, combinedMask, this->descriptorOffset, currentTimepoint); // Compute the warped floating image descriptors - GetMINDSSCImageDescriptor(this->warpedImage, + GetMindSscImageDescriptor(this->warpedImage, this->warpedFloatingImageDescriptor, combinedMask, this->descriptorOffset, @@ -633,7 +580,7 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { } - for (int desc_index = 0; desc_index < this->descriptor_number; ++desc_index) { + for (int desc_index = 0; desc_index < this->descriptorNumber; ++desc_index) { // Compute the warped image descriptors gradient reg_getImageGradient_symDiff(this->warpedFloatingImageDescriptor, this->warpedGradient, @@ -644,7 +591,7 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { // Compute the gradient of the ssd for the forward transformation switch (referenceImageDescriptor->datatype) { case NIFTI_TYPE_FLOAT32: - reg_getVoxelBasedSSDGradient(this->referenceImageDescriptor, + reg_getVoxelBasedSsdGradient(this->referenceImageDescriptor, this->warpedFloatingImageDescriptor, this->warpedGradient, this->voxelBasedGradient, @@ -655,7 +602,7 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { nullptr); break; case NIFTI_TYPE_FLOAT64: - reg_getVoxelBasedSSDGradient(this->referenceImageDescriptor, + reg_getVoxelBasedSsdGradient(this->referenceImageDescriptor, this->warpedFloatingImageDescriptor, this->warpedGradient, this->voxelBasedGradient, @@ -681,31 +628,31 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { reg_tools_removeNanFromMask(this->floatingImage, combinedMask); reg_tools_removeNanFromMask(this->warpedImageBw, combinedMask); - if (this->mind_type == MIND_TYPE) { - GetMINDImageDescriptor(this->floatingImage, + if (this->mindType == MIND_TYPE) { + GetMindImageDescriptor(this->floatingImage, this->floatingImageDescriptor, combinedMask, this->descriptorOffset, currentTimepoint); - GetMINDImageDescriptor(this->warpedImageBw, + GetMindImageDescriptor(this->warpedImageBw, this->warpedReferenceImageDescriptor, combinedMask, this->descriptorOffset, currentTimepoint); - } else if (this->mind_type == MINDSSC_TYPE) { - GetMINDSSCImageDescriptor(this->floatingImage, + } else if (this->mindType == MINDSSC_TYPE) { + GetMindSscImageDescriptor(this->floatingImage, this->floatingImageDescriptor, combinedMask, this->descriptorOffset, currentTimepoint); - GetMINDSSCImageDescriptor(this->warpedImageBw, + GetMindSscImageDescriptor(this->warpedImageBw, this->warpedReferenceImageDescriptor, combinedMask, this->descriptorOffset, currentTimepoint); } - for (int desc_index = 0; desc_index < this->descriptor_number; ++desc_index) { + for (int desc_index = 0; desc_index < this->descriptorNumber; ++desc_index) { reg_getImageGradient_symDiff(this->warpedReferenceImageDescriptor, this->warpedGradientBw, combinedMask, @@ -715,7 +662,7 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { // Compute the gradient of the nmi for the backward transformation switch (floatingImage->datatype) { case NIFTI_TYPE_FLOAT32: - reg_getVoxelBasedSSDGradient(this->floatingImageDescriptor, + reg_getVoxelBasedSsdGradient(this->floatingImageDescriptor, this->warpedReferenceImageDescriptor, this->warpedGradientBw, this->voxelBasedGradientBw, @@ -726,7 +673,7 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { nullptr); break; case NIFTI_TYPE_FLOAT64: - reg_getVoxelBasedSSDGradient(this->floatingImageDescriptor, + reg_getVoxelBasedSsdGradient(this->floatingImageDescriptor, this->warpedReferenceImageDescriptor, this->warpedGradientBw, this->voxelBasedGradientBw, @@ -747,7 +694,7 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { } /* *************************************************************** */ reg_mindssc::reg_mindssc(): reg_mind() { - this->mind_type = MINDSSC_TYPE; + this->mindType = MINDSSC_TYPE; #ifndef NDEBUG reg_print_msg_debug("reg_mindssc constructor called"); #endif diff --git a/reg-lib/cpu/_reg_mind.h b/reg-lib/cpu/_reg_mind.h index cf09a4a8..9eb88336 100644 --- a/reg-lib/cpu/_reg_mind.h +++ b/reg-lib/cpu/_reg_mind.h @@ -30,18 +30,20 @@ class reg_mind: public reg_ssd { /// @brief Initialise the reg_mind object virtual void InitialiseMeasure(nifti_image *refImg, - nifti_image *floImg, - int *refMask, - nifti_image *warpedImg, - nifti_image *warpedGrad, - nifti_image *voxelBasedGrad, - nifti_image *localWeightSim = nullptr, - int *floMask = nullptr, - nifti_image *warpedImgBw = nullptr, - nifti_image *warpedGradBw = nullptr, - nifti_image *voxelBasedGradBw = nullptr) override; - /// @brief Returns the mind based measure of similarity value - virtual double GetSimilarityMeasureValue() override; + nifti_image *floImg, + int *refMask, + nifti_image *warpedImg, + nifti_image *warpedGrad, + nifti_image *voxelBasedGrad, + nifti_image *localWeightSim = nullptr, + int *floMask = nullptr, + nifti_image *warpedImgBw = nullptr, + nifti_image *warpedGradBw = nullptr, + nifti_image *voxelBasedGradBw = nullptr) override; + /// @brief Returns the forward mind-based measure of similarity value + virtual double GetSimilarityMeasureValueFw() override; + /// @brief Returns the backward mind-based measure of similarity value + virtual double GetSimilarityMeasureValueBw() override; /// @brief Compute the voxel based gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override; virtual void SetDescriptorOffset(int); @@ -55,8 +57,8 @@ class reg_mind: public reg_ssd { double timePointWeightDescriptor[255] = {0}; int descriptorOffset; - int mind_type; - int descriptor_number; + int mindType; + int descriptorNumber; }; /* *************************************************************** */ /// @brief MIND-SSC measure of similarity class @@ -69,16 +71,16 @@ class reg_mindssc: public reg_mind { }; /* *************************************************************** */ extern "C++" -void GetMINDImageDescriptor(nifti_image *inputImgPtr, - nifti_image *MINDImgPtr, - int *mask, - int descriptorOffset, - int currentTimepoint); +void GetMindImageDescriptor(const nifti_image *inputImage, + nifti_image *mindImage, + const int *mask, + const int& descriptorOffset, + const int& currentTimepoint); /* *************************************************************** */ extern "C++" -void GetMINDSSCImageDescriptor(nifti_image *inputImgPtr, - nifti_image *MINDSSCImgPtr, - int *mask, - int descriptorOffset, - int currentTimepoint); +void GetMindSscImageDescriptor(const nifti_image *inputImage, + nifti_image *mindSscImage, + const int *mask, + const int& descriptorOffset, + const int& currentTimepoint); /* *************************************************************** */ diff --git a/reg-lib/cpu/_reg_nmi.cpp b/reg-lib/cpu/_reg_nmi.cpp index 23288d73..4036cf08 100755 --- a/reg-lib/cpu/_reg_nmi.cpp +++ b/reg-lib/cpu/_reg_nmi.cpp @@ -213,7 +213,7 @@ void reg_getNMIValue(const nifti_image *referenceImage, const unsigned short *floatingBinNumber, const unsigned short *totalBinNumber, double **jointHistogramLog, - double **jointhistogramPro, + double **jointHistogramPro, double **entropyValues, const int *referenceMask) { // Create pointers to the image data arrays @@ -230,7 +230,7 @@ void reg_getNMIValue(const nifti_image *referenceImage, reg_print_msg_debug(text); #endif // Define some pointers to the current histograms - double *jointHistoProPtr = jointhistogramPro[t]; + double *jointHistoProPtr = jointHistogramPro[t]; double *jointHistoLogPtr = jointHistogramLog[t]; // Empty the joint histogram memset(jointHistoProPtr, 0, totalBinNumber[t] * sizeof(double)); @@ -355,71 +355,65 @@ void reg_getNMIValue(const nifti_image *referenceImage, } // iterate over all time point in the reference image } /* *************************************************************** */ -double reg_nmi::GetSimilarityMeasureValue() { - // Check that all the specified image are of the same datatype - if (this->referenceImage->datatype != NIFTI_TYPE_FLOAT32 && this->referenceImage->datatype != NIFTI_TYPE_FLOAT64) { - reg_print_fct_error("reg_nmi::GetSimilarityMeasureValue()"); - reg_print_msg_error("Input images are expected to be of floating precision type"); - reg_exit(); - } - if (this->warpedImage->datatype != this->referenceImage->datatype) { - reg_print_fct_error("reg_nmi::GetSimilarityMeasureValue()"); - reg_print_msg_error("Both input images are expected to have the same type"); - reg_exit(); - } +double GetSimilarityMeasureValue(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const double *timePointWeight, + const unsigned short *referenceBinNumber, + const unsigned short *floatingBinNumber, + const unsigned short *totalBinNumber, + double **jointHistogramLog, + double **jointHistogramPro, + double **entropyValues, + const int *referenceMask, + const int& referenceTimePoint) { std::visit([&](auto&& refImgDataType) { using RefImgDataType = std::decay_t; - reg_getNMIValue(this->referenceImage, - this->warpedImage, - this->timePointWeight, - this->referenceBinNumber, - this->floatingBinNumber, - this->totalBinNumber, - this->jointHistogramLog, - this->jointHistogramPro, - this->entropyValues, - this->referenceMask); - }, NiftiImage::getFloatingDataType(this->referenceImage)); - - if (this->isSymmetric) { - // Check that all the specified image are of the same datatype - if (this->floatingImage->datatype != NIFTI_TYPE_FLOAT32 && this->floatingImage->datatype != NIFTI_TYPE_FLOAT64) { - reg_print_fct_error("reg_nmi::GetSimilarityMeasureValue()"); - reg_print_msg_error("Input images are expected to be of floating precision type"); - reg_exit(); - } - if (this->floatingImage->datatype != this->warpedImageBw->datatype) { - reg_print_fct_error("reg_nmi::GetSimilarityMeasureValue()"); - reg_print_msg_error("Both input images are expected to have the same type"); - reg_exit(); - } - std::visit([&](auto&& floImgDataType) { - using FloImgDataType = std::decay_t; - reg_getNMIValue(this->floatingImage, - this->warpedImageBw, - this->timePointWeight, - this->floatingBinNumber, - this->referenceBinNumber, - this->totalBinNumber, - this->jointHistogramLogBw, - this->jointHistogramProBw, - this->entropyValuesBw, - this->floatingMask); - }, NiftiImage::getFloatingDataType(this->floatingImage)); - } + reg_getNMIValue(referenceImage, + warpedImage, + timePointWeight, + referenceBinNumber, + floatingBinNumber, + totalBinNumber, + jointHistogramLog, + jointHistogramPro, + entropyValues, + referenceMask); + }, NiftiImage::getFloatingDataType(referenceImage)); - double nmiFw = 0, nmiBw = 0; - for (int t = 0; t < this->referenceTimePoint; ++t) { - if (this->timePointWeight[t] > 0) { - nmiFw += timePointWeight[t] * (this->entropyValues[t][0] + this->entropyValues[t][1]) / this->entropyValues[t][2]; - if (this->isSymmetric) - nmiBw += timePointWeight[t] * (this->entropyValuesBw[t][0] + this->entropyValuesBw[t][1]) / this->entropyValuesBw[t][2]; - } + double nmi = 0; + for (int t = 0; t < referenceTimePoint; ++t) { + if (timePointWeight[t] > 0) + nmi += timePointWeight[t] * (entropyValues[t][0] + entropyValues[t][1]) / entropyValues[t][2]; } -#ifndef NDEBUG - reg_print_msg_debug("reg_nmi::GetSimilarityMeasureValue called"); -#endif - return nmiFw + nmiBw; + return nmi; +} +/* *************************************************************** */ +double reg_nmi::GetSimilarityMeasureValueFw() { + return ::GetSimilarityMeasureValue(this->referenceImage, + this->warpedImage, + this->timePointWeight, + this->referenceBinNumber, + this->floatingBinNumber, + this->totalBinNumber, + this->jointHistogramLog, + this->jointHistogramPro, + this->entropyValues, + this->referenceMask, + this->referenceTimePoint); +} +/* *************************************************************** */ +double reg_nmi::GetSimilarityMeasureValueBw() { + return ::GetSimilarityMeasureValue(this->floatingImage, + this->warpedImageBw, + this->timePointWeight, + this->floatingBinNumber, + this->referenceBinNumber, + this->totalBinNumber, + this->jointHistogramLogBw, + this->jointHistogramProBw, + this->entropyValuesBw, + this->floatingMask, + this->referenceTimePoint); } /* *************************************************************** */ template diff --git a/reg-lib/cpu/_reg_nmi.h b/reg-lib/cpu/_reg_nmi.h index 78cd06ad..3f66e70e 100755 --- a/reg-lib/cpu/_reg_nmi.h +++ b/reg-lib/cpu/_reg_nmi.h @@ -38,8 +38,10 @@ class reg_nmi: public reg_measure { nifti_image *warpedImgBw = nullptr, nifti_image *warpedGradBw = nullptr, nifti_image *voxelBasedGradBw = nullptr) override; - /// @brief Returns the nmi value - virtual double GetSimilarityMeasureValue() override; + /// @brief Returns the nmi value forwards + virtual double GetSimilarityMeasureValueFw() override; + /// @brief Returns the nmi value backwards + virtual double GetSimilarityMeasureValueBw() override; /// @brief Compute the voxel based nmi gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override; @@ -84,7 +86,7 @@ void reg_getNMIValue(const nifti_image *referenceImage, const unsigned short *floatingBinNumber, const unsigned short *totalBinNumber, double **jointHistogramLog, - double **jointhistogramPro, + double **jointHistogramPro, double **entropyValues, const int *referenceMask); /* *************************************************************** */ @@ -213,8 +215,10 @@ class reg_multichannel_nmi: public reg_measure { /// @brief reg_multichannel_nmi class destructor virtual ~reg_multichannel_nmi() {} - /// @brief Returns the nmi value - virtual double GetSimilarityMeasureValue() override { return 0; } + /// @brief Returns the nmi value forwards + virtual double GetSimilarityMeasureValueFw() override { return 0; } + /// @brief Returns the nmi value backwards + virtual double GetSimilarityMeasureValueBw() override { return 0; } /// @brief Compute the voxel based nmi gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override { diff --git a/reg-lib/cpu/_reg_ssd.cpp b/reg-lib/cpu/_reg_ssd.cpp index ac3a3a4b..9287934a 100755 --- a/reg-lib/cpu/_reg_ssd.cpp +++ b/reg-lib/cpu/_reg_ssd.cpp @@ -80,7 +80,7 @@ void reg_ssd::InitialiseMeasure(nifti_image *refImg, #endif #ifndef NDEBUG char text[255]; - reg_print_msg_debug("reg_ssd::InitialiseMeasure()."); + reg_print_msg_debug("reg_ssd::InitialiseMeasure()"); for (int i = 0; i < this->referenceImage->nt; ++i) { sprintf(text, "Weight for timepoint %i: %f", i, this->timePointWeight[i]); reg_print_msg_debug(text); @@ -98,13 +98,13 @@ void reg_ssd::SetNormaliseTimepoint(int timepoint, bool normalise) { } /* *************************************************************** */ template -double reg_getSSDValue(nifti_image *referenceImage, - nifti_image *warpedImage, - double *timePointWeight, - nifti_image *jacobianDetImage, - int *mask, +double reg_getSsdValue(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const double *timePointWeight, + const nifti_image *jacobianDetImage, + const int *mask, float *currentValue, - nifti_image *localWeightSimImage) { + const nifti_image *localWeightSim) { #ifdef _WIN32 long voxel; const long voxelNumber = (long)NiftiImage::calcVoxelNumber(referenceImage, 3); @@ -113,34 +113,34 @@ double reg_getSSDValue(nifti_image *referenceImage, const size_t voxelNumber = NiftiImage::calcVoxelNumber(referenceImage, 3); #endif // Create pointers to the reference and warped image data - DataType *referencePtr = static_cast(referenceImage->data); - DataType *warpedPtr = static_cast(warpedImage->data); + const DataType *referencePtr = static_cast(referenceImage->data); + const DataType *warpedPtr = static_cast(warpedImage->data); // Create a pointer to the Jacobian determinant image if defined - DataType *jacDetPtr = nullptr; + const DataType *jacDetPtr = nullptr; if (jacobianDetImage != nullptr) jacDetPtr = static_cast(jacobianDetImage->data); // Create a pointer to the local weight image if defined - DataType *localWeightPtr = nullptr; - if (localWeightSimImage != nullptr) - localWeightPtr = static_cast(localWeightSimImage->data); + const DataType *localWeightPtr = nullptr; + if (localWeightSim != nullptr) + localWeightPtr = static_cast(localWeightSim->data); - double SSD_global = 0; + double ssdGlobal = 0; double refValue, warValue, diff; // Loop over the different time points for (int time = 0; time < referenceImage->nt; ++time) { if (timePointWeight[time] > 0) { // Create pointers to the current time point of the reference and warped images - DataType *currentRefPtr = &referencePtr[time * voxelNumber]; - DataType *currentWarPtr = &warpedPtr[time * voxelNumber]; + const DataType *currentRefPtr = &referencePtr[time * voxelNumber]; + const DataType *currentWarPtr = &warpedPtr[time * voxelNumber]; - double SSD_local = 0., n = 0.; + double ssdLocal = 0, n = 0; #ifdef _OPENMP #pragma omp parallel for default(none) \ shared(referenceImage, warpedImage, currentRefPtr, currentWarPtr, mask, \ jacobianDetImage, jacDetPtr, voxelNumber, localWeightPtr) \ private(refValue, warValue, diff) \ - reduction(+:SSD_local) \ + reduction(+:ssdLocal) \ reduction(+:n) #endif for (voxel = 0; voxel < voxelNumber; ++voxel) { @@ -158,108 +158,76 @@ double reg_getSSDValue(nifti_image *referenceImage, #endif // Jacobian determinant modulation of the ssd if required if (jacDetPtr != nullptr) { - SSD_local += diff * jacDetPtr[voxel]; + ssdLocal += diff * jacDetPtr[voxel]; n += jacDetPtr[voxel]; } else if (localWeightPtr != nullptr) { - SSD_local += diff * localWeightPtr[voxel]; + ssdLocal += diff * localWeightPtr[voxel]; n += localWeightPtr[voxel]; } else { - SSD_local += diff; + ssdLocal += diff; n += 1.0; } } } } - SSD_local *= timePointWeight[time]; - currentValue[time] = static_cast(-SSD_local); - SSD_global -= SSD_local / n; + ssdLocal *= timePointWeight[time]; + currentValue[time] = static_cast(-ssdLocal); + ssdGlobal -= ssdLocal / n; } } - return SSD_global; + return ssdGlobal; } -template double reg_getSSDValue(nifti_image*, nifti_image*, double*, nifti_image*, int*, float*, nifti_image*); -template double reg_getSSDValue(nifti_image*, nifti_image*, double*, nifti_image*, int*, float*, nifti_image*); /* *************************************************************** */ -double reg_ssd::GetSimilarityMeasureValue() { - // Check that all the specified image are of the same datatype - if (this->warpedImage->datatype != this->referenceImage->datatype) { - reg_print_fct_error("reg_ssd::GetSimilarityMeasureValue"); - reg_print_msg_error("Both input images are expected to have the same type"); - reg_exit(); - } - double SSDValue = 0; - switch (this->referenceImage->datatype) { - case NIFTI_TYPE_FLOAT32: - SSDValue = reg_getSSDValue(this->referenceImage, - this->warpedImage, - this->timePointWeight, - nullptr, // TODO this->forwardJacDetImagePointer, - this->referenceMask, - this->currentValue, - this->localWeightSim); - break; - case NIFTI_TYPE_FLOAT64: - SSDValue = reg_getSSDValue(this->referenceImage, - this->warpedImage, - this->timePointWeight, - nullptr, // TODO this->forwardJacDetImagePointer, - this->referenceMask, - this->currentValue, - this->localWeightSim); - break; - default: - reg_print_fct_error("reg_ssd::GetSimilarityMeasureValue"); - reg_print_msg_error("Warped pixel type unsupported"); - reg_exit(); - } - - // Backward computation - if (this->isSymmetric) { - // Check that all the specified image are of the same datatype - if (this->warpedImageBw->datatype != this->floatingImage->datatype) { - reg_print_fct_error("reg_ssd::GetSimilarityMeasureValue"); - reg_print_msg_error("Both input images are expected to have the same type"); - reg_exit(); - } - switch (this->floatingImage->datatype) { - case NIFTI_TYPE_FLOAT32: - SSDValue += reg_getSSDValue(this->floatingImage, - this->warpedImageBw, - this->timePointWeight, - nullptr, // TODO this->backwardJacDetImagePointer, - this->floatingMask, - this->currentValue, - nullptr); - break; - case NIFTI_TYPE_FLOAT64: - SSDValue += reg_getSSDValue(this->floatingImage, - this->warpedImageBw, - this->timePointWeight, - nullptr, // TODO this->backwardJacDetImagePointer, - this->floatingMask, - this->currentValue, - nullptr); - break; - default: - reg_print_fct_error("reg_ssd::GetSimilarityMeasureValue"); - reg_print_msg_error("Warped pixel type unsupported"); - reg_exit(); - } - } - return SSDValue; +double GetSimilarityMeasureValue(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const double *timePointWeight, + const nifti_image *jacobianDetImage, + const int *mask, + float *currentValue, + const nifti_image *localWeightSim) { + return std::visit([&](auto&& refImgDataType) { + using RefImgDataType = std::decay_t; + return reg_getSsdValue(referenceImage, + warpedImage, + timePointWeight, + jacobianDetImage, + mask, + currentValue, + localWeightSim); + }, NiftiImage::getFloatingDataType(referenceImage)); +} +/* *************************************************************** */ +double reg_ssd::GetSimilarityMeasureValueFw() { + return ::GetSimilarityMeasureValue(this->referenceImage, + this->warpedImage, + this->timePointWeight, + nullptr, // TODO this->forwardJacDetImagePointer, + this->referenceMask, + this->currentValue, + this->localWeightSim); +} +/* *************************************************************** */ +double reg_ssd::GetSimilarityMeasureValueBw() { + return ::GetSimilarityMeasureValue(this->floatingImage, + this->warpedImageBw, + this->timePointWeight, + nullptr, // TODO this->backwardJacDetImagePointer, + this->floatingMask, + this->currentValue, + nullptr); } /* *************************************************************** */ template -void reg_getVoxelBasedSSDGradient(nifti_image *referenceImage, - nifti_image *warpedImage, - nifti_image *warpedGradient, +void reg_getVoxelBasedSsdGradient(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const nifti_image *warpedGradient, nifti_image *measureGradientImage, - nifti_image *jacobianDetImage, - int *mask, - int currentTimepoint, - double timepointWeight, - nifti_image *localWeightSimImage) { + const nifti_image *jacobianDetImage, + const int *mask, + const int& currentTimepoint, + const double& timepointWeight, + const nifti_image *localWeightSim) { if (currentTimepoint < 0 || currentTimepoint >= referenceImage->nt) { reg_print_fct_error("reg_getVoxelBasedSSDGradient"); reg_print_msg_error("The specified active timepoint is not defined in the ref/war images"); @@ -274,33 +242,33 @@ void reg_getVoxelBasedSSDGradient(nifti_image *referenceImage, const size_t voxelNumber = NiftiImage::calcVoxelNumber(referenceImage, 3); #endif // Pointers to the image data - DataType *refImagePtr = static_cast(referenceImage->data); - DataType *currentRefPtr = &refImagePtr[currentTimepoint * voxelNumber]; - DataType *warImagePtr = static_cast(warpedImage->data); - DataType *currentWarPtr = &warImagePtr[currentTimepoint * voxelNumber]; + const DataType *refImagePtr = static_cast(referenceImage->data); + const DataType *currentRefPtr = &refImagePtr[currentTimepoint * voxelNumber]; + const DataType *warImagePtr = static_cast(warpedImage->data); + const DataType *currentWarPtr = &warImagePtr[currentTimepoint * voxelNumber]; // Pointers to the spatial gradient of the warped image - DataType *spatialGradPtrX = static_cast(warpedGradient->data); - DataType *spatialGradPtrY = &spatialGradPtrX[voxelNumber]; - DataType *spatialGradPtrZ = nullptr; + const DataType *spatialGradPtrX = static_cast(warpedGradient->data); + const DataType *spatialGradPtrY = &spatialGradPtrX[voxelNumber]; + const DataType *spatialGradPtrZ = nullptr; if (referenceImage->nz > 1) spatialGradPtrZ = &spatialGradPtrY[voxelNumber]; // Pointers to the measure of similarity gradient - DataType *measureGradPtrX = static_cast(measureGradientImage->data); + DataType *measureGradPtrX = static_cast(measureGradientImage->data); DataType *measureGradPtrY = &measureGradPtrX[voxelNumber]; DataType *measureGradPtrZ = nullptr; if (referenceImage->nz > 1) measureGradPtrZ = &measureGradPtrY[voxelNumber]; // Create a pointer to the Jacobian determinant values if defined - DataType *jacDetPtr = nullptr; + const DataType *jacDetPtr = nullptr; if (jacobianDetImage != nullptr) - jacDetPtr = static_cast(jacobianDetImage->data); + jacDetPtr = static_cast(jacobianDetImage->data); // Create a pointer to the local weight image if defined - DataType *localWeightPtr = nullptr; - if (localWeightSimImage != nullptr) - localWeightPtr = static_cast(localWeightSimImage->data); + const DataType *localWeightPtr = nullptr; + if (localWeightSim != nullptr) + localWeightPtr = static_cast(localWeightSim->data); // find number of active voxels and correct weight double activeVoxel_num = 0; @@ -310,7 +278,7 @@ void reg_getVoxelBasedSSDGradient(nifti_image *referenceImage, activeVoxel_num += 1.0; } } - double adjusted_weight = timepointWeight / activeVoxel_num; + double adjustedWeight = timepointWeight / activeVoxel_num; double refValue, warValue, common; @@ -319,13 +287,13 @@ void reg_getVoxelBasedSSDGradient(nifti_image *referenceImage, shared(referenceImage, warpedImage, currentRefPtr, currentWarPtr, \ mask, jacDetPtr, spatialGradPtrX, spatialGradPtrY, spatialGradPtrZ, \ measureGradPtrX, measureGradPtrY, measureGradPtrZ, voxelNumber, \ - localWeightPtr, adjusted_weight) \ + localWeightPtr, adjustedWeight) \ private(refValue, warValue, common) #endif for (voxel = 0; voxel < voxelNumber; voxel++) { if (mask[voxel] > -1) { - refValue = (double)(currentRefPtr[voxel] * referenceImage->scl_slope + referenceImage->scl_inter); - warValue = (double)(currentWarPtr[voxel] * warpedImage->scl_slope + warpedImage->scl_inter); + refValue = currentRefPtr[voxel] * referenceImage->scl_slope + referenceImage->scl_inter; + warValue = currentWarPtr[voxel] * warpedImage->scl_slope + warpedImage->scl_inter; if (refValue == refValue && warValue == warValue) { #ifdef MRF_USE_SAD common = refValue > warValue ? -1.f : 1.f; @@ -338,25 +306,21 @@ void reg_getVoxelBasedSSDGradient(nifti_image *referenceImage, else if (localWeightPtr != nullptr) common *= localWeightPtr[voxel]; - common *= adjusted_weight; + common *= adjustedWeight; if (spatialGradPtrX[voxel] == spatialGradPtrX[voxel]) - measureGradPtrX[voxel] += (DataType)(common * spatialGradPtrX[voxel]); + measureGradPtrX[voxel] += static_cast(common * spatialGradPtrX[voxel]); if (spatialGradPtrY[voxel] == spatialGradPtrY[voxel]) - measureGradPtrY[voxel] += (DataType)(common * spatialGradPtrY[voxel]); + measureGradPtrY[voxel] += static_cast(common * spatialGradPtrY[voxel]); if (measureGradPtrZ != nullptr) { if (spatialGradPtrZ[voxel] == spatialGradPtrZ[voxel]) - measureGradPtrZ[voxel] += (DataType)(common * spatialGradPtrZ[voxel]); + measureGradPtrZ[voxel] += static_cast(common * spatialGradPtrZ[voxel]); } } } } } -template void reg_getVoxelBasedSSDGradient -(nifti_image*, nifti_image*, nifti_image*, nifti_image*, nifti_image*, int*, int, double, nifti_image*); -template void reg_getVoxelBasedSSDGradient -(nifti_image*, nifti_image*, nifti_image*, nifti_image*, nifti_image*, int*, int, double, nifti_image*); /* *************************************************************** */ void reg_ssd::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { // Check if the specified time point exists and is active @@ -376,7 +340,7 @@ void reg_ssd::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { // Compute the gradient of the ssd for the forward transformation switch (dtype) { case NIFTI_TYPE_FLOAT32: - reg_getVoxelBasedSSDGradient(this->referenceImage, + reg_getVoxelBasedSsdGradient(this->referenceImage, this->warpedImage, this->warpedGradient, this->voxelBasedGradient, @@ -387,7 +351,7 @@ void reg_ssd::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { this->localWeightSim); break; case NIFTI_TYPE_FLOAT64: - reg_getVoxelBasedSSDGradient(this->referenceImage, + reg_getVoxelBasedSsdGradient(this->referenceImage, this->warpedImage, this->warpedGradient, this->voxelBasedGradient, @@ -415,7 +379,7 @@ void reg_ssd::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { // Compute the gradient of the nmi for the backward transformation switch (dtype) { case NIFTI_TYPE_FLOAT32: - reg_getVoxelBasedSSDGradient(this->floatingImage, + reg_getVoxelBasedSsdGradient(this->floatingImage, this->warpedImageBw, this->warpedGradientBw, this->voxelBasedGradientBw, @@ -426,7 +390,7 @@ void reg_ssd::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { nullptr); break; case NIFTI_TYPE_FLOAT64: - reg_getVoxelBasedSSDGradient(this->floatingImage, + reg_getVoxelBasedSsdGradient(this->floatingImage, this->warpedImageBw, this->warpedGradientBw, this->voxelBasedGradientBw, diff --git a/reg-lib/cpu/_reg_ssd.h b/reg-lib/cpu/_reg_ssd.h index 5492f60c..43dbefe3 100755 --- a/reg-lib/cpu/_reg_ssd.h +++ b/reg-lib/cpu/_reg_ssd.h @@ -16,7 +16,6 @@ #include "_reg_measure.h" -/* *************************************************************** */ /* *************************************************************** */ /// @brief SSD measure of similarity class class reg_ssd: public reg_measure { @@ -40,8 +39,10 @@ class reg_ssd: public reg_measure { nifti_image *voxelBasedGradBw = nullptr) override; /// @brief Define if the specified time point should be normalised void SetNormaliseTimepoint(int timepoint, bool normalise); - /// @brief Returns the ssd value - virtual double GetSimilarityMeasureValue() override; + /// @brief Returns the ssd value forwards + virtual double GetSimilarityMeasureValueFw() override; + /// @brief Returns the ssd value backwards + virtual double GetSimilarityMeasureValueBw() override; /// @brief Compute the voxel based ssd gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override; /// @brief Here @@ -56,12 +57,11 @@ class reg_ssd: public reg_measure { bool normaliseTimePoint[255]; }; /* *************************************************************** */ - /** @brief Computes and returns the SSD between two input images * @param referenceImage First input image to use to compute the metric * @param warpedImage Second input image to use to compute the metric * @param activeTimePoint Specified which time point volumes have to be considered - * @param jacobianDeterminantImage Image that contains the Jacobian + * @param jacobianDetImage Image that contains the Jacobian * determinant of a transformation at every voxel position. This * image is used to modulate the SSD. The argument is ignored if the * pointer is set to nullptr @@ -70,22 +70,22 @@ class reg_ssd: public reg_measure { * @return Returns the computed sum squared difference */ extern "C++" template -double reg_getSSDValue(nifti_image *referenceImage, - nifti_image *warpedImage, - double *timePointWeight, - nifti_image *jacobianDeterminantImage, - int *mask, +double reg_getSsdValue(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const double *timePointWeight, + const nifti_image *jacobianDetImage, + const int *mask, float *currentValue, - nifti_image *localWeightImage); - + const nifti_image *localWeightSim); +/* *************************************************************** */ /** @brief Compute a voxel based gradient of the sum squared difference. * @param referenceImage First input image to use to compute the metric * @param warpedImage Second input image to use to compute the metric * @param activeTimePoint Specified which time point volumes have to be considered - * @param warpedImageGradient Spatial gradient of the input warped image - * @param ssdGradientImage Output image that will be updated with the + * @param warpedGradient Spatial gradient of the input warped image + * @param measureGradientImage Output image that will be updated with the * value of the SSD gradient - * @param jacobianDeterminantImage Image that contains the Jacobian + * @param jacobianDetImage Image that contains the Jacobian * determinant of a transformation at every voxel position. This * image is used to modulate the SSD. The argument is ignored if the * pointer is set to nullptr @@ -93,12 +93,13 @@ double reg_getSSDValue(nifti_image *referenceImage, * should be considered. If set to nullptr, all voxels are considered */ extern "C++" template -void reg_getVoxelBasedSSDGradient(nifti_image *referenceImage, - nifti_image *warpedImage, - nifti_image *warpedImageGradient, - nifti_image *ssdGradientImage, - nifti_image *jacobianDeterminantImage, - int *mask, - int currentTimepoint, - double timepointWeight, - nifti_image *localWeightImage); +void reg_getVoxelBasedSsdGradient(const nifti_image *referenceImage, + const nifti_image *warpedImage, + const nifti_image *warpedGradient, + nifti_image *measureGradientImage, + const nifti_image *jacobianDetImage, + const int *mask, + const int& currentTimepoint, + const double& timepointWeight, + const nifti_image *localWeightSim); +/* *************************************************************** */ diff --git a/reg-lib/cuda/_reg_measure_gpu.h b/reg-lib/cuda/_reg_measure_gpu.h index d91c39d6..1ff52195 100755 --- a/reg-lib/cuda/_reg_measure_gpu.h +++ b/reg-lib/cuda/_reg_measure_gpu.h @@ -99,7 +99,8 @@ class reg_lncc_gpu: public reg_lncc, public reg_measure_gpu { public: /// @brief reg_lncc class constructor reg_lncc_gpu() { - fprintf(stderr, "[ERROR] CUDA CANNOT BE USED WITH LNCC YET\n"); + reg_print_fct_error("reg_lncc_gpu::reg_lncc_gpu"); + reg_print_msg_error("CUDA CANNOT BE USED WITH LNCC YET"); reg_exit(); } /// @brief reg_lncc class destructor @@ -127,8 +128,10 @@ class reg_lncc_gpu: public reg_lncc, public reg_measure_gpu { float4 *warpedGradBwCuda = nullptr, nifti_image *voxelBasedGradBw = nullptr, float4 *voxelBasedGradBwCuda = nullptr) override {} - /// @brief Returns the lncc value - virtual double GetSimilarityMeasureValue() override { return 0; } + /// @brief Returns the lncc value forwards + virtual double GetSimilarityMeasureValueFw() override { return 0; } + /// @brief Returns the lncc value backwards + virtual double GetSimilarityMeasureValueBw() override { return 0; } /// @brief Compute the voxel based lncc gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override {} }; @@ -166,8 +169,10 @@ class reg_kld_gpu: public reg_kld, public reg_measure_gpu { float4 *warpedGradBwCuda = nullptr, nifti_image *voxelBasedGradBw = nullptr, float4 *voxelBasedGradBwCuda = nullptr) override {} - /// @brief Returns the kld value - virtual double GetSimilarityMeasureValue() override { return 0; } + /// @brief Returns the kld value forwards + virtual double GetSimilarityMeasureValueFw() override { return 0; } + /// @brief Returns the kld value backwards + virtual double GetSimilarityMeasureValueBw() override { return 0; } /// @brief Compute the voxel based kld gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override {} }; @@ -205,8 +210,10 @@ class reg_dti_gpu: public reg_dti, public reg_measure_gpu { float4 *warpedGradBwCuda = nullptr, nifti_image *voxelBasedGradBw = nullptr, float4 *voxelBasedGradBwCuda = nullptr) override {} - /// @brief Returns the dti value - virtual double GetSimilarityMeasureValue() override { return 0; } + /// @brief Returns the dti value forwards + virtual double GetSimilarityMeasureValueFw() override { return 0; } + /// @brief Returns the dti value backwards + virtual double GetSimilarityMeasureValueBw() override { return 0; } /// @brief Compute the voxel based dti gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override {} }; diff --git a/reg-lib/cuda/_reg_nmi_gpu.cu b/reg-lib/cuda/_reg_nmi_gpu.cu index 5efd0391..2e55b78b 100755 --- a/reg-lib/cuda/_reg_nmi_gpu.cu +++ b/reg-lib/cuda/_reg_nmi_gpu.cu @@ -63,48 +63,67 @@ void reg_nmi_gpu::InitialiseMeasure(nifti_image *refImg, cudaArray *refImgCuda, #endif } /* *************************************************************** */ -double reg_nmi_gpu::GetSimilarityMeasureValue() { - // The NMI computation is performed into the host for now - // The relevant images have to be transferred from the device to the host - cudaCommon_transferFromDeviceToNifti(this->warpedImage, this->warpedImageCuda); - reg_getNMIValue(this->referenceImage, - this->warpedImage, - this->timePointWeight, - this->referenceBinNumber, - this->floatingBinNumber, - this->totalBinNumber, - this->jointHistogramLog, - this->jointHistogramPro, - this->entropyValues, - this->referenceMask); +double GetSimilarityMeasureValue(const nifti_image *referenceImage, + nifti_image *warpedImage, + const float *warpedImageCuda, + const double *timePointWeight, + const unsigned short *referenceBinNumber, + const unsigned short *floatingBinNumber, + const unsigned short *totalBinNumber, + double **jointHistogramLog, + double **jointHistogramPro, + double **entropyValues, + const int *referenceMask, + const int& referenceTimePoint) { + // The NMI computation is performed on the host for now + cudaCommon_transferFromDeviceToNifti(warpedImage, warpedImageCuda); + reg_getNMIValue(referenceImage, + warpedImage, + timePointWeight, + referenceBinNumber, + floatingBinNumber, + totalBinNumber, + jointHistogramLog, + jointHistogramPro, + entropyValues, + referenceMask); - if (this->isSymmetric) { - cudaCommon_transferFromDeviceToNifti(this->warpedImageBw, this->warpedImageBwCuda); - reg_getNMIValue(this->floatingImage, - this->warpedImageBw, - this->timePointWeight, - this->floatingBinNumber, - this->referenceBinNumber, - this->totalBinNumber, - this->jointHistogramLogBw, - this->jointHistogramProBw, - this->entropyValuesBw, - this->floatingMask); - } - - double nmiFw = 0, nmiBw = 0; - for (int t = 0; t < this->referenceTimePoint; ++t) { - if (this->timePointWeight[t] > 0) { - nmiFw += timePointWeight[t] * (this->entropyValues[t][0] + this->entropyValues[t][1]) / this->entropyValues[t][2]; - if (this->isSymmetric) - nmiBw += timePointWeight[t] * (this->entropyValuesBw[t][0] + this->entropyValuesBw[t][1]) / this->entropyValuesBw[t][2]; - } + double nmi = 0; + for (int t = 0; t < referenceTimePoint; ++t) { + if (timePointWeight[t] > 0) + nmi += timePointWeight[t] * (entropyValues[t][0] + entropyValues[t][1]) / entropyValues[t][2]; } - -#ifndef NDEBUG - reg_print_msg_debug("reg_nmi_gpu::GetSimilarityMeasureValue called"); -#endif - return nmiFw + nmiBw; + return nmi; +} +/* *************************************************************** */ +double reg_nmi_gpu::GetSimilarityMeasureValueFw() { + return ::GetSimilarityMeasureValue(this->referenceImage, + this->warpedImage, + this->warpedImageCuda, + this->timePointWeight, + this->referenceBinNumber, + this->floatingBinNumber, + this->totalBinNumber, + this->jointHistogramLog, + this->jointHistogramPro, + this->entropyValues, + this->referenceMask, + this->referenceTimePoint); +} +/* *************************************************************** */ +double reg_nmi_gpu::GetSimilarityMeasureValueBw() { + return ::GetSimilarityMeasureValue(this->floatingImage, + this->warpedImageBw, + this->warpedImageBwCuda, + this->timePointWeight, + this->floatingBinNumber, + this->referenceBinNumber, + this->totalBinNumber, + this->jointHistogramLogBw, + this->jointHistogramProBw, + this->entropyValuesBw, + this->floatingMask, + this->referenceTimePoint); } /* *************************************************************** */ /// Called when we only have one target and one source image @@ -201,7 +220,7 @@ void reg_nmi_gpu::GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) { this->referenceBinNumber[0]); } #ifndef NDEBUG - reg_print_msg_debug("reg_nmi_gpu::GetVoxelBasedSimilarityMeasureGradient called\n"); + reg_print_msg_debug("reg_nmi_gpu::GetVoxelBasedSimilarityMeasureGradient called"); #endif } /* *************************************************************** */ diff --git a/reg-lib/cuda/_reg_nmi_gpu.h b/reg-lib/cuda/_reg_nmi_gpu.h index ff24a676..2b55270b 100755 --- a/reg-lib/cuda/_reg_nmi_gpu.h +++ b/reg-lib/cuda/_reg_nmi_gpu.h @@ -47,8 +47,10 @@ class reg_nmi_gpu: public reg_nmi, public reg_measure_gpu { float4 *warpedGradBwCuda = nullptr, nifti_image *voxelBasedGradBw = nullptr, float4 *voxelBasedGradBwCuda = nullptr) override; - /// @brief Returns the nmi value - virtual double GetSimilarityMeasureValue() override; + /// @brief Returns the nmi value forwards + virtual double GetSimilarityMeasureValueFw() override; + /// @brief Returns the nmi value backwards + virtual double GetSimilarityMeasureValueBw() override; /// @brief Compute the voxel based nmi gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override; }; @@ -82,8 +84,10 @@ class reg_multichannel_nmi_gpu: public reg_multichannel_nmi, public reg_measure_ reg_multichannel_nmi_gpu() {} /// @brief reg_multichannel_nmi_gpu class destructor virtual ~reg_multichannel_nmi_gpu() {} - /// @brief Returns the nmi value - virtual double GetSimilarityMeasureValue() override { return 0; } + /// @brief Returns the nmi value forwards + virtual double GetSimilarityMeasureValueFw() override { return 0; } + /// @brief Returns the nmi value backwards + virtual double GetSimilarityMeasureValueBw() override { return 0; } /// @brief Compute the voxel based nmi gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override {} }; diff --git a/reg-lib/cuda/_reg_ssd_gpu.cu b/reg-lib/cuda/_reg_ssd_gpu.cu index 1ea2ba08..dc62ea53 100755 --- a/reg-lib/cuda/_reg_ssd_gpu.cu +++ b/reg-lib/cuda/_reg_ssd_gpu.cu @@ -98,7 +98,7 @@ double reg_getSSDValue_gpu(const nifti_image *referenceImage, return ssd; } /* *************************************************************** */ -double reg_ssd_gpu::GetSimilarityMeasureValue() { +double reg_ssd_gpu::GetSimilarityMeasureValueFw() { const double SSDValue = reg_getSSDValue_gpu(this->referenceImage, this->referenceImageCuda, this->warpedImageCuda, @@ -107,6 +107,10 @@ double reg_ssd_gpu::GetSimilarityMeasureValue() { return -SSDValue; } /* *************************************************************** */ +double reg_ssd_gpu::GetSimilarityMeasureValueBw() { + return 0; +} +/* *************************************************************** */ void reg_getVoxelBasedSSDGradient_gpu(const nifti_image *referenceImage, const cudaArray *referenceImageCuda, const float *warpedCuda, diff --git a/reg-lib/cuda/_reg_ssd_gpu.h b/reg-lib/cuda/_reg_ssd_gpu.h index 34764df3..c0a994be 100755 --- a/reg-lib/cuda/_reg_ssd_gpu.h +++ b/reg-lib/cuda/_reg_ssd_gpu.h @@ -48,8 +48,10 @@ class reg_ssd_gpu: public reg_ssd, public reg_measure_gpu { float4 *warpedGradBwCuda = nullptr, nifti_image *voxelBasedGradBw = nullptr, float4 *voxelBasedGradBwCuda = nullptr) override; - /// @brief Returns the ssd value - virtual double GetSimilarityMeasureValue() override; + /// @brief Returns the ssd value forwards + virtual double GetSimilarityMeasureValueFw() override; + /// @brief Returns the ssd value backwards + virtual double GetSimilarityMeasureValueBw() override; /// @brief Compute the voxel based ssd gradient virtual void GetVoxelBasedSimilarityMeasureGradient(int currentTimepoint) override; };