From a148f146a44e5ee30281a89063df1a7be7e35503 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Onur=20=C3=9Clgen?= Date: Tue, 27 Aug 2024 22:42:59 +0100 Subject: [PATCH] Fix Floor(), Ceil(), and Round() functions --- niftyreg_build_version.txt | 2 +- reg-apps/reg_aladin.cpp | 2 +- reg-apps/reg_benchmark.cpp | 60 +++++++-------- reg-apps/reg_f3d.cpp | 2 +- reg-apps/reg_ppcnr.cpp | 2 +- reg-apps/reg_resample.cpp | 2 +- reg-apps/reg_tools.cpp | 16 ++-- reg-apps/reg_transform.cpp | 14 ++-- reg-lib/_reg_base.cpp | 8 +- reg-lib/cpu/Maths.hpp | 22 +++--- reg-lib/cpu/_reg_blockMatching.cpp | 6 +- reg-lib/cpu/_reg_localTrans.cpp | 66 ++++++++-------- reg-lib/cpu/_reg_localTrans_jac.cpp | 40 +++++----- reg-lib/cpu/_reg_localTrans_regul.cpp | 12 +-- reg-lib/cpu/_reg_resampling.cpp | 76 +++++++++---------- reg-lib/cpu/_reg_ssd.cpp | 12 +-- reg-lib/cpu/_reg_tools.cpp | 8 +- reg-lib/cuda/CudaLocalTransformation.cu | 26 +++---- .../cuda/CudaLocalTransformationKernels.cu | 32 ++++---- reg-lib/cuda/CudaResampling.cu | 6 +- reg-lib/cuda/CudaToolsKernels.cu | 2 +- reg-lib/cuda/resampleKernel.cu | 10 +-- .../reg_test_voxelCentricToNodeCentric.cpp | 2 +- 23 files changed, 214 insertions(+), 214 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index e828e5d0..5910394b 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -434 +435 diff --git a/reg-apps/reg_aladin.cpp b/reg-apps/reg_aladin.cpp index 6cf515a4..ff41859b 100755 --- a/reg-apps/reg_aladin.cpp +++ b/reg-apps/reg_aladin.cpp @@ -455,7 +455,7 @@ int main(int argc, char **argv) { time_t end; time(&end); - const int minutes = Floor((end - start) / 60.0f); + const int minutes = Floor((end - start) / 60.0f); const int seconds = static_cast(end - start) - 60 * minutes; NR_VERBOSE_APP("Registration performed in " << minutes << " min " << seconds << " sec"); NR_VERBOSE_APP("Have a good day!"); diff --git a/reg-apps/reg_benchmark.cpp b/reg-apps/reg_benchmark.cpp index dd439f62..8ff35f9a 100644 --- a/reg-apps/reg_benchmark.cpp +++ b/reg-apps/reg_benchmark.cpp @@ -122,9 +122,9 @@ int main(int argc, char **argv) // A control point image is created dim_img[0]=5; - dim_img[1]=Floor(targetImage->nx*targetImage->dx/gridSpacing)+4; - dim_img[2]=Floor(targetImage->ny*targetImage->dy/gridSpacing)+4; - dim_img[3]=Floor(targetImage->nz*targetImage->dz/gridSpacing)+4; + dim_img[1]=Floor(targetImage->nx*targetImage->dx/gridSpacing)+4; + dim_img[2]=Floor(targetImage->ny*targetImage->dy/gridSpacing)+4; + dim_img[3]=Floor(targetImage->nz*targetImage->dz/gridSpacing)+4; dim_img[5]=3; dim_img[4]=dim_img[6]=dim_img[7]=1; nifti_image *controlPointImage = nifti_make_new_nim(dim_img, NIFTI_TYPE_FLOAT32, true); @@ -245,7 +245,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf( "CPU - %i affine deformation field computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i affine deformation field computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -261,7 +261,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i affine deformation field computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i affine deformation field computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -295,7 +295,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i spline deformation field computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i spline deformation field computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -314,7 +314,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i spline deformation field computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i spline deformation field computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -345,7 +345,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i scaling-and-squaring - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i scaling-and-squarings - %i min %i sec\n", maxIt, minutes, seconds); @@ -362,7 +362,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i scaling-and-squaring - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i scaling-and-squarings - %i min %i sec\n", maxIt, minutes, seconds); @@ -395,7 +395,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i linear interpolation computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i linear interpolation computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -416,7 +416,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i linear interpolation computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i linear interpolation computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -447,7 +447,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i spatial gradient computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i spatial gradient computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -466,7 +466,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i spatial gradient computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i spatial gradient computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -525,7 +525,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i voxel-based NMI gradient computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i voxel-based NMI gradient computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -558,7 +558,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i voxel-based NMI gradient computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i voxel-based NMI gradient computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -588,9 +588,9 @@ int main(int argc, char **argv) maxIt=10000 / dimension; // maxIt=1; int smoothingRadius[3]; - smoothingRadius[0] = Floor( 2.0*controlPointImage->dx/targetImage->dx ); - smoothingRadius[1] = Floor( 2.0*controlPointImage->dy/targetImage->dy ); - smoothingRadius[2] = Floor( 2.0*controlPointImage->dz/targetImage->dz ); + smoothingRadius[0] = Floor( 2.0*controlPointImage->dx/targetImage->dx ); + smoothingRadius[1] = Floor( 2.0*controlPointImage->dy/targetImage->dy ); + smoothingRadius[2] = Floor( 2.0*controlPointImage->dz/targetImage->dz ); time(&start); for(int i=0; i(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i node-based NMI gradient computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i node-based NMI gradient computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -620,7 +620,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i node-based NMI gradient computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i node-based NMI gradient computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -650,7 +650,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i BE computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i BE computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -665,7 +665,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i BE computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i BE computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -690,7 +690,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i BE gradient computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i BE gradient computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -708,7 +708,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i BE gradient computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i BE gradient computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -729,7 +729,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i |Jac| penalty term computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i |Jac| penalty term computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -743,7 +743,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i |Jac| penalty term computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i |Jac| penalty term computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -764,7 +764,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i Approx. |Jac| penalty term computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i Approx. |Jac| penalty term computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -778,7 +778,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i Approx. |Jac| penalty term computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i Approx. |Jac| penalty term computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -831,7 +831,7 @@ int main(int argc, char **argv) } time(&end); cpuTime=(end-start); - minutes = Floor(float(cpuTime)/60.0f); + minutes = Floor(float(cpuTime)/60.0f); seconds = (int)(cpuTime - 60*minutes); printf("CPU - %i block matching computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "CPU - %i block matching computations - %i min %i sec\n", maxIt, minutes, seconds); @@ -852,7 +852,7 @@ int main(int argc, char **argv) } time(&end); gpuTime=(end-start); - minutes = Floor(float(gpuTime)/60.0f); + minutes = Floor(float(gpuTime)/60.0f); seconds = (int)(gpuTime - 60*minutes); printf("GPU - %i block matching computations - %i min %i sec\n", maxIt, minutes, seconds); fprintf(outputFile, "GPU - %i block matching computations - %i min %i sec\n", maxIt, minutes, seconds); diff --git a/reg-apps/reg_f3d.cpp b/reg-apps/reg_f3d.cpp index e0b3fe48..e8b16d53 100755 --- a/reg-apps/reg_f3d.cpp +++ b/reg-apps/reg_f3d.cpp @@ -711,7 +711,7 @@ int main(int argc, char **argv) { time_t end; time(&end); - const int minutes = Floor((end - start) / 60.0f); + const int minutes = Floor((end - start) / 60.0f); const int seconds = static_cast(end - start) - 60 * minutes; NR_VERBOSE_APP("Registration performed in " << minutes << " min " << seconds << " sec"); NR_VERBOSE_APP("Have a good day!"); diff --git a/reg-apps/reg_ppcnr.cpp b/reg-apps/reg_ppcnr.cpp index e3c664f3..d3ad209a 100755 --- a/reg-apps/reg_ppcnr.cpp +++ b/reg-apps/reg_ppcnr.cpp @@ -963,7 +963,7 @@ int main(int argc, char **argv) time_t end; time( &end ); - int minutes = Floor(float(end-start)/60.0f); + int minutes = Floor(float(end-start)/60.0f); int seconds = (int)(end-start - 60*minutes); NR_COUT << "* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n"; if(flag->locality) diff --git a/reg-apps/reg_resample.cpp b/reg-apps/reg_resample.cpp index bfcfe963..fc2590ed 100755 --- a/reg-apps/reg_resample.cpp +++ b/reg-apps/reg_resample.cpp @@ -483,7 +483,7 @@ int main(int argc, char **argv) param->interpolation, param->paddingValue, jacobian, - (char)Round(param->PSF_Algorithm)); + Round(param->PSF_Algorithm)); NR_DEBUG("PSF resampling completed"); free(jacobian); } diff --git a/reg-apps/reg_tools.cpp b/reg-apps/reg_tools.cpp index 033077e5..a4172fa0 100755 --- a/reg-apps/reg_tools.cpp +++ b/reg-apps/reg_tools.cpp @@ -499,8 +499,8 @@ int main(int argc, char **argv) reg_tools_changeDatatype(image); nifti_image *normImage = nifti_dup(*image); HeapSort(static_cast(normImage->data), normImage->nvox); - float minValue = static_cast(normImage->data)[Floor(03*(int)normImage->nvox/100)]; - float maxValue = static_cast(normImage->data)[Floor(97*(int)normImage->nvox/100)]; + float minValue = static_cast(normImage->data)[Floor(0.03*normImage->nvox)]; + float maxValue = static_cast(normImage->data)[Floor(0.97*normImage->nvox)]; reg_tools_subtractValueFromImage(image,normImage,minValue); reg_tools_divideValueToImage(normImage,normImage,maxValue-minValue); if(flag->outputImageFlag) @@ -803,10 +803,10 @@ int main(int argc, char **argv) // Define the size of the new image int newDim[8]; for(size_t i=0; i<8; ++i) newDim[i]=image->dim[i]; - newDim[1]=Ceil((float)image->dim[1]*image->pixdim[1]/param->pixdimX); - newDim[2]=Ceil((float)image->dim[2]*image->pixdim[2]/param->pixdimY); + newDim[1]=Ceil((float)image->dim[1]*image->pixdim[1]/param->pixdimX); + newDim[2]=Ceil((float)image->dim[2]*image->pixdim[2]/param->pixdimY); if(image->nz>1) - newDim[3]=Ceil((float)image->dim[3]*image->pixdim[3]/param->pixdimZ); + newDim[3]=Ceil((float)image->dim[3]*image->pixdim[3]/param->pixdimZ); // Create the new image nifti_image *newImg=nifti_make_new_nim(newDim,image->datatype,true); newImg->pixdim[1]=newImg->dx=param->pixdimX; @@ -954,7 +954,7 @@ int main(int argc, char **argv) for(int y=0; yny; ++y){ for(int x=0; xnx; ++x){ size_t outIndex = ((z*image->ny+y)*image->nx+x)*image->nt*image->nu+t; - outPtr[outIndex] = Round(*inPtr); + outPtr[outIndex] = Round(*inPtr); ++inPtr; } } @@ -997,8 +997,8 @@ int main(int argc, char **argv) float value = *inPtr * 255.f; size_t outIndex = ((z*image->ny+y)*image->nx+x)*3; if (value > 0) - outPtr[outIndex] = static_cast(Round(value>255?255:value)); - else outPtr[outIndex+1] = static_cast(Round(-value<-255?-255:-value)); + outPtr[outIndex] = Round(value > 255 ? 255 : value); + else outPtr[outIndex + 1] = Round(-value < -255 ? -255 : -value); outPtr[outIndex+2] = 0; ++inPtr; } diff --git a/reg-apps/reg_transform.cpp b/reg-apps/reg_transform.cpp index cb19fd66..5b796154 100755 --- a/reg-apps/reg_transform.cpp +++ b/reg-apps/reg_transform.cpp @@ -398,7 +398,7 @@ int main(int argc, char **argv) { if (affineTransformation != nullptr) { reg_affine_getDeformationField(affineTransformation, outputTransformationImage); } else { - switch (Round(inputTransformationImage->intent_p1)) { + switch (Round(inputTransformationImage->intent_p1)) { case DEF_FIELD: NR_INFO("The specified transformation is a deformation field:"); NR_INFO(inputTransformationImage->fname); @@ -468,7 +468,7 @@ int main(int argc, char **argv) { } // Save the generated transformation reg_io_WriteImageFile(outputTransformationImage, param->outputTransName); - switch (Round(outputTransformationImage->intent_p1)) { + switch (Round(outputTransformationImage->intent_p1)) { case DEF_FIELD: NR_INFO("The deformation field has been saved as:"); NR_INFO(param->outputTransName); @@ -593,7 +593,7 @@ int main(int argc, char **argv) { output1TransImage->data = calloc(output1TransImage->nvox, output1TransImage->nbyper); if (affine1Trans != nullptr) { reg_affine_getDeformationField(affine1Trans, output1TransImage); - } else switch (Round(input1TransImage->intent_p1)) { + } else switch (Round(input1TransImage->intent_p1)) { case LIN_SPLINE_GRID: case CUB_SPLINE_GRID: NR_INFO("Transformation 1 is a spline parametrisation:"); @@ -659,7 +659,7 @@ int main(int argc, char **argv) { reg_affine_getDeformationField(affine2Trans, output2TransImage); reg_defField_compose(output2TransImage, output1TransImage, nullptr); } else { - switch (Round(input2TransImage->intent_p1)) { + switch (Round(input2TransImage->intent_p1)) { case LIN_SPLINE_GRID: case CUB_SPLINE_GRID: NR_INFO("Transformation 2 is a spline parametrisation:"); @@ -821,7 +821,7 @@ int main(int argc, char **argv) { if (affineTransformation != nullptr) { reg_affine_getDeformationField(affineTransformation, deformationFieldImage); } else if (inputTransformationImage != nullptr) { - switch (Round(inputTransformationImage->intent_p1)) { + switch (Round(inputTransformationImage->intent_p1)) { case DEF_FIELD: NR_INFO("The specified transformation is a deformation field:"); NR_INFO(inputTransformationImage->fname); @@ -982,7 +982,7 @@ int main(int argc, char **argv) { NR_ERROR("Error when reading the input image: " << param->inputTransName); return EXIT_FAILURE; } - switch (Round(inputTransImage->intent_p1)) { + switch (Round(inputTransImage->intent_p1)) { case LIN_SPLINE_GRID: case CUB_SPLINE_GRID: reg_getDisplacementFromDeformation(inputTransImage); @@ -1099,7 +1099,7 @@ int main(int argc, char **argv) { outputTransImage->scl_inter = 0.f; outputTransImage->data = malloc(outputTransImage->nvox * outputTransImage->nbyper); // Invert the provided - switch (Round(inputTransImage->intent_p1)) { + switch (Round(inputTransImage->intent_p1)) { case DEF_FIELD: reg_defFieldInvert(inputTransImage, outputTransImage, 1.0e-6f); memset(outputTransImage->descrip, 0, 80); diff --git a/reg-lib/_reg_base.cpp b/reg-lib/_reg_base.cpp index 4eb441ef..d5c7624f 100644 --- a/reg-lib/_reg_base.cpp +++ b/reg-lib/_reg_base.cpp @@ -405,9 +405,9 @@ void reg_base::Initialise() { HeapSort(refDataPtr, tmpReference->nvox); // Update the reference threshold values if no value has been setup by the user if (referenceThresholdLow[0] == std::numeric_limits::lowest()) - referenceThresholdLow[0] = refDataPtr[Round((float)tmpReference->nvox * 0.02f)]; + referenceThresholdLow[0] = refDataPtr[Round(tmpReference->nvox * 0.02)]; if (referenceThresholdUp[0] == std::numeric_limits::max()) - referenceThresholdUp[0] = refDataPtr[Round((float)tmpReference->nvox * 0.98f)]; + referenceThresholdUp[0] = refDataPtr[Round(tmpReference->nvox * 0.98)]; // Create a copy of the floating image to extract the robust range NiftiImage tmpFloating = inputFloating; @@ -417,9 +417,9 @@ void reg_base::Initialise() { HeapSort(floDataPtr, tmpFloating->nvox); // Update the floating threshold values if no value has been setup by the user if (floatingThresholdLow[0] == std::numeric_limits::lowest()) - floatingThresholdLow[0] = floDataPtr[Round((float)tmpFloating->nvox * 0.02f)]; + floatingThresholdLow[0] = floDataPtr[Round(tmpFloating->nvox * 0.02)]; if (floatingThresholdUp[0] == std::numeric_limits::max()) - floatingThresholdUp[0] = floDataPtr[Round((float)tmpFloating->nvox * 0.98f)]; + floatingThresholdUp[0] = floDataPtr[Round(tmpFloating->nvox * 0.98)]; } // FINEST LEVEL OF REGISTRATION diff --git a/reg-lib/cpu/Maths.hpp b/reg-lib/cpu/Maths.hpp index 56782eda..58dbe38a 100644 --- a/reg-lib/cpu/Maths.hpp +++ b/reg-lib/cpu/Maths.hpp @@ -59,19 +59,19 @@ template DEVICE inline T Cube(const T& x) { return x * x * x; } -template -DEVICE inline int Floor(const T& x) { - const int i = static_cast(x); - return i - (x < i); +template +DEVICE inline RetT Floor(const T& x) { + const int64_t i = static_cast(x); + return static_cast(i - (x < i)); } -template -DEVICE inline int Ceil(const T& x) { - const int i = static_cast(x); - return i + (x > i); +template +DEVICE inline RetT Ceil(const T& x) { + const int64_t i = static_cast(x); + return static_cast(i + (x > i)); } -template -DEVICE inline int Round(const T& x) { - return static_cast(x + (x >= 0 ? 0.5 : -0.5)); +template +DEVICE inline RetT Round(const T& x) { + return static_cast(static_cast(x + (x >= 0 ? 0.5 : -0.5))); } /* *************************************************************** */ DEVICE inline void Divide(const int num, const int denom, int& quot, int& rem) { diff --git a/reg-lib/cpu/_reg_blockMatching.cpp b/reg-lib/cpu/_reg_blockMatching.cpp index fce081f5..7f48c1f2 100755 --- a/reg-lib/cpu/_reg_blockMatching.cpp +++ b/reg-lib/cpu/_reg_blockMatching.cpp @@ -232,10 +232,10 @@ void initialise_block_matching_method(nifti_image * reference, } params->voxelCaptureRange = 3; - params->blockNumber[0] = Ceil((double)reference->nx / (double)BLOCK_WIDTH); - params->blockNumber[1] = Ceil((double)reference->ny / (double)BLOCK_WIDTH); + params->blockNumber[0] = Ceil((double)reference->nx / (double)BLOCK_WIDTH); + params->blockNumber[1] = Ceil((double)reference->ny / (double)BLOCK_WIDTH); if (reference->nz > 1) { - params->blockNumber[2] = Ceil((double)reference->nz / (double)BLOCK_WIDTH); + params->blockNumber[2] = Ceil((double)reference->nz / (double)BLOCK_WIDTH); params->dim = 3; } else { diff --git a/reg-lib/cpu/_reg_localTrans.cpp b/reg-lib/cpu/_reg_localTrans.cpp index 418b310d..0939a91e 100755 --- a/reg-lib/cpu/_reg_localTrans.cpp +++ b/reg-lib/cpu/_reg_localTrans.cpp @@ -24,9 +24,9 @@ void reg_createControlPointGrid(NiftiImage& controlPointGridImage, const float *spacing) { // Define the control point grid dimensions vector dims{ - Ceil(referenceImage->nx * referenceImage->dx / spacing[0] + 3.f), - Ceil(referenceImage->ny * referenceImage->dy / spacing[1] + 3.f), - referenceImage->nz > 1 ? Ceil(referenceImage->nz * referenceImage->dz / spacing[2] + 3.f) : 1, + Ceil(referenceImage->nx * referenceImage->dx / spacing[0] + 3.f), + Ceil(referenceImage->ny * referenceImage->dy / spacing[1] + 3.f), + referenceImage->nz > 1 ? Ceil(referenceImage->nz * referenceImage->dz / spacing[2] + 3.f) : 1, 1, referenceImage->nz > 1 ? 3 : 2 }; @@ -277,9 +277,9 @@ void reg_createSymmetricControlPointGrids(NiftiImage& forwardGridImage, // Compute the dimension of the control point grids const vector dims{ - Ceil((maxPosition[0] - minPosition[0]) / spacing[0] + 3.f), - Ceil((maxPosition[1] - minPosition[1]) / spacing[1] + 3.f), - referenceImage->nz > 1 ? Ceil((maxPosition[2] - minPosition[2]) / spacing[2] + 3.f) : 1, + Ceil((maxPosition[0] - minPosition[0]) / spacing[0] + 3.f), + Ceil((maxPosition[1] - minPosition[1]) / spacing[1] + 3.f), + referenceImage->nz > 1 ? Ceil((maxPosition[2] - minPosition[2]) / spacing[2] + 3.f) : 1, 1, referenceImage->nz > 1 ? 3 : 2 }; @@ -452,17 +452,17 @@ void reg_linear_spline_getDeformationField3D(nifti_image *splineControlPoint, referenceMatrix_real_to_voxel.m[2][3]; // The spline coefficients are computed - xPre = Floor(voxel[0]); + xPre = Floor(voxel[0]); xBasis[1] = voxel[0] - static_cast(xPre); if (xBasis[1] < 0) xBasis[1] = 0; //rounding error xBasis[0] = 1.f - xBasis[1]; - yPre = Floor(voxel[1]); + yPre = Floor(voxel[1]); yBasis[1] = voxel[1] - static_cast(yPre); if (yBasis[1] < 0) yBasis[1] = 0; //rounding error yBasis[0] = 1.f - yBasis[1]; - zPre = Floor(voxel[2]); + zPre = Floor(voxel[2]); zBasis[1] = voxel[2] - static_cast(zPre); if (zBasis[1] < 0) zBasis[1] = 0; //rounding error zBasis[0] = 1.f - zBasis[1]; @@ -643,13 +643,13 @@ void reg_cubic_spline_getDeformationField2D(nifti_image *splineControlPoint, + referenceMatrix_real_to_voxel->m[1][3]; // The spline coefficients are computed - xPre = Floor(xVoxel); + xPre = Floor(xVoxel); basis = xVoxel - static_cast(xPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, xBasis); else get_SplineBasisValues(basis, xBasis); - yPre = Floor(yVoxel); + yPre = Floor(yVoxel); basis = yVoxel - static_cast(yPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, yBasis); @@ -970,19 +970,19 @@ void reg_cubic_spline_getDeformationField3D(nifti_image *splineControlPoint, referenceMatrix_real_to_voxel.m[2][3]; // The spline coefficients are computed - xPre = Floor(voxel[0]); + xPre = Floor(voxel[0]); basis = voxel[0] - static_cast(xPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, xBasis); else get_SplineBasisValues(basis, xBasis); - yPre = Floor(voxel[1]); + yPre = Floor(voxel[1]); basis = voxel[1] - static_cast(yPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, yBasis); else get_SplineBasisValues(basis, yBasis); - zPre = Floor(voxel[2]); + zPre = Floor(voxel[2]); basis = voxel[2] - static_cast(zPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, zBasis); @@ -1614,9 +1614,9 @@ void reg_voxelCentricToNodeCentric(nifti_image *nodeImage, // linear interpolation is performed DataType basisX[2], basisY[2], basisZ[2] = { 0, 0 }; int pre[3] = { - Floor(voxelCoord[0]), - Floor(voxelCoord[1]), - Floor(voxelCoord[2]) + Floor(voxelCoord[0]), + Floor(voxelCoord[1]), + Floor(voxelCoord[2]) }; basisX[1] = voxelCoord[0] - static_cast(pre[0]); basisX[0] = static_cast(1) - basisX[1]; @@ -1731,8 +1731,8 @@ void reg_spline_refineControlPointGrid2D(nifti_image *splineControlPoint, splineControlPoint->dy = splineControlPoint->pixdim[2] = splineControlPoint->dy / 2.0f; splineControlPoint->dz = 1.0f; if (referenceImage != nullptr) { - splineControlPoint->dim[1] = splineControlPoint->nx = Ceil(referenceImage->nx * referenceImage->dx / splineControlPoint->dx + 3.f); - splineControlPoint->dim[2] = splineControlPoint->ny = Ceil(referenceImage->ny * referenceImage->dy / splineControlPoint->dy + 3.f); + splineControlPoint->dim[1] = splineControlPoint->nx = Ceil(referenceImage->nx * referenceImage->dx / splineControlPoint->dx + 3.f); + splineControlPoint->dim[2] = splineControlPoint->ny = Ceil(referenceImage->ny * referenceImage->dy / splineControlPoint->dy + 3.f); } else { splineControlPoint->dim[1] = splineControlPoint->nx = (oldDim[1] - 3) * 2 + 3; splineControlPoint->dim[2] = splineControlPoint->ny = (oldDim[2] - 3) * 2 + 3; @@ -1819,9 +1819,9 @@ void reg_spline_refineControlPointGrid3D(nifti_image *splineControlPoint, nifti_ splineControlPoint->dz = splineControlPoint->pixdim[3] = splineControlPoint->dz / 2.0f; if (referenceImage != nullptr) { - splineControlPoint->dim[1] = splineControlPoint->nx = Ceil(referenceImage->nx * referenceImage->dx / splineControlPoint->dx + 3.f); - splineControlPoint->dim[2] = splineControlPoint->ny = Ceil(referenceImage->ny * referenceImage->dy / splineControlPoint->dy + 3.f); - splineControlPoint->dim[3] = splineControlPoint->nz = Ceil(referenceImage->nz * referenceImage->dz / splineControlPoint->dz + 3.f); + splineControlPoint->dim[1] = splineControlPoint->nx = Ceil(referenceImage->nx * referenceImage->dx / splineControlPoint->dx + 3.f); + splineControlPoint->dim[2] = splineControlPoint->ny = Ceil(referenceImage->ny * referenceImage->dy / splineControlPoint->dy + 3.f); + splineControlPoint->dim[3] = splineControlPoint->nz = Ceil(referenceImage->nz * referenceImage->dz / splineControlPoint->dz + 3.f); } else { splineControlPoint->dim[1] = splineControlPoint->nx = (oldDim[1] - 3) * 2 + 3; splineControlPoint->dim[2] = splineControlPoint->ny = (oldDim[2] - 3) * 2 + 3; @@ -2292,8 +2292,8 @@ void reg_defField_compose2D(const nifti_image *deformationField, df_real2Voxel->m[1][3]; // Linear interpolation to compute the new deformation - pre[0] = Floor(voxelX); - pre[1] = Floor(voxelY); + pre[0] = Floor(voxelX); + pre[1] = Floor(voxelY); relX[1] = voxelX - static_cast(pre[0]); relX[0] = 1.f - relX[1]; relY[1] = voxelY - static_cast(pre[1]); @@ -2399,9 +2399,9 @@ void reg_defField_compose3D(const nifti_image *deformationField, df_real2Voxel.m[2][3]; // Linear interpolation to compute the new deformation - pre[0] = Floor(voxel[0]); - pre[1] = Floor(voxel[1]); - pre[2] = Floor(voxel[2]); + pre[0] = Floor(voxel[0]); + pre[1] = Floor(voxel[1]); + pre[2] = Floor(voxel[2]); relX[1] = voxel[0] - static_cast(pre[0]); relX[0] = 1.f - relX[1]; relY[1] = voxel[1] - static_cast(pre[1]); @@ -3110,13 +3110,13 @@ void reg_spline_cppComposition_2D(nifti_image *grid1, + matrix_real_to_voxel1->m[1][3]; // The spline coefficients are computed - int xPre = Floor(xVoxel); + int xPre = Floor(xVoxel); basis = xVoxel - static_cast(xPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, xBasis); else get_SplineBasisValues(basis, xBasis); - int yPre = Floor(yVoxel); + int yPre = Floor(yVoxel); basis = yVoxel - static_cast(yPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, yBasis); @@ -3322,19 +3322,19 @@ void reg_spline_cppComposition_3D(nifti_image *grid1, + matrix_real_to_voxel1->m[2][3]; // The spline coefficients are computed - xPre = Floor(xVoxel); + xPre = Floor(xVoxel); basis = xVoxel - static_cast(xPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, xBasis); else get_SplineBasisValues(basis, xBasis); - yPre = Floor(yVoxel); + yPre = Floor(yVoxel); basis = yVoxel - static_cast(yPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, yBasis); else get_SplineBasisValues(basis, yBasis); - zPre = Floor(zVoxel); + zPre = Floor(zVoxel); basis = zVoxel - static_cast(zPre--); if (basis < 0) basis = 0; //rounding error if (bspline) get_BSplineBasisValues(basis, zBasis); @@ -3528,7 +3528,7 @@ void reg_defField_getDeformationFieldFromFlowField(nifti_image *flowField, squaringNumber = squaringNumber < 6 ? 6 : squaringNumber; // Set the number of squaring step in the flow field if (fabs(flowField->intent_p2) != squaringNumber) { - NR_WARN("Changing from " << Round(fabs(flowField->intent_p2)) << " to " << abs(squaringNumber) << + NR_WARN("Changing from " << Round(fabs(flowField->intent_p2)) << " to " << abs(squaringNumber) << " squaring step (equivalent to scaling down by " << (int)pow(2.0f, squaringNumber) << ")"); } // Update the number of squaring step required diff --git a/reg-lib/cpu/_reg_localTrans_jac.cpp b/reg-lib/cpu/_reg_localTrans_jac.cpp index 303057cb..e023ed9f 100755 --- a/reg-lib/cpu/_reg_localTrans_jac.cpp +++ b/reg-lib/cpu/_reg_localTrans_jac.cpp @@ -163,9 +163,9 @@ void reg_linear_spline_jacobian3D(nifti_image *splineControlPoint, // Compute the position in the grid Mat44Mul(transformation,imageCoord,gridCoord); // Compute the anterior node coord - pre[0]=Floor(gridCoord[0]); - pre[1]=Floor(gridCoord[1]); - pre[2]=Floor(gridCoord[2]); + pre[0]=Floor(gridCoord[0]); + pre[1]=Floor(gridCoord[1]); + pre[2]=Floor(gridCoord[2]); int controlPoint_index=(pre[2]*splineControlPoint->ny+pre[1])*splineControlPoint->nx+pre[0]; jacobianMatrix.m[0][0] = (coeffPtrX[controlPoint_index+1] - coeffPtrX[controlPoint_index]); @@ -378,8 +378,8 @@ void reg_cubic_spline_jacobian2D(nifti_image *splineControlPoint, // Compute the position in the grid Mat44Mul(transformation,imageCoord,gridCoord); // Compute the anterior node coord - pre[0]=Floor(gridCoord[0]); - pre[1]=Floor(gridCoord[1]); + pre[0]=Floor(gridCoord[0]); + pre[1]=Floor(gridCoord[1]); // Compute the basis values and their first derivatives basis = gridCoord[0] - pre[0]; get_BSplineBasisValues(basis, xBasis, xFirst); @@ -788,9 +788,9 @@ void reg_cubic_spline_jacobian3D(nifti_image *splineControlPoint, // Compute the position in the grid Mat44Mul(transformation,imageCoord,gridCoord); // Compute the anterior node coord - pre[0]=Floor(gridCoord[0]); - pre[1]=Floor(gridCoord[1]); - pre[2]=Floor(gridCoord[2]); + pre[0]=Floor(gridCoord[0]); + pre[1]=Floor(gridCoord[1]); + pre[2]=Floor(gridCoord[2]); // Compute the basis values and their first derivatives basis = gridCoord[0] - pre[0]; get_BSplineBasisValues(basis, xBasis, xFirst); @@ -1476,7 +1476,7 @@ void reg_spline_jacobianDetGradient2D(nifti_image *splineControlPoint, // Loop over all the control points in the surrounding area - for(pixelY=Ceil((y-3)*gridVoxelSpacing[1]); pixelY<=Ceil((y+1)*gridVoxelSpacing[1]); pixelY++) + for(pixelY=Ceil((y-3)*gridVoxelSpacing[1]); pixelY<=Ceil((y+1)*gridVoxelSpacing[1]); pixelY++) { if(pixelY>-1 && pixelYny) { @@ -1485,9 +1485,9 @@ void reg_spline_jacobianDetGradient2D(nifti_image *splineControlPoint, basis=(DataType)pixelY/gridVoxelSpacing[1]-(DataType)yPre; get_BSplineBasisValue(basis,y-yPre,yBasis,yFirst); - jacIndex = pixelY*referenceImage->nx+Ceil((x-3)*gridVoxelSpacing[0]); + jacIndex = pixelY*referenceImage->nx+Ceil((x-3)*gridVoxelSpacing[0]); - for(pixelX=Ceil((x-3)*gridVoxelSpacing[0]); pixelX<=Ceil((x+1)*gridVoxelSpacing[0]); pixelX++) + for(pixelX=Ceil((x-3)*gridVoxelSpacing[0]); pixelX<=Ceil((x+1)*gridVoxelSpacing[0]); pixelX++) { if(pixelX>-1 && pixelXnx && (yFirst!=0 || yBasis!=0)) { @@ -1740,7 +1740,7 @@ void reg_spline_jacobianDetGradient3D(nifti_image *splineControlPoint, jacobianConstraint[0]=jacobianConstraint[1]=jacobianConstraint[2]=0.; // Loop over all the control points in the surrounding area - for(pixelZ=Ceil((z-3)*gridVoxelSpacing[2]); pixelZ<=Ceil((z+1)*gridVoxelSpacing[2]); pixelZ++) + for(pixelZ=Ceil((z-3)*gridVoxelSpacing[2]); pixelZ<=Ceil((z+1)*gridVoxelSpacing[2]); pixelZ++) { if(pixelZ>-1 && pixelZnz) { @@ -1749,7 +1749,7 @@ void reg_spline_jacobianDetGradient3D(nifti_image *splineControlPoint, basis=(DataType)pixelZ/gridVoxelSpacing[2]-(DataType)zPre; get_BSplineBasisValue(basis,z-zPre,zBasis,zFirst); - for(pixelY=Ceil((y-3)*gridVoxelSpacing[1]); pixelY<=Ceil((y+1)*gridVoxelSpacing[1]); pixelY++) + for(pixelY=Ceil((y-3)*gridVoxelSpacing[1]); pixelY<=Ceil((y+1)*gridVoxelSpacing[1]); pixelY++) { if(pixelY>-1 && pixelYny && (zFirst!=0 || zBasis!=0)) { @@ -1758,9 +1758,9 @@ void reg_spline_jacobianDetGradient3D(nifti_image *splineControlPoint, basis=(DataType)pixelY/gridVoxelSpacing[1]-(DataType)yPre; get_BSplineBasisValue(basis,y-yPre,yBasis,yFirst); - jacIndex = (pixelZ*referenceImage->ny+pixelY)*referenceImage->nx+Ceil((x-3)*gridVoxelSpacing[0]); + jacIndex = (pixelZ*referenceImage->ny+pixelY)*referenceImage->nx+Ceil((x-3)*gridVoxelSpacing[0]); - for(pixelX=Ceil((x-3)*gridVoxelSpacing[0]); pixelX<=Ceil((x+1)*gridVoxelSpacing[0]); pixelX++) + for(pixelX=Ceil((x-3)*gridVoxelSpacing[0]); pixelX<=Ceil((x+1)*gridVoxelSpacing[0]); pixelX++) { if(pixelX>-1 && pixelXnx && (yFirst!=0 || yBasis!=0)) { @@ -2068,12 +2068,12 @@ double reg_spline_correctFolding2D(nifti_image *splineControlPoint, // Loop over all the control points in the surrounding area - for(pixelY=Ceil((y-3)*gridVoxelSpacing[1]); pixelY((y-3)*gridVoxelSpacing[1]); pixelY((y+1)*gridVoxelSpacing[1]); pixelY++) { if(pixelY>-1 && pixelYny) { - for(pixelX=Ceil((x-3)*gridVoxelSpacing[0]); pixelX((x-3)*gridVoxelSpacing[0]); pixelX((x+1)*gridVoxelSpacing[0]); pixelX++) { if(pixelX>-1 && pixelXnx) { @@ -2340,17 +2340,17 @@ double reg_spline_correctFolding3D(nifti_image *splineControlPoint, correctFolding=false; // Loop over all the control points in the surrounding area - for(pixelZ=Ceil((z-3)*gridVoxelSpacing[2]); pixelZ((z-3)*gridVoxelSpacing[2]); pixelZ((z+1)*gridVoxelSpacing[2]); pixelZ++) { if(pixelZ>-1 && pixelZnz) { - for(pixelY=Ceil((y-3)*gridVoxelSpacing[1]); pixelY((y-3)*gridVoxelSpacing[1]); pixelY((y+1)*gridVoxelSpacing[1]); pixelY++) { if(pixelY>-1 && pixelYny) { - for(pixelX=Ceil((x-3)*gridVoxelSpacing[0]); pixelX((x-3)*gridVoxelSpacing[0]); pixelX((x+1)*gridVoxelSpacing[0]); pixelX++) { if(pixelX>-1 && pixelXnx) { diff --git a/reg-lib/cpu/_reg_localTrans_regul.cpp b/reg-lib/cpu/_reg_localTrans_regul.cpp index 907b3f3b..d1175e4b 100755 --- a/reg-lib/cpu/_reg_localTrans_regul.cpp +++ b/reg-lib/cpu/_reg_localTrans_regul.cpp @@ -892,9 +892,9 @@ double reg_spline_getLandmarkDistance_core(const nifti_image *controlPointImage, Mat44Mul(*gridRealToVox, refPosition, defPosition); // Extract the corresponding nodes - previous[0] = Floor(defPosition[0]) - 1; - previous[1] = Floor(defPosition[1]) - 1; - previous[2] = Floor(defPosition[2]) - 1; + previous[0] = Floor(defPosition[0]) - 1; + previous[1] = Floor(defPosition[1]) - 1; + previous[2] = Floor(defPosition[2]) - 1; // Check that the specified landmark belongs to the input image if (previous[0] > -1 && previous[0] + 3 < controlPointImage->nx && previous[1] > -1 && previous[1] + 3 < controlPointImage->ny && @@ -1005,9 +1005,9 @@ void reg_spline_getLandmarkDistanceGradient_core(const nifti_image *controlPoint Mat44Mul(*gridRealToVox, refPosition, defPosition); if (imageDim == 2) defPosition[2] = 0; // Extract the corresponding nodes - previous[0] = Floor(defPosition[0]) - 1; - previous[1] = Floor(defPosition[1]) - 1; - previous[2] = Floor(defPosition[2]) - 1; + previous[0] = Floor(defPosition[0]) - 1; + previous[1] = Floor(defPosition[1]) - 1; + previous[2] = Floor(defPosition[2]) - 1; // Check that the specified landmark belongs to the input image if (previous[0] > -1 && previous[0] + 3 < controlPointImage->nx && previous[1] > -1 && previous[1] + 3 < controlPointImage->ny && diff --git a/reg-lib/cpu/_reg_resampling.cpp b/reg-lib/cpu/_reg_resampling.cpp index 4a6f9447..dec53254 100755 --- a/reg-lib/cpu/_reg_resampling.cpp +++ b/reg-lib/cpu/_reg_resampling.cpp @@ -397,9 +397,9 @@ void ResampleImage3D(const nifti_image *floatingImage, // real -> voxel; floating space Mat44Mul(*floatingIJKMatrix, world, position); - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); - previous[2] = Floor(position[2]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); + previous[2] = Floor(position[2]); relative[0] = static_cast(position[0]) - static_cast(previous[0]); relative[1] = static_cast(position[1]) - static_cast(previous[1]); @@ -547,8 +547,8 @@ void ResampleImage2D(const nifti_image *floatingImage, // real -> voxel; floating space Mat44Mul(*floatingIJKMatrix, world, position); - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); relative[0] = static_cast(position[0]) - static_cast(previous[0]); relative[1] = static_cast(position[1]) - static_cast(previous[1]); @@ -818,13 +818,13 @@ void ResampleImage3D_PSF_Sinc(const nifti_image *floatingImage, // Interpolate (trilinearly) the deformation field for non-integer positions float scalling = 1.0f; - currentAPre = (float)Floor(currentA + (shiftSamp[0] / warpedImage->pixdim[1]) * scalling); + currentAPre = Floor(currentA + (shiftSamp[0] / warpedImage->pixdim[1]) * scalling); currentARel = currentA + (shiftSamp[0] / warpedImage->pixdim[1] * scalling) - (float)(currentAPre); - currentBPre = (float)Floor(currentB + (shiftSamp[1] / warpedImage->pixdim[2])); + currentBPre = Floor(currentB + (shiftSamp[1] / warpedImage->pixdim[2])); currentBRel = currentB + (shiftSamp[1] / warpedImage->pixdim[2] * scalling) - (float)(currentBPre); - currentCPre = (float)Floor(currentC + (shiftSamp[2] / warpedImage->pixdim[3] * scalling)); + currentCPre = Floor(currentC + (shiftSamp[2] / warpedImage->pixdim[3] * scalling)); currentCRel = currentC + (shiftSamp[2] / warpedImage->pixdim[3] * scalling) - (float)(currentCPre); // Interpolate the PSF world coordinates @@ -870,9 +870,9 @@ void ResampleImage3D_PSF_Sinc(const nifti_image *floatingImage, // real -> voxel; floating space Mat44Mul(*floatingIJKMatrix, psfWorld, position); - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); - previous[2] = Floor(position[2]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); + previous[2] = Floor(position[2]); relative[0] = position[0] - static_cast(previous[0]); relative[1] = position[1] - static_cast(previous[1]); @@ -1191,13 +1191,13 @@ void ResampleImage3D_PSF(const nifti_image *floatingImage, if (psfWeight != 0.f) { // If the relative weight is above 0 // Interpolate (trilinearly) the deformation field for non-integer positions - currentAPre = (size_t)(currentA + (size_t)Floor(psf_xyz[0] / (float)warpedImage->pixdim[1])); + currentAPre = currentA + Floor(psf_xyz[0] / (float)warpedImage->pixdim[1]); currentARel = (float)currentA + (float)(psf_xyz[0] / (float)warpedImage->pixdim[1]) - (float)(currentAPre); - currentBPre = (size_t)(currentB + (size_t)Floor(psf_xyz[1] / (float)warpedImage->pixdim[2])); + currentBPre = currentB + Floor(psf_xyz[1] / (float)warpedImage->pixdim[2]); currentBRel = (float)currentB + (float)(psf_xyz[1] / (float)warpedImage->pixdim[2]) - (float)(currentBPre); - currentCPre = (size_t)(currentC + (size_t)Floor(psf_xyz[2] / (float)warpedImage->pixdim[3])); + currentCPre = currentC + Floor(psf_xyz[2] / (float)warpedImage->pixdim[3]); currentCRel = (float)currentC + (float)(psf_xyz[2] / (float)warpedImage->pixdim[3]) - (float)(currentCPre); // Interpolate the PSF world coordinates @@ -1242,9 +1242,9 @@ void ResampleImage3D_PSF(const nifti_image *floatingImage, // real -> voxel; floating space Mat44Mul(*floatingIJKMatrix, psfWorld, position); - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); - previous[2] = Floor(position[2]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); + previous[2] = Floor(position[2]); relative[0] = position[0] - static_cast(previous[0]); relative[1] = position[1] - static_cast(previous[1]); @@ -1453,10 +1453,10 @@ void reg_bilinearResampleGradient(const nifti_image *floatingImage, floating_mm_to_voxel->m[1][3]; // Extract the floating value using bilinear interpolation - anteIntX[0] = Floor(xFloCoord); - anteIntX[1] = Ceil(xFloCoord); - anteIntY[0] = Floor(yFloCoord); - anteIntY[1] = Ceil(yFloCoord); + anteIntX[0] = Floor(xFloCoord); + anteIntX[1] = Ceil(xFloCoord); + anteIntY[0] = Floor(yFloCoord); + anteIntY[1] = Ceil(yFloCoord); val_x = 0; val_y = 0; basisX[1] = fabs(xFloCoord - (DataType)anteIntX[0]); @@ -1633,12 +1633,12 @@ void reg_trilinearResampleGradient(const nifti_image *floatingImage, floating_mm_to_voxel->m[2][3]; // Extract the floating value using bilinear interpolation - anteIntX[0] = Floor(xFloCoord); - anteIntX[1] = Ceil(xFloCoord); - anteIntY[0] = Floor(yFloCoord); - anteIntY[1] = Ceil(yFloCoord); - anteIntZ[0] = Floor(zFloCoord); - anteIntZ[1] = Ceil(zFloCoord); + anteIntX[0] = Floor(xFloCoord); + anteIntX[1] = Ceil(xFloCoord); + anteIntY[0] = Floor(yFloCoord); + anteIntY[1] = Ceil(yFloCoord); + anteIntZ[0] = Floor(zFloCoord); + anteIntZ[1] = Ceil(zFloCoord); val_x = 0; val_y = 0; val_z = 0; @@ -1859,9 +1859,9 @@ void TrilinearImageGradient(const nifti_image *floatingImage, /* real -> voxel; floating space */ Mat44Mul(*floatingIJKMatrix, world, position); - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); - previous[2] = Floor(position[2]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); + previous[2] = Floor(position[2]); // basis values along the x axis relative = position[0] - (FieldType)previous[0]; xBasis[0] = (FieldType)(1.0 - relative); @@ -2026,8 +2026,8 @@ void BilinearImageGradient(const nifti_image *floatingImage, position[0] = world[0] * floatingIJKMatrix->m[0][0] + world[1] * floatingIJKMatrix->m[0][1] + floatingIJKMatrix->m[0][3]; position[1] = world[0] * floatingIJKMatrix->m[1][0] + world[1] * floatingIJKMatrix->m[1][1] + floatingIJKMatrix->m[1][3]; - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); // basis values along the x axis relative = position[0] - (FieldType)previous[0]; relative = relative > 0 ? relative : 0; @@ -2136,9 +2136,9 @@ void CubicSplineImageGradient3D(const nifti_image *floatingImage, /* real -> voxel; floating space */ Mat44Mul(*floatingIJKMatrix, world, position); - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); - previous[2] = Floor(position[2]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); + previous[2] = Floor(position[2]); // basis values along the x axis relative = position[0] - (FieldType)previous[0]; @@ -2273,8 +2273,8 @@ void CubicSplineImageGradient2D(const nifti_image *floatingImage, position[0] = world[0] * floatingIJKMatrix->m[0][0] + world[1] * floatingIJKMatrix->m[0][1] + floatingIJKMatrix->m[0][3]; position[1] = world[0] * floatingIJKMatrix->m[1][0] + world[1] * floatingIJKMatrix->m[1][1] + floatingIJKMatrix->m[1][3]; - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); // basis values along the x axis relative = position[0] - (FieldType)previous[0]; relative = relative > 0 ? relative : 0; @@ -2530,7 +2530,7 @@ nifti_image* reg_makeIsotropic(nifti_image *img, int inter) { for (size_t i = 0; i < 8; ++i) newDim[i] = img->dim[i]; for (size_t i = 1; i < 4; ++i) { if (i < static_cast(img->dim[0] + 1)) - newDim[i] = Ceil(img->dim[i] * img->pixdim[i] / smallestPixDim); + newDim[i] = Ceil(img->dim[i] * img->pixdim[i] / smallestPixDim); } // Create the new image nifti_image *newImg = nifti_make_new_nim(newDim, img->datatype, true); diff --git a/reg-lib/cpu/_reg_ssd.cpp b/reg-lib/cpu/_reg_ssd.cpp index 2a4bddfb..50968e3c 100755 --- a/reg-lib/cpu/_reg_ssd.cpp +++ b/reg-lib/cpu/_reg_ssd.cpp @@ -353,9 +353,9 @@ void GetDiscretisedValueSSD_core3D(nifti_image *controlPointGridImage, // Compute the block size const int blockSize[3] = { - Ceil(controlPointGridImage->dx / refImage->dx), - Ceil(controlPointGridImage->dy / refImage->dy), - Ceil(controlPointGridImage->dz / refImage->dz), + Ceil(controlPointGridImage->dx / refImage->dx), + Ceil(controlPointGridImage->dy / refImage->dy), + Ceil(controlPointGridImage->dz / refImage->dz), }; int voxelBlockNumber = blockSize[0] * blockSize[1] * blockSize[2]; int voxelBlockNumber_t = blockSize[0] * blockSize[1] * blockSize[2] * refImage->nt; @@ -405,9 +405,9 @@ void GetDiscretisedValueSSD_core3D(nifti_image *controlPointGridImage, // Compute the corresponding image voxel position Mat44Mul(grid2img_vox, gridVox, imageVox); - imageVox[0] = static_cast(Round(imageVox[0])); - imageVox[1] = static_cast(Round(imageVox[1])); - imageVox[2] = static_cast(Round(imageVox[2])); + imageVox[0] = Round(imageVox[0]); + imageVox[1] = Round(imageVox[1]); + imageVox[2] = Round(imageVox[2]); //INIT for (idBlock = 0; idBlock < voxelBlockNumber_t; idBlock++) { diff --git a/reg-lib/cpu/_reg_tools.cpp b/reg-lib/cpu/_reg_tools.cpp index 21aa5869..d3581339 100755 --- a/reg-lib/cpu/_reg_tools.cpp +++ b/reg-lib/cpu/_reg_tools.cpp @@ -1368,7 +1368,7 @@ void reg_downsampleImage(nifti_image *image, int type, bool *downsampleAxis) { int oldDim[4]; for (int i = 1; i < 4; i++) { oldDim[i] = image->dim[i]; - if (image->dim[i] > 1 && downsampleAxis[i]) image->dim[i] = Ceil(image->dim[i] / 2.0); + if (image->dim[i] > 1 && downsampleAxis[i]) image->dim[i] = Ceil(image->dim[i] / 2.0); if (image->pixdim[i] > 0 && downsampleAxis[i]) image->pixdim[i] = image->pixdim[i] * 2.0f; } image->nx = image->dim[1]; @@ -1446,9 +1446,9 @@ void reg_downsampleImage(nifti_image *image, int type, bool *downsampleAxis) { z * image->qto_xyz.m[2][2] + image->qto_xyz.m[2][3]; // Extract the position in voxel in the old image; - position[0] = Round(real[0] * real2Voxel_qform.m[0][0] + real[1] * real2Voxel_qform.m[0][1] + real[2] * real2Voxel_qform.m[0][2] + real2Voxel_qform.m[0][3]); - position[1] = Round(real[0] * real2Voxel_qform.m[1][0] + real[1] * real2Voxel_qform.m[1][1] + real[2] * real2Voxel_qform.m[1][2] + real2Voxel_qform.m[1][3]); - position[2] = Round(real[0] * real2Voxel_qform.m[2][0] + real[1] * real2Voxel_qform.m[2][1] + real[2] * real2Voxel_qform.m[2][2] + real2Voxel_qform.m[2][3]); + position[0] = Round(real[0] * real2Voxel_qform.m[0][0] + real[1] * real2Voxel_qform.m[0][1] + real[2] * real2Voxel_qform.m[0][2] + real2Voxel_qform.m[0][3]); + position[1] = Round(real[0] * real2Voxel_qform.m[1][0] + real[1] * real2Voxel_qform.m[1][1] + real[2] * real2Voxel_qform.m[1][2] + real2Voxel_qform.m[1][3]); + position[2] = Round(real[0] * real2Voxel_qform.m[2][0] + real[1] * real2Voxel_qform.m[2][1] + real[2] * real2Voxel_qform.m[2][2] + real2Voxel_qform.m[2][3]); if (oldDim[3] == 1) position[2] = 0; // Nearest neighbour is used as downsampling ratio is constant intensity = std::numeric_limits::quiet_NaN(); diff --git a/reg-lib/cuda/CudaLocalTransformation.cu b/reg-lib/cuda/CudaLocalTransformation.cu index 8e901204..10e3ef92 100644 --- a/reg-lib/cuda/CudaLocalTransformation.cu +++ b/reg-lib/cuda/CudaLocalTransformation.cu @@ -275,7 +275,7 @@ void ComputeApproxJacobianValues(const nifti_image *controlPointImage, // The Jacobian matrix is computed for every control point if (controlPointImage->nz > 1) { const unsigned blocks = blockSize->GetApproxJacobianValues3d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)controlPointNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)controlPointNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); GetApproxJacobianValues3d<<>>(jacobianMatricesCuda, jacobianDetCuda, *controlPointTexture, @@ -283,7 +283,7 @@ void ComputeApproxJacobianValues(const nifti_image *controlPointImage, NR_CUDA_CHECK_KERNEL(gridDims, blockDims); } else { const unsigned blocks = blockSize->GetApproxJacobianValues2d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)controlPointNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)controlPointNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); GetApproxJacobianValues2d<<>>(jacobianMatricesCuda, jacobianDetCuda, *controlPointTexture, @@ -311,7 +311,7 @@ void ComputeJacobianValues(const nifti_image *controlPointImage, // The Jacobian matrix is computed for every voxel if (controlPointImage->nz > 1) { const unsigned blocks = blockSize->GetJacobianValues3d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)voxelNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)voxelNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); // 8 floats of shared memory are allocated per thread @@ -322,7 +322,7 @@ void ComputeJacobianValues(const nifti_image *controlPointImage, NR_CUDA_CHECK_KERNEL(gridDims, blockDims); } else { const unsigned blocks = blockSize->GetJacobianValues2d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)voxelNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)voxelNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); GetJacobianValues2d<<>>(jacobianMatricesCuda, jacobianDetCuda, *controlPointTexture, @@ -360,7 +360,7 @@ double GetJacobianPenaltyTerm(const nifti_image *referenceImage, // The Jacobian determinant are squared and logged (might not be english but will do) const unsigned blocks = CudaContext::GetBlockSize()->LogSquaredValues; - const unsigned grids = (unsigned)Ceil(sqrtf((float)jacNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)jacNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); LogSquaredValues<<>>(jacobianDetCuda, (unsigned)jacNumber); @@ -412,7 +412,7 @@ void GetJacobianPenaltyTermGradient(const nifti_image *referenceImage, if (approx) { if (controlPointImage->nz > 1) { const unsigned blocks = blockSize->ComputeApproxJacGradient3d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)controlPointNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)controlPointNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); ComputeApproxJacGradient3d<<>>(transGradientCuda, *jacobianDeterminantTexture, @@ -421,7 +421,7 @@ void GetJacobianPenaltyTermGradient(const nifti_image *referenceImage, NR_CUDA_CHECK_KERNEL(gridDims, blockDims); } else { const unsigned blocks = blockSize->ComputeApproxJacGradient2d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)controlPointNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)controlPointNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); ComputeApproxJacGradient2d<<>>(transGradientCuda, *jacobianDeterminantTexture, @@ -436,7 +436,7 @@ void GetJacobianPenaltyTermGradient(const nifti_image *referenceImage, controlPointImage->dz / referenceImage->dz); if (controlPointImage->nz > 1) { const unsigned blocks = blockSize->ComputeJacGradient3d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)controlPointNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)controlPointNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); ComputeJacGradient3d<<>>(transGradientCuda, *jacobianDeterminantTexture, @@ -446,7 +446,7 @@ void GetJacobianPenaltyTermGradient(const nifti_image *referenceImage, NR_CUDA_CHECK_KERNEL(gridDims, blockDims); } else { const unsigned blocks = blockSize->ComputeJacGradient2d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)controlPointNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)controlPointNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); ComputeJacGradient2d<<>>(transGradientCuda, *jacobianDeterminantTexture, @@ -491,7 +491,7 @@ double CorrectFolding(const nifti_image *referenceImage, NR_CUDA_SAFE_CALL(cudaMalloc(&jacobianDet2Cuda, jacobianDetSize)); NR_CUDA_SAFE_CALL(cudaMemcpy(jacobianDet2Cuda, jacobianDetCuda, jacobianDetSize, cudaMemcpyDeviceToDevice)); const unsigned blocks = blockSize->LogSquaredValues; - const unsigned grids = (unsigned)Ceil(sqrtf((float)jacNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)jacNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); LogSquaredValues<<>>(jacobianDet2Cuda, (unsigned)jacNumber); @@ -520,7 +520,7 @@ double CorrectFolding(const nifti_image *referenceImage, auto jacobianMatricesTexture = Cuda::CreateTextureObject(jacobianMatricesCuda, 9 * jacNumber, cudaChannelFormatKindFloat, 1); if (approx) { const unsigned blocks = blockSize->ApproxCorrectFolding3d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)controlPointNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)controlPointNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); ApproxCorrectFolding3d<<>>(controlPointImageCuda, *jacobianDeterminantTexture, @@ -533,7 +533,7 @@ double CorrectFolding(const nifti_image *referenceImage, controlPointImage->dy / referenceImage->dy, controlPointImage->dz / referenceImage->dz); const unsigned blocks = blockSize->CorrectFolding3d; - const unsigned grids = (unsigned)Ceil(sqrtf((float)controlPointNumber / (float)blocks)); + const unsigned grids = Ceil(sqrtf((float)controlPointNumber / (float)blocks)); const dim3 gridDims(grids, grids, 1); const dim3 blockDims(blocks, 1, 1); CorrectFolding3d<<>>(controlPointImageCuda, *jacobianDeterminantTexture, @@ -691,7 +691,7 @@ void GetDeformationFieldFromFlowField(nifti_image *flowField, squaringNumber = squaringNumber < 6 ? 6 : squaringNumber; // Set the number of squaring step in the flow field if (fabs(flowField->intent_p2) != squaringNumber) - NR_WARN("Changing from " << Round(fabs(flowField->intent_p2)) << " to " << abs(squaringNumber) << + NR_WARN("Changing from " << Round(fabs(flowField->intent_p2)) << " to " << abs(squaringNumber) << " squaring step (equivalent to scaling down by " << (int)pow(2.0f, squaringNumber) << ")"); // Update the number of squaring step required flowField->intent_p2 = static_cast(flowField->intent_p2 >= 0 ? squaringNumber : -squaringNumber); diff --git a/reg-lib/cuda/CudaLocalTransformationKernels.cu b/reg-lib/cuda/CudaLocalTransformationKernels.cu index c3c344be..ba380b20 100644 --- a/reg-lib/cuda/CudaLocalTransformationKernels.cu +++ b/reg-lib/cuda/CudaLocalTransformationKernels.cu @@ -304,7 +304,7 @@ __device__ void GetDeformationField3d(float4 *deformationField, yVoxel < 0 || yVoxel >= referenceImageDims.y || zVoxel < 0 || zVoxel >= referenceImageDims.z) return; - nodePre = { Floor(xVoxel), Floor(yVoxel), Floor(zVoxel) }; + nodePre = { Floor(xVoxel), Floor(yVoxel), Floor(zVoxel) }; basis = { xVoxel - float(nodePre.x--), yVoxel - float(nodePre.y--), zVoxel - float(nodePre.z--) }; } else { // starting deformation field is blank - !composition const auto [x, y, z] = IndexToDims(index, referenceImageDims); @@ -367,7 +367,7 @@ __device__ void GetDeformationField2d(float4 *deformationField, if (xVoxel < 0 || xVoxel >= referenceImageDims.x || yVoxel < 0 || yVoxel >= referenceImageDims.y) return; - nodePre = { Floor(xVoxel), Floor(yVoxel) }; + nodePre = { Floor(xVoxel), Floor(yVoxel) }; basis = { xVoxel - float(nodePre.x--), yVoxel - float(nodePre.y--) }; } else { // starting deformation field is blank - !composition const auto [x, y, z] = IndexToDims(index, referenceImageDims); @@ -556,7 +556,7 @@ __global__ void GetJacobianValues2d(float *jacobianMatrices, const auto [x, y, z] = IndexToDims(tid, referenceImageDims); // the "nearest previous" node is determined [0,0,0] - const int2 nodePre = { Floor((float)x / controlPointSpacing.x), Floor((float)y / controlPointSpacing.y) }; + const int2 nodePre = { Floor((float)x / controlPointSpacing.x), Floor((float)y / controlPointSpacing.y) }; float xBasis[4], yBasis[4], xFirst[4], yFirst[4], relative; @@ -624,9 +624,9 @@ __global__ void GetJacobianValues3d(float *jacobianMatrices, // the "nearest previous" node is determined [0,0,0] const int3 nodePre = { - Floor((float)x / controlPointSpacing.x), - Floor((float)y / controlPointSpacing.y), - Floor((float)z / controlPointSpacing.z) + Floor((float)x / controlPointSpacing.x), + Floor((float)y / controlPointSpacing.y), + Floor((float)z / controlPointSpacing.z) }; extern __shared__ float yFirst[]; @@ -872,14 +872,14 @@ __global__ void ComputeJacGradient2d(float4 *gradient, const auto [x, y, z] = IndexToDims(tid, controlPointImageDims); float2 jacobianGradient{}; - for (int pixelY = Ceil((y - 3) * controlPointVoxelSpacing.y); pixelY <= Ceil((y + 1) * controlPointVoxelSpacing.y); ++pixelY) { + for (int pixelY = Ceil((y - 3) * controlPointVoxelSpacing.y); pixelY <= Ceil((y + 1) * controlPointVoxelSpacing.y); ++pixelY) { if (-1 < pixelY && pixelY < referenceImageDims.y) { const int yPre = (int)((float)pixelY / controlPointVoxelSpacing.y); float basis = (float)pixelY / controlPointVoxelSpacing.y - (float)yPre; float yBasis, yFirst; GetBSplineBasisValue(basis, y - yPre, &yBasis, &yFirst); - for (int pixelX = Ceil((x - 3) * controlPointVoxelSpacing.x); pixelX <= Ceil((x + 1) * controlPointVoxelSpacing.x); ++pixelX) { + for (int pixelX = Ceil((x - 3) * controlPointVoxelSpacing.x); pixelX <= Ceil((x + 1) * controlPointVoxelSpacing.x); ++pixelX) { if (-1 < pixelX && pixelX < referenceImageDims.x && (yFirst != 0.f || yBasis != 0.f)) { const int xPre = (int)((float)pixelX / controlPointVoxelSpacing.x); basis = (float)pixelX / controlPointVoxelSpacing.x - (float)xPre; @@ -925,21 +925,21 @@ __global__ void ComputeJacGradient3d(float4 *gradient, const auto [x, y, z] = IndexToDims(tid, controlPointImageDims); float3 jacobianGradient{}; - for (int pixelZ = Ceil((z - 3) * controlPointVoxelSpacing.z); pixelZ <= Ceil((z + 1) * controlPointVoxelSpacing.z); ++pixelZ) { + for (int pixelZ = Ceil((z - 3) * controlPointVoxelSpacing.z); pixelZ <= Ceil((z + 1) * controlPointVoxelSpacing.z); ++pixelZ) { if (-1 < pixelZ && pixelZ < referenceImageDims.z) { const int zPre = (int)((float)pixelZ / controlPointVoxelSpacing.z); float basis = (float)pixelZ / controlPointVoxelSpacing.z - (float)zPre; float zBasis, zFirst; GetBSplineBasisValue(basis, z - zPre, &zBasis, &zFirst); - for (int pixelY = Ceil((y - 3) * controlPointVoxelSpacing.y); pixelY <= Ceil((y + 1) * controlPointVoxelSpacing.y); ++pixelY) { + for (int pixelY = Ceil((y - 3) * controlPointVoxelSpacing.y); pixelY <= Ceil((y + 1) * controlPointVoxelSpacing.y); ++pixelY) { if (-1 < pixelY && pixelY < referenceImageDims.y && (zFirst != 0.f || zBasis != 0.f)) { const int yPre = (int)((float)pixelY / controlPointVoxelSpacing.y); basis = (float)pixelY / controlPointVoxelSpacing.y - (float)yPre; float yBasis, yFirst; GetBSplineBasisValue(basis, y - yPre, &yBasis, &yFirst); - for (int pixelX = Ceil((x - 3) * controlPointVoxelSpacing.x); pixelX <= Ceil((x + 1) * controlPointVoxelSpacing.x); ++pixelX) { + for (int pixelX = Ceil((x - 3) * controlPointVoxelSpacing.x); pixelX <= Ceil((x + 1) * controlPointVoxelSpacing.x); ++pixelX) { if (-1 < pixelX && pixelX < referenceImageDims.x && (yFirst != 0.f || yBasis != 0.f)) { const int xPre = (int)((float)pixelX / controlPointVoxelSpacing.x); basis = (float)pixelX / controlPointVoxelSpacing.x - (float)xPre; @@ -1063,11 +1063,11 @@ __global__ void CorrectFolding3d(float4 *controlPointGrid, const auto [x, y, z] = IndexToDims(tid, controlPointImageDims); float3 foldingCorrection{}; - for (int pixelZ = Ceil((z - 3) * controlPointVoxelSpacing.z); pixelZ < Ceil((z + 1) * controlPointVoxelSpacing.z); ++pixelZ) { + for (int pixelZ = Ceil((z - 3) * controlPointVoxelSpacing.z); pixelZ < Ceil((z + 1) * controlPointVoxelSpacing.z); ++pixelZ) { if (-1 < pixelZ && pixelZ < referenceImageDims.z) { - for (int pixelY = Ceil((y - 3) * controlPointVoxelSpacing.y); pixelY < Ceil((y + 1) * controlPointVoxelSpacing.y); ++pixelY) { + for (int pixelY = Ceil((y - 3) * controlPointVoxelSpacing.y); pixelY < Ceil((y + 1) * controlPointVoxelSpacing.y); ++pixelY) { if (-1 < pixelY && pixelY < referenceImageDims.y) { - for (int pixelX = Ceil((x - 3) * controlPointVoxelSpacing.x); pixelX < Ceil((x + 1) * controlPointVoxelSpacing.x); ++pixelX) { + for (int pixelX = Ceil((x - 3) * controlPointVoxelSpacing.x); pixelX < Ceil((x + 1) * controlPointVoxelSpacing.x); ++pixelX) { if (-1 < pixelX && pixelX < referenceImageDims.x) { int jacIndex = (pixelZ * referenceImageDims.y + pixelY) * referenceImageDims.x + pixelX; float detJac = tex1Dfetch(jacobianDeterminantTexture, jacIndex); @@ -1141,7 +1141,7 @@ __device__ void DefFieldComposeKernel(float4 *deformationField, }; // Linear interpolation - const int3 ante = { Floor(voxelPosition.x), Floor(voxelPosition.y), Floor(voxelPosition.z) }; + const int3 ante = { Floor(voxelPosition.x), Floor(voxelPosition.y), Floor(voxelPosition.z) }; float relX[2], relY[2], relZ[2]; relX[1] = voxelPosition.x - (float)ante.x; relX[0] = 1.f - relX[1]; relY[1] = voxelPosition.y - (float)ante.y; relY[0] = 1.f - relY[1]; @@ -1173,7 +1173,7 @@ __device__ void DefFieldComposeKernel(float4 *deformationField, }; // Linear interpolation - const int2 ante = { Floor(voxelPosition.x), Floor(voxelPosition.y) }; + const int2 ante = { Floor(voxelPosition.x), Floor(voxelPosition.y) }; float relX[2], relY[2]; relX[1] = voxelPosition.x - (float)ante.x; relX[0] = 1.f - relX[1]; relY[1] = voxelPosition.y - (float)ante.y; relY[0] = 1.f - relY[1]; diff --git a/reg-lib/cuda/CudaResampling.cu b/reg-lib/cuda/CudaResampling.cu index 58c33998..37722ed1 100644 --- a/reg-lib/cuda/CudaResampling.cu +++ b/reg-lib/cuda/CudaResampling.cu @@ -49,12 +49,12 @@ __inline__ __device__ void TransformInterpolate(const mat44 matrix, const float4 } // Compute the linear interpolation - previous.x = Floor(voxelDeformation[0]); - previous.y = Floor(voxelDeformation[1]); + previous.x = Floor(voxelDeformation[0]); + previous.y = Floor(voxelDeformation[1]); InterpLinearKernel(voxelDeformation[0] - static_cast(previous.x), xBasis); InterpLinearKernel(voxelDeformation[1] - static_cast(previous.y), yBasis); if constexpr (is3d) { - previous.z = Floor(voxelDeformation[2]); + previous.z = Floor(voxelDeformation[2]); InterpLinearKernel(voxelDeformation[2] - static_cast(previous.z), zBasis); } } diff --git a/reg-lib/cuda/CudaToolsKernels.cu b/reg-lib/cuda/CudaToolsKernels.cu index f502ac4f..c0bba323 100644 --- a/reg-lib/cuda/CudaToolsKernels.cu +++ b/reg-lib/cuda/CudaToolsKernels.cu @@ -30,7 +30,7 @@ __device__ void VoxelCentricToNodeCentricKernel(float4 *nodeImageCuda, // Linear interpolation float basisX[2], basisY[2], basisZ[2], interpolatedValue[3]{}; - const int pre[3] = { Floor(voxelCoord[0]), Floor(voxelCoord[1]), Floor(voxelCoord[2]) }; + const int pre[3] = { Floor(voxelCoord[0]), Floor(voxelCoord[1]), Floor(voxelCoord[2]) }; basisX[1] = voxelCoord[0] - static_cast(pre[0]); basisX[0] = 1.f - basisX[1]; basisY[1] = voxelCoord[1] - static_cast(pre[1]); diff --git a/reg-lib/cuda/resampleKernel.cu b/reg-lib/cuda/resampleKernel.cu index 50bcb91c..b4df65fc 100644 --- a/reg-lib/cuda/resampleKernel.cu +++ b/reg-lib/cuda/resampleKernel.cu @@ -214,8 +214,8 @@ __global__ void ResampleImage2D(float* floatingImage, // real -> voxel; floating space reg_mat44_mul_cuda(sourceIJKMatrix, world, position); - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); relative[0] = (double)(position[0]) - (double)(previous[0]); relative[1] = (double)(position[1]) - (double)(previous[1]); @@ -307,9 +307,9 @@ __global__ void ResampleImage3D(float* floatingImage, // real -> voxel; floating space reg_mat44_mul_cuda(sourceIJKMatrix, world, position); - previous[0] = Floor(position[0]); - previous[1] = Floor(position[1]); - previous[2] = Floor(position[2]); + previous[0] = Floor(position[0]); + previous[1] = Floor(position[1]); + previous[2] = Floor(position[2]); relative[0] = (double)(position[0]) - (double)(previous[0]); relative[1] = (double)(position[1]) - (double)(previous[1]); diff --git a/reg-test/reg_test_voxelCentricToNodeCentric.cpp b/reg-test/reg_test_voxelCentricToNodeCentric.cpp index b588989a..8a5d099f 100644 --- a/reg-test/reg_test_voxelCentricToNodeCentric.cpp +++ b/reg-test/reg_test_voxelCentricToNodeCentric.cpp @@ -192,7 +192,7 @@ class VoxelCentricToNodeCentricTest { Mat44Mul(transformation, nodeCoord, voxelCoord); // Linear interpolation DataType basisX[2], basisY[2], basisZ[2]; - const int pre[3] = { Floor(voxelCoord[0]), Floor(voxelCoord[1]), Floor(voxelCoord[2]) }; + const int pre[3] = { Floor(voxelCoord[0]), Floor(voxelCoord[1]), Floor(voxelCoord[2]) }; basisX[1] = voxelCoord[0] - static_cast(pre[0]); basisX[0] = static_cast(1) - basisX[1]; basisY[1] = voxelCoord[1] - static_cast(pre[1]);