diff --git a/README.txt b/README.txt index f9ad9667..264a4073 100644 --- a/README.txt +++ b/README.txt @@ -15,7 +15,8 @@ presented by Ourselin et al.[1]. The symmetric versions of the rigid and affine registration have been presented in Modat et al.[2]. The non-linear registration is based on the work is based on the work initially presented by Rueckert et al.[3]. The current implementation has been presented -in Modat et al.[4]. +in Modat et al.[4]. Sliding-region handling was added as it was presented +at WBIR 2018 in Eiben et al.[5]. Ourselin et al.[1] presented an algorithm called Aladin, which is based on a block-matching approach and a Trimmed Least Square (TLS) scheme. Firstly, @@ -43,7 +44,10 @@ an objective function composed from the Normalised Mutual Information (NMI) and the Bending-Energy (BE) is used. The objective function value is optimised using the analytical derivative of both, the NMI and the BE within a conjugate gradient scheme. The symmetric version of the algorithm takes advantage of -stationary velocity field parametrisation. +stationary velocity field parametrisation. Optional sliding region handling was +implemented using two separate control-point grids, one for each region [5]. +The boundary between the regions is defined as the zero-crossing of a signed +distance map that needs to be given in the floating image space. reg f3d is the command to perform non-linear registration. A third program, called reg resample, is been embedded in the package. It @@ -118,6 +122,11 @@ Imaging, 18(8), 712–721. doi:10.1109/42.796284 [4] Modat, et al. (2010). Fast free-form deformation using graphics processing units. Computer Methods And Programs In Biomedicine,98(3), 278–284. doi:10.1016/j.cmpb.2009.09.002 +[5] Eiben et al. (2018). +Eiben, et al. (2018) Statistical Motion Mask and Sliding Registration. +Biomedical Image Registration, WBIR, 13-23 +doi:10.1007/978-3-319-92258-4_2 + ############################################################################## ############################################################################## diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 38b10c1b..39f5b693 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -68 +71 diff --git a/reg-apps/reg_f3d.cpp b/reg-apps/reg_f3d.cpp index 76135cef..7ea3bd47 100755 --- a/reg-apps/reg_f3d.cpp +++ b/reg-apps/reg_f3d.cpp @@ -13,6 +13,7 @@ #include "_reg_ReadWriteMatrix.h" #include "_reg_f3d2.h" #include "reg_f3d.h" +#include "_reg_f3d_sli.h" #include //#include //DOES NOT WORK ON WINDOWS ! @@ -49,7 +50,7 @@ void Usage(char *exec) reg_print_info(exec, "***************"); reg_print_info(exec, "*** Initial transformation options (One option will be considered):"); reg_print_info(exec, "\t-aff \t\tFilename which contains an affine transformation (Affine*Reference=Floating)"); - reg_print_info(exec, "\t-incpp \tFilename ofloatf control point grid input"); + reg_print_info(exec, "\t-incpp \tFilename of control point grid input"); reg_print_info(exec, "\t\t\t\tThe coarse spacing is defined by this file."); reg_print_info(exec, ""); reg_print_info(exec, "*** Output options:"); @@ -131,6 +132,15 @@ void Usage(char *exec) reg_print_info(exec, "\t-fmask \tFilename of a mask image in the floating space"); reg_print_info(exec, ""); + //sliding regions registration options + reg_print_info(exec, "*** Sliding regions options:"); + reg_print_info(exec, "\t-sli \t\t\tUse two transformations and a penalty term to allow for sliding regions."); + reg_print_info(exec, "\t\t\t\tFor more details see Eiben et al., Statistical motion mask and Sliding image registration, WBIR, 2018" ); + reg_print_info(exec, "\t-dmap \tFilename of a distance map image in the floating space, defining the two sliding regions"); + reg_print_info(exec, "\t-go \t\tWeight of the gap-overlap penalty term [0.1]"); + reg_print_info(exec, "\t-incpp2 \tFilename of control point grid input for region 2"); + reg_print_info(exec, ""); + // reg_print_info(exec, "*** Platform options:"); //#if defined(_USE_CUDA) && defined(_USE_OPENCL) // reg_print_info(exec, "\t-platf \t\tChoose platform: CPU=0 | Cuda=1 | OpenCL=2 [0]"); @@ -297,8 +307,7 @@ int main(int argc, char **argv) //\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ // Check the type of registration object to create reg_f3d *REG=NULL; - float *referenceLandmark=NULL; - float *floatingLandmark=NULL; + reg_f3d_sli *REG_SLI = NULL; for(int i=1; i(referenceImage->nt,floatingImage->nt); break; } + if (strcmp(argv[i], "-sli") == 0 || strcmp(argv[i], "--sli") == 0) + { + REG_SLI = new reg_f3d_sli(referenceImage->nt, floatingImage->nt); + REG = static_cast *>(REG_SLI); + break; + } } if(REG==NULL) REG=new reg_f3d(referenceImage->nt,floatingImage->nt); @@ -328,6 +343,10 @@ int main(int argc, char **argv) bool useMeanLNCC=false; int refBinNumber=0; int floBinNumber=0; + float *referenceLandmark = NULL; + float *floatingLandmark = NULL; + nifti_image *distMapImage = NULL; + nifti_image *inputCPPImageR2 = NULL; /* read the input parameter */ for(int i=1; iUseBCHUpdate(atoi(argv[++i])); } + //sliding region registration options + else if (strcmp(argv[i], "-dmap") == 0 || strcmp(argv[i], "--dmap") == 0) + { + distMapImage = reg_io_ReadImageFile(argv[++i]); + if (distMapImage == NULL) + { + reg_print_msg_error("Error when reading the distance map image"); + reg_print_msg_error(argv[i - 1]); + return EXIT_FAILURE; + } + if (REG_SLI != NULL) + { + REG_SLI->SetDistanceMapImage(distMapImage); + } + else + { + reg_print_msg_error("Sliding regions registrations (-sli flag) must be used when specifying a distance map image"); + return EXIT_FAILURE; + } + } + else if (strcmp(argv[i], "-go") == 0 || strcmp(argv[i], "--go") == 0) + { + if (REG_SLI != NULL) + { + REG_SLI->SetGapOverlapWeight(atof(argv[++i])); + } + else + { + reg_print_msg_error("Sliding regions registrations (-sli flag) must be used when specifying the gap-overlap penalty term weight"); + return EXIT_FAILURE; + } + } + else if (strcmp(argv[i], "-incpp2") == 0 || strcmp(argv[i], "--incpp2") == 0) + { + inputCPPImageR2 = reg_io_ReadImageFile(argv[++i]); + if (inputCPPImageR2 == NULL) + { + reg_print_msg_error("Error when reading the input control point grid image for region 2:"); + reg_print_msg_error(argv[i - 1]); + return EXIT_FAILURE; + } + if (REG_SLI != NULL) + { + REG_SLI->SetRegion2ControlPointGridImage(inputCPPImageR2); + } + else + { + reg_print_msg_error("Sliding regions registrations (-sli flag) must be used when specifying a input control point grid for region 2"); + return EXIT_FAILURE; + } + } + else if(strcmp(argv[i], "-omp")==0 || strcmp(argv[i], "--omp")==0) { #if defined (_OPENMP) @@ -802,7 +873,8 @@ int main(int argc, char **argv) strcmp(argv[i], "-Version")!=0 && strcmp(argv[i], "-V")!=0 && strcmp(argv[i], "-v")!=0 && strcmp(argv[i], "--v")!=0 && strcmp(argv[i], "-gpu")!=0 && strcmp(argv[i], "--gpu")!=0 && - strcmp(argv[i], "-vel")!=0 && strcmp(argv[i], "-sym")!=0) + strcmp(argv[i], "-vel")!=0 && strcmp(argv[i], "-sym")!=0 && + strcmp(argv[i], "-sli") != 0) { reg_print_msg_error("\tParameter unknown:"); reg_print_msg_error(argv[i]); @@ -836,44 +908,119 @@ int main(int argc, char **argv) REG->Run(); // Save the control point image - nifti_image *outputControlPointGridImage = REG->GetControlPointPositionImage(); if(outputCPPImageName==NULL) outputCPPImageName=(char *)"outputCPP.nii"; - memset(outputControlPointGridImage->descrip, 0, 80); - strcpy (outputControlPointGridImage->descrip,"Control point position from NiftyReg (reg_f3d)"); - if(strcmp("NiftyReg F3D2", REG->GetExecutableName())==0) - strcpy (outputControlPointGridImage->descrip,"Velocity field grid from NiftyReg (reg_f3d2)"); - reg_io_WriteImageFile(outputControlPointGridImage,outputCPPImageName); - nifti_image_free(outputControlPointGridImage); - outputControlPointGridImage=NULL; + //check if sliding region registration + if (REG_SLI == NULL) + { + //if not get output image, set description, save, and clear + nifti_image *outputControlPointGridImage = REG->GetControlPointPositionImage(); + memset(outputControlPointGridImage->descrip, 0, 80); + strcpy(outputControlPointGridImage->descrip, "Control point position from NiftyReg (reg_f3d)"); + if (strcmp("NiftyReg F3D2", REG->GetExecutableName()) == 0) + strcpy(outputControlPointGridImage->descrip, "Velocity field grid from NiftyReg (reg_f3d2)"); + reg_io_WriteImageFile(outputControlPointGridImage, outputCPPImageName); + nifti_image_free(outputControlPointGridImage); + outputControlPointGridImage = NULL; - // Save the backward control point image - if(REG->GetSymmetricStatus()) + // and save backwards CPP image if exists + if (REG->GetSymmetricStatus()) + { + // _backward is added to the forward control point grid image name + std::string b(outputCPPImageName); + if (b.find(".nii.gz") != std::string::npos) + b.replace(b.find(".nii.gz"), 7, "_backward.nii.gz"); + else if (b.find(".nii") != std::string::npos) + b.replace(b.find(".nii"), 4, "_backward.nii"); + else if (b.find(".hdr") != std::string::npos) + b.replace(b.find(".hdr"), 4, "_backward.hdr"); + else if (b.find(".img.gz") != std::string::npos) + b.replace(b.find(".img.gz"), 7, "_backward.img.gz"); + else if (b.find(".img") != std::string::npos) + b.replace(b.find(".img"), 4, "_backward.img"); + else if (b.find(".png") != std::string::npos) + b.replace(b.find(".png"), 4, "_backward.png"); + else if (b.find(".nrrd") != std::string::npos) + b.replace(b.find(".nrrd"), 5, "_backward.nrrd"); + else b.append("_backward.nii"); + nifti_image *outputBackwardControlPointGridImage = REG->GetBackwardControlPointPositionImage(); + memset(outputBackwardControlPointGridImage->descrip, 0, 80); + strcpy(outputBackwardControlPointGridImage->descrip, "Backward Control point position from NiftyReg (reg_f3d)"); + if (strcmp("NiftyReg F3D2", REG->GetExecutableName()) == 0) + strcpy(outputBackwardControlPointGridImage->descrip, "Backward velocity field grid from NiftyReg (reg_f3d2)"); + reg_io_WriteImageFile(outputBackwardControlPointGridImage, b.c_str()); + nifti_image_free(outputBackwardControlPointGridImage); + outputBackwardControlPointGridImage = NULL; + } + } + else { - // _backward is added to the forward control point grid image name - std::string b(outputCPPImageName); - if(b.find( ".nii.gz") != std::string::npos) - b.replace(b.find( ".nii.gz"),7,"_backward.nii.gz"); - else if(b.find( ".nii") != std::string::npos) - b.replace(b.find( ".nii"),4,"_backward.nii"); - else if(b.find( ".hdr") != std::string::npos) - b.replace(b.find( ".hdr"),4,"_backward.hdr"); - else if(b.find( ".img.gz") != std::string::npos) - b.replace(b.find( ".img.gz"),7,"_backward.img.gz"); - else if(b.find( ".img") != std::string::npos) - b.replace(b.find( ".img"),4,"_backward.img"); - else if(b.find( ".png") != std::string::npos) - b.replace(b.find( ".png"),4,"_backward.png"); - else if(b.find( ".nrrd") != std::string::npos) - b.replace(b.find( ".nrrd"),5,"_backward.nrrd"); - else b.append("_backward.nii"); - nifti_image *outputBackwardControlPointGridImage = REG->GetBackwardControlPointPositionImage(); - memset(outputBackwardControlPointGridImage->descrip, 0, 80); - strcpy (outputBackwardControlPointGridImage->descrip,"Backward Control point position from NiftyReg (reg_f3d)"); - if(strcmp("NiftyReg F3D2", REG->GetExecutableName())==0) - strcpy (outputBackwardControlPointGridImage->descrip,"Backward velocity field grid from NiftyReg (reg_f3d2)"); - reg_io_WriteImageFile(outputBackwardControlPointGridImage,b.c_str()); - nifti_image_free(outputBackwardControlPointGridImage); - outputBackwardControlPointGridImage=NULL; + //if sliding registration create two versions of CPP name, + //one with _region1 appended and the other with _region2 appended + std::string region1Name(outputCPPImageName); + std::string region2Name(outputCPPImageName); + if (region1Name.find(".nii.gz") != std::string::npos) + { + region1Name.replace(region1Name.find(".nii.gz"), 7, "_region1.nii.gz"); + region2Name.replace(region2Name.find(".nii.gz"), 7, "_region2.nii.gz"); + } + else if (region1Name.find(".nii") != std::string::npos) + { + region1Name.replace(region1Name.find(".nii"), 4, "_region1.nii"); + region2Name.replace(region2Name.find(".nii"), 4, "_region2.nii"); + } + else if (region1Name.find(".hdr") != std::string::npos) + { + region1Name.replace(region1Name.find(".hdr"), 4, "_region1.hdr"); + region2Name.replace(region2Name.find(".hdr"), 4, "_region2.hdr"); + } + else if (region1Name.find(".img.gz") != std::string::npos) + { + region1Name.replace(region1Name.find(".img.gz"), 7, "_region1.img.gz"); + region2Name.replace(region2Name.find(".img.gz"), 7, "_region2.img.gz"); + } + else if (region1Name.find(".img") != std::string::npos) + { + region1Name.replace(region1Name.find(".img"), 4, "_region1.img"); + region2Name.replace(region2Name.find(".img"), 4, "_region2.img"); + } + else if (region1Name.find(".png") != std::string::npos) + { + region1Name.replace(region1Name.find(".png"), 4, "_region1.png"); + region2Name.replace(region2Name.find(".png"), 4, "_region2.png"); + } + else if (region1Name.find(".nrrd") != std::string::npos) + { + region1Name.replace(region1Name.find(".nrrd"), 5, "_region1.nrrd"); + region2Name.replace(region2Name.find(".nrrd"), 5, "_region2.nrrd"); + } + else + { + region1Name.append("_region1.nii"); + region2Name.append("_region2.nii"); + } + + //now get output CPP image for region 1 + nifti_image *region1OutputCPPImage = REG_SLI->GetControlPointPositionImage(); + //set description for CPP image for region 1 + memset(region1OutputCPPImage->descrip, 0, 80); + strcpy(region1OutputCPPImage->descrip, "Control point position for Region 1 from NiftyReg (reg_f3d_sli)"); + //save CPP image for region 1 + reg_io_WriteImageFile(region1OutputCPPImage, region1Name.c_str()); + //clear CPP image for region 1 + nifti_image_free(region1OutputCPPImage); + region1OutputCPPImage = NULL; + + //and same for region 2... + //get output CPP image for region 2 + nifti_image *region2OutputCPPImage = REG_SLI->GetRegion2ControlPointPositionImage(); + //set description for CPP image for region 2 + memset(region2OutputCPPImage->descrip, 0, 80); + strcpy(region2OutputCPPImage->descrip, "Control point position for Region 2 from NiftyReg (reg_f3d_sli)"); + //save CPP image for region 2 + reg_io_WriteImageFile(region2OutputCPPImage, region2Name.c_str()); + //clear CPP image for region 2 + nifti_image_free(region2OutputCPPImage); + region2OutputCPPImage = NULL; } // Save the warped image(s) @@ -936,6 +1083,8 @@ int main(int argc, char **argv) if(inputCCPImage!=NULL) nifti_image_free(inputCCPImage); if(referenceMaskImage!=NULL) nifti_image_free(referenceMaskImage); if(floatingMaskImage!=NULL) nifti_image_free(floatingMaskImage); + if (distMapImage != NULL) nifti_image_free(distMapImage); + if (inputCPPImageR2 != NULL) nifti_image_free(inputCPPImageR2); #ifdef NDEBUG if(verbose) diff --git a/reg-lib/CMakeLists.txt b/reg-lib/CMakeLists.txt index 9c14fce6..0889be86 100755 --- a/reg-lib/CMakeLists.txt +++ b/reg-lib/CMakeLists.txt @@ -204,6 +204,8 @@ set(_reg_f3d_files _reg_f3d2.cpp _reg_f3d_sym.h _reg_f3d_sym.cpp + _reg_f3d_sli.h + _reg_f3d_sli.cpp ) set(_reg_f3d_libraries _reg_localTrans @@ -224,6 +226,7 @@ install(FILES _reg_base.h DESTINATION include) install(FILES _reg_f3d.h DESTINATION include) install(FILES _reg_f3d2.h DESTINATION include) install(FILES _reg_f3d_sym.h DESTINATION include) +install(FILES _reg_f3d_sli.h DESTINATION include) install(FILES cpu/_reg_optimiser.cpp cpu/_reg_optimiser.h DESTINATION include) set(NIFTYREG_LIBRARIES "${NIFTYREG_LIBRARIES};_reg_f3d") #----------------------------------------------------------------------------- diff --git a/reg-lib/_reg_base.cpp b/reg-lib/_reg_base.cpp index 11258600..40f51e60 100644 --- a/reg-lib/_reg_base.cpp +++ b/reg-lib/_reg_base.cpp @@ -1532,11 +1532,6 @@ void reg_base::Run() this->currentMask = this->maskPyramid[0]; } - // Allocate image that depends on the reference image - this->AllocateWarped(); - this->AllocateDeformationField(); - this->AllocateWarpedGradient(); - // The grid is refined if necessary T maxStepSize=this->InitialiseCurrentLevel(); T currentSize = maxStepSize; @@ -1544,7 +1539,12 @@ void reg_base::Run() this->DisplayCurrentLevelParameters(); - // Allocate image that are required to compute the gradient + // Allocate image that depends on the reference image + this->AllocateWarped(); + this->AllocateDeformationField(); + this->AllocateWarpedGradient(); + + // Allocate image that are required to compute the gradient this->AllocateVoxelBasedMeasureGradient(); this->AllocateTransformationGradient(); diff --git a/reg-lib/_reg_f3d.cpp b/reg-lib/_reg_f3d.cpp index eb36f032..29447c1e 100644 --- a/reg-lib/_reg_f3d.cpp +++ b/reg-lib/_reg_f3d.cpp @@ -543,6 +543,13 @@ void reg_f3d::GetSimilarityMeasureGradient() { this->GetVoxelBasedGradient(); + + /*char im_name[100]; + sprintf(im_name, "debug_vox_grad_lev_%i_iter_%i.nii", this->currentLevel, this->optimiser->GetCurrentIterationNumber()); + reg_io_WriteImageFile(this->voxelBasedMeasureGradient, im_name);*/ + + + int kernel_type=CUBIC_SPLINE_KERNEL; // The voxel based NMI gradient is convolved with a spline kernel // Convolution along the x axis @@ -593,6 +600,12 @@ void reg_f3d::GetSimilarityMeasureGradient() false, // no update &reorientation ); + + + /*sprintf(im_name, "debug_cp_grad_lev_%i_iter_%i.nii", this->currentLevel, this->optimiser->GetCurrentIterationNumber()); + reg_io_WriteImageFile(this->transformationGradient, im_name); + reg_exit(0);*/ + #ifndef NDEBUG reg_print_fct_debug("reg_f3d::GetSimilarityMeasureGradient"); #endif @@ -1023,10 +1036,10 @@ nifti_image **reg_f3d::GetWarpedImage() this->currentFloating = this->inputFloating; this->currentMask=NULL; - reg_base::AllocateWarped(); - reg_base::AllocateDeformationField(); - reg_base::WarpFloatingImage(3); // cubic spline interpolation - reg_base::ClearDeformationField(); + this->AllocateWarped(); + this->AllocateDeformationField(); + this->WarpFloatingImage(3); // cubic spline interpolation + this->ClearDeformationField(); nifti_image **warpedImage= (nifti_image **)malloc(2*sizeof(nifti_image *)); warpedImage[0]=nifti_copy_nim_info(this->warped); diff --git a/reg-lib/_reg_f3d_sli.cpp b/reg-lib/_reg_f3d_sli.cpp new file mode 100644 index 00000000..470f2e56 --- /dev/null +++ b/reg-lib/_reg_f3d_sli.cpp @@ -0,0 +1,1694 @@ +/* +* _reg_f3d_sli.cpp +* +* +* Created by Jamie McClelland on 20/08/2017. +* Copyright (c) 2017, University College London. All rights reserved. +* Centre for Medical Image Computing (CMIC) +* See the LICENSE.txt file in the nifty_reg root folder +* +*/ + + +#ifndef _REG_F3D_SLI_CPP +#define _REG_F3D_SLI_CPP + +#include "_reg_f3d_sli.h" + +/* *************************************************************** */ +/* *************************************************************** */ +template +reg_f3d_sli::reg_f3d_sli(int refTimePoint, int floTimePoint) + :reg_f3d::reg_f3d(refTimePoint, floTimePoint) +{ + this->executableName = (char *)"NiftyReg F3D Sliding regions"; + + this->inputRegion2ControlPointGrid = NULL; + this->region2ControlPointGrid = NULL; + this->region2DeformationFieldImage = NULL; + this->region2VoxelBasedMeasureGradientImage = NULL; + this->region2TransformationGradient = NULL; + + this->region1DeformationFieldImage = NULL; + this->region1VoxelBasedMeasureGradientImage = NULL; + + this->inputDistanceMap = NULL; + this->distanceMapPyramid = NULL; + this->currentDistanceMap = NULL; + + this->warpedDistanceMapRegion1 = NULL; + this->warpedDistanceMapRegion2 = NULL; + this->warpedDistanceMapGradientRegion1 = NULL; + this->warpedDistanceMapGradientRegion2 = NULL; + + this->gapOverlapGradientWRTDefFieldRegion1 = NULL; + this->gapOverlapGradientWRTDefFieldRegion2 = NULL; + + this->gapOverlapWeight = 0.1; + this->currentWGO = 0; + this->bestWGO = 0; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::reg_f3d_sli"); +#endif +} +/* *************************************************************** */ +template +reg_f3d_sli::~reg_f3d_sli() +{ + if (this->region2ControlPointGrid != NULL) + { + nifti_image_free(this->region2ControlPointGrid); + this->region2ControlPointGrid = NULL; + } + + if (this->distanceMapPyramid != NULL) + { + if (this->usePyramid) + { + for (unsigned int n = 0; n < this->levelToPerform; n++) + { + if (this->distanceMapPyramid[n] != NULL) + { + nifti_image_free(this->distanceMapPyramid[n]); + this->distanceMapPyramid[n] = NULL; + } + } + } + else + { + if (this->distanceMapPyramid[0] != NULL) + { + nifti_image_free(this->distanceMapPyramid[0]); + this->distanceMapPyramid[0] = NULL; + } + } + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::~reg_f3d_sli"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +template +void reg_f3d_sli::GetDeformationField() +{ + + //get deformation field for region1 + //do not use a mask as gap-overlap penalty term calculted for voxels outside mask + reg_spline_getDeformationField(this->controlPointGrid, + this->region1DeformationFieldImage, + NULL, //no mask + false, //composition + true); //bspline + //get deformation field for region2 + reg_spline_getDeformationField(this->region2ControlPointGrid, + this->region2DeformationFieldImage, + NULL, //no mask + false, //composition + true); //bspline + + //warp distance map using region1 def field + reg_resampleImage(this->currentDistanceMap, + this->warpedDistanceMapRegion1, + this->region1DeformationFieldImage, + NULL, //no mask + this->interpolation, + std::numeric_limits::quiet_NaN()); //set padding value to NaN + //warp distance map using region2 def field + reg_resampleImage(this->currentDistanceMap, + this->warpedDistanceMapRegion2, + this->region2DeformationFieldImage, + NULL, //no mask + this->interpolation, + std::numeric_limits::quiet_NaN()); //set padding value to NaN + + //loop over voxels and set combined deformation field (deformationFieldImage) + //using appropriate region, based on warped distance maps + //combined def field only needs to be set within the mask + size_t numVox = this->region1DeformationFieldImage->nx * + this->region1DeformationFieldImage->ny * + this->region1DeformationFieldImage->nz; + //pointers to deformation fields + T *region1DFPtrX = static_cast(this->region1DeformationFieldImage->data); + T *region1DFPtrY = ®ion1DFPtrX[numVox]; + T *region1DFPtrZ = NULL; + T *region2DFPtrX = static_cast(this->region2DeformationFieldImage->data); + T *region2DFPtrY = ®ion2DFPtrX[numVox]; + T *region2DFPtrZ = NULL; + T *combinedDFPtrX = static_cast(this->deformationFieldImage->data); + T *combinedDFPtrY = &combinedDFPtrX[numVox]; + T *combinedDFPtrZ = NULL; + //pointers to warped distance maps + T *warpedDMR1Ptr = static_cast(this->warpedDistanceMapRegion1->data); + T *warpedDMR2Ptr = static_cast(this->warpedDistanceMapRegion2->data); + //are images 3D? + if (this->region1DeformationFieldImage->nz > 1) + { + //extra pointers required for 3D + region1DFPtrZ = ®ion1DFPtrY[numVox]; + region2DFPtrZ = ®ion2DFPtrY[numVox]; + combinedDFPtrZ = &combinedDFPtrY[numVox]; + } + + //loop over voxels + for (size_t n = 0; n < numVox; n++) + { + //check mask exists and in mask + if (this->currentMask == NULL || this->currentMask[n] > -1) + { + //warped distance maps (WDMs) will contain NaN values if the transform + //maps the voxel outside the extent of the distance map so need to check + //for NaN values + if (warpedDMR1Ptr[n] != warpedDMR1Ptr[n]) + { + if (warpedDMR2Ptr[n] != warpedDMR2Ptr[n]) + { + //both WDMs are NaN so set combined def field to NaN + combinedDFPtrX[n] = std::numeric_limits::quiet_NaN(); + combinedDFPtrY[n] = std::numeric_limits::quiet_NaN(); + if (combinedDFPtrZ != NULL) + combinedDFPtrZ[n] = std::numeric_limits::quiet_NaN(); + } + else + { + //check if region2 transform maps into region1, i.e. if region2 WDM < 0 + if (warpedDMR2Ptr[n] < 0) + { + //set combined def field to NaN + combinedDFPtrX[n] = std::numeric_limits::quiet_NaN(); + combinedDFPtrY[n] = std::numeric_limits::quiet_NaN(); + if (combinedDFPtrZ != NULL) + combinedDFPtrZ[n] = std::numeric_limits::quiet_NaN(); + } + else + { + //set combined def field to region2 def field + combinedDFPtrX[n] = region2DFPtrX[n]; + combinedDFPtrY[n] = region2DFPtrY[n]; + if (combinedDFPtrZ != NULL) + combinedDFPtrZ[n] = region2DFPtrZ[n]; + } + } + }//if (warpedDMR1Ptr[n] != warpedDMR1Ptr[n]) + else + { + //region1 WDM is not NaN, but still need to check region2 WDM + if (warpedDMR2Ptr[n] != warpedDMR2Ptr[n]) + { + //region2 WDM is NaN so check if region1 transform maps into region2, i.e. if region1 WDM >= 0 + if (warpedDMR1Ptr[n] >= 0) + { + //set combined def field to NaN + combinedDFPtrX[n] = std::numeric_limits::quiet_NaN(); + combinedDFPtrY[n] = std::numeric_limits::quiet_NaN(); + if (combinedDFPtrZ != NULL) + combinedDFPtrZ[n] = std::numeric_limits::quiet_NaN(); + } + else + { + //set combined def field to region1 def field + combinedDFPtrX[n] = region1DFPtrX[n]; + combinedDFPtrY[n] = region1DFPtrY[n]; + if (combinedDFPtrZ != NULL) + combinedDFPtrZ[n] = region1DFPtrZ[n]; + } + } + else + { + //region1 WDM and region2 WDM are both not NaN + //so if sum of WDMs < 0 set combined def field to region1 def field + //if >= 0 set combined def field to region2 def field + if ((warpedDMR1Ptr[n] + warpedDMR2Ptr[n]) < 0) + { + combinedDFPtrX[n] = region1DFPtrX[n]; + combinedDFPtrY[n] = region1DFPtrY[n]; + if (combinedDFPtrZ != NULL) + combinedDFPtrZ[n] = region1DFPtrZ[n]; + } + else + { + combinedDFPtrX[n] = region2DFPtrX[n]; + combinedDFPtrY[n] = region2DFPtrY[n]; + if (combinedDFPtrZ != NULL) + combinedDFPtrZ[n] = region2DFPtrZ[n]; + } + }//else (warpedDMR2Ptr[n] != warpedDMR2Ptr[n]) + }//else (warpedDMR1Ptr[n] != warpedDMR1Ptr[n]) + }//if (this->currentMask[n] > -1) + //not in mask so set combined def field to NaN + }//for (size_t n = 0; n < numVox; n++) + + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetDeformationField()"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::AllocateDeformationField() +{ + //clear any previously allocated deformation fields + this->ClearDeformationField(); + + //call method from reg_base to allocate combined deformation field + reg_base::AllocateDeformationField(); + + //now allocate def fields for regions 1 and 2 using header info from combined def field + this->region1DeformationFieldImage = nifti_copy_nim_info(this->deformationFieldImage); + this->region1DeformationFieldImage->data = (void *)calloc(this->region1DeformationFieldImage->nvox, + this->region1DeformationFieldImage->nbyper); + this->region2DeformationFieldImage = nifti_copy_nim_info(this->deformationFieldImage); + this->region2DeformationFieldImage->data = (void *)calloc(this->region2DeformationFieldImage->nvox, + this->region2DeformationFieldImage->nbyper); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::AllocateDeformationField"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::ClearDeformationField() +{ + //call method from reg_base to clear combined def field + reg_base::ClearDeformationField(); + + //now clear def fields for regions 1 and 2 + if (this->region1DeformationFieldImage != NULL) + { + nifti_image_free(this->region1DeformationFieldImage); + this->region1DeformationFieldImage = NULL; + } + if (this->region2DeformationFieldImage != NULL) + { + nifti_image_free(this->region2DeformationFieldImage); + this->region2DeformationFieldImage = NULL; + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ClearDeformationField"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::AllocateWarped() +{ + //clear any previously allocated warped images + this->ClearWarped(); + + //call method from reg_base to allocate warped floating image + reg_base::AllocateWarped(); + + //Allocate warped distance maps for region 1 and region 2 + //use header info from warped image, but update some info using header from current distance map + this->warpedDistanceMapRegion1 = nifti_copy_nim_info(this->warped); + this->warpedDistanceMapRegion2 = nifti_copy_nim_info(this->warped); + this->warpedDistanceMapRegion1->dim[0] = this->warpedDistanceMapRegion1->ndim = + this->warpedDistanceMapRegion2->dim[0] = this->warpedDistanceMapRegion2->ndim = + this->currentDistanceMap->ndim; + this->warpedDistanceMapRegion1->dim[4] = this->warpedDistanceMapRegion1->nt = + this->warpedDistanceMapRegion2->dim[4] = this->warpedDistanceMapRegion2->nt = + this->currentDistanceMap->nt; + this->warpedDistanceMapRegion1->nvox = this->warpedDistanceMapRegion2->nvox = + this->warpedDistanceMapRegion1->nx * + this->warpedDistanceMapRegion1->ny * + this->warpedDistanceMapRegion1->nz * + this->warpedDistanceMapRegion1->nt; + this->warpedDistanceMapRegion1->datatype = this->warpedDistanceMapRegion2->datatype = this->currentDistanceMap->datatype; + this->warpedDistanceMapRegion1->nbyper = this->warpedDistanceMapRegion2->nbyper = this->currentDistanceMap->nbyper; + //now allocate memory for warped distance maps data + this->warpedDistanceMapRegion1->data = (void *)calloc(this->warpedDistanceMapRegion1->nvox, this->warpedDistanceMapRegion1->nbyper); + this->warpedDistanceMapRegion2->data = (void *)calloc(this->warpedDistanceMapRegion2->nvox, this->warpedDistanceMapRegion2->nbyper); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::AllocateWarped"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::ClearWarped() +{ + //call method from reg_base to clear warped floating image + reg_base::ClearWarped(); + + //now clear warped distance maps + if (this->warpedDistanceMapRegion1 != NULL) + { + nifti_image_free(this->warpedDistanceMapRegion1); + this->warpedDistanceMapRegion1 = NULL; + } + if (this->warpedDistanceMapRegion2 != NULL) + { + nifti_image_free(this->warpedDistanceMapRegion2); + this->warpedDistanceMapRegion2 = NULL; + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ClearWarped"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +template +double reg_f3d_sli::GetObjectiveFunctionValue() +{ + //call method from reg_f3d to calculate objective function value for similarity + //measure and other penalty terms + double objFuncValue = reg_f3d::GetObjectiveFunctionValue(); + + //calculate weighted gap-overlap penalty term + this->currentWGO = this->ComputeGapOverlapPenaltyTerm(); + +#ifndef NDEBUG + char text[255]; + sprintf(text, " | (wGO) %g", this->currentWGO); + reg_print_msg_debug(text); + reg_print_fct_debug("reg_f3d::GetObjectiveFunctionValue"); +#endif + + //return objective function value with weighted gap-overlap value subtracted + return objFuncValue - this->currentWGO; +} +/* *************************************************************** */ +template +double reg_f3d_sli::ComputeGapOverlapPenaltyTerm() +{ + //NOTE - this method assumes the current warped distance maps (WDMs) have already + //been calculated by calling the GetDeformationField() method prior to calling this + //method. The GetDeformtionField method will usually be called when warping the image + //to calculate the image similarities, so this prevents re-calculating the WDMs + //unnecessarily, but if the image similarities all have a weight of 0 and therefore + //the warped image is not calculated, the GetDeformationField() method must still be + //called. + + //NOTE2 - the gap-overlap penalty term is calculated at all voxels within the reference + //image, even if they are outside the mask or have a NaN value in the reference or + //warped image - this is to ensure the transformations for the 2 regions are free of + //gaps and overlaps, even in areas where the images are not being used to drive the + //registration + + if (this->gapOverlapWeight <= 0) + return 0.; + + //loop over all voxels and sum up gap-overlap penalty term values from each voxel. + //the gap-overlap penalty term is defined as -WDM1*WDM2 if WDM1*WDM2<0 (i.e. the + //WDMs point to different regions, indicating a gap or overlap), and 0 otherwise + double gapOverlapTotal = 0.; + double gapOverlapValue = 0.; + + //pointers to warped distance maps + T *warpedDMR1Ptr = static_cast(this->warpedDistanceMapRegion1->data); + T *warpedDMR2Ptr = static_cast(this->warpedDistanceMapRegion2->data); + + size_t numVox = this->warpedDistanceMapRegion1->nx * + this->warpedDistanceMapRegion1->ny * + this->warpedDistanceMapRegion1->nz; + for (size_t n = 0; n < numVox; n++) + { + gapOverlapValue = warpedDMR1Ptr[n] * warpedDMR2Ptr[n]; + //if NaN value in either WDM then gapOverlapValue = NaN, so will fail + //test for less than 0 + if (gapOverlapValue < 0) + gapOverlapTotal -= gapOverlapValue; + } + + //normalise by the number of voxels and return weighted value + //gapOverlapTotal /= double(numVox); + return double(this->gapOverlapWeight) * gapOverlapTotal; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ComputeGapOverlapPenaltyTerm()"); +#endif +} +/* *************************************************************** */ +template +double reg_f3d_sli ::ComputeBendingEnergyPenaltyTerm() +{ + //check if penalty term used, i.e. weight > 0 + if (this->bendingEnergyWeight <= 0) return 0.; + + //calculate the bending energy penalty term for region 1 + double region1PenaltyTerm = reg_f3d::ComputeBendingEnergyPenaltyTerm(); + + //calculate the bending energy penalty term for region 2 + double region2PenaltyTerm = this->bendingEnergyWeight * reg_spline_approxBendingEnergy(this->region2ControlPointGrid); +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ComputeBendingEnergyPenaltyTerm"); +#endif + return region1PenaltyTerm + region2PenaltyTerm; +} +/* *************************************************************** */ +template +double reg_f3d_sli::ComputeLinearEnergyPenaltyTerm() +{ + //check if penalty term used, i.e. weight > 0 + if (this->linearEnergyWeight <= 0) return 0.; + + //calculate the linear energy penalty term for region 1 + double region1PenaltyTerm = reg_f3d::ComputeLinearEnergyPenaltyTerm(); + + //calculate the bending energy penalty term for region 2 + double region2PenaltyTerm = this->linearEnergyWeight * reg_spline_approxLinearEnergy(this->region2ControlPointGrid); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ComputeLinearEnergyPenaltyTerm"); +#endif + return region1PenaltyTerm + region2PenaltyTerm; +} +/* *************************************************************** */ +template +double reg_f3d_sli::ComputeJacobianBasedPenaltyTerm(int type) +{ + //check if penalty term used, i.e. weight > 0 + if (this->jacobianLogWeight <= 0) return 0.; + + //Jacobian penalty term not currently implemented for sliding region registrations + //so throw error + reg_print_fct_error("reg_f3d_sli::ComputeJacobianBasedPenaltyTerm"); + reg_print_msg_error("Jacobian penalty term not currently implemented for sliding region registrations"); + reg_exit(); +} +/* *************************************************************** */ +template +double reg_f3d_sli::ComputeLandmarkDistancePenaltyTerm() +{ + //check if penalty term used, i.e. weight > 0 + if (this->landmarkRegWeight <= 0) return 0.; + + //Landmark penalty term not currently implemented for sliding region registrations + //so throw error + reg_print_fct_error("reg_f3d_sli::ComputeLandmarkDistancePenaltyTerm"); + reg_print_msg_error("Landmark distance penalty term not currently implemented for sliding region registrations"); + reg_exit(); +} +/* *************************************************************** */ +/* *************************************************************** */ +template +void reg_f3d_sli::GetObjectiveFunctionGradient() +{ + //note - cannot call method from reg_f3d as objective function gradient will + //be smoothed before the gap-overlap gradient is added to it, so need to + //reproduce code here + + //check if gradient is approximated + if (!this->useApproxGradient) + { + // Compute the gradient of the similarity measure + if (this->similarityWeight>0) + { + this->WarpFloatingImage(this->interpolation); + this->GetSimilarityMeasureGradient(); + } + else + { + this->SetGradientImageToZero(); + } + // Compute the penalty term gradients if required + this->GetBendingEnergyGradient(); + this->GetJacobianBasedGradient(); + this->GetLinearEnergyGradient(); + this->GetLandmarkDistanceGradient(); + //include the gap-penalty term gradient + this->GetGapOverlapGradient(); + } + else + { + this->GetApproximatedGradient(); + } + + //increment the optimiser iteration number + this->optimiser->IncrementCurrentIterationNumber(); + + // Smooth the gradient if required + this->SmoothGradient(); +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetObjectiveFunctionGradient"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::GetSimilarityMeasureGradient() +{ + //get voxel-based similairty gradient + this->GetVoxelBasedGradient(); + + //The voxel based gradient images for each region are filled with zeros + reg_tools_multiplyValueToImage(this->region1VoxelBasedMeasureGradientImage, + this->region1VoxelBasedMeasureGradientImage, + 0.f); + reg_tools_multiplyValueToImage(this->region2VoxelBasedMeasureGradientImage, + this->region2VoxelBasedMeasureGradientImage, + 0.f); + + //pointers to warped distance maps + T *warpedDMR1Ptr = static_cast(this->warpedDistanceMapRegion1->data); + T *warpedDMR2Ptr = static_cast(this->warpedDistanceMapRegion2->data); + //pointers to voxel-based similarity gradients + size_t numVox = this->voxelBasedMeasureGradient->nx * + this->voxelBasedMeasureGradient->ny * + this->voxelBasedMeasureGradient->nz; + T *combinedVBMGPtrX = static_cast(this->voxelBasedMeasureGradient->data); + T *combinedVBMGPtrY = &combinedVBMGPtrX[numVox]; + T *combinedVBMGPtrZ = NULL; + T *region1VBMGPtrX = static_cast(this->region1VoxelBasedMeasureGradientImage->data); + T *region1VBMGPtrY = ®ion1VBMGPtrX[numVox]; + T *region1VBMGPtrZ = NULL; + T *region2VBMGPtrX = static_cast(this->region2VoxelBasedMeasureGradientImage->data); + T *region2VBMGPtrY = ®ion2VBMGPtrX[numVox]; + T *region2VBMGPtrZ = NULL; + //check for 3D + if (this->voxelBasedMeasureGradient->nz > 1) + { + combinedVBMGPtrZ = &combinedVBMGPtrY[numVox]; + region1VBMGPtrZ = ®ion1VBMGPtrY[numVox]; + region2VBMGPtrZ = ®ion2VBMGPtrY[numVox]; + } + //loop over voxels and split voxel-based gradient between two regions + //based on warped distance maps (WDMs). + //Note - GetDeformationField() will be called before this method, so + //WDMs will have already been calculated + for (size_t n = 0; n < numVox; n++) + { + //is in mask? + if (this->currentMask[n] > -1) + { + //need to check for NaNs in WDMs + //if WDM1 is NaN and WDM2 >= 0 (indicating region2 transform maps into region 2) + //then copy voxel-based gradient value in to region2VoxelBasedMeasureGradientImage + if (warpedDMR1Ptr[n] != warpedDMR1Ptr[n] && warpedDMR2Ptr[n] >= 0) + { + region2VBMGPtrX[n] = combinedVBMGPtrX[n]; + region2VBMGPtrY[n] = combinedVBMGPtrY[n]; + if (combinedVBMGPtrZ != NULL) + { + region2VBMGPtrZ[n] = combinedVBMGPtrZ[n]; + } + } + //if WDM2 is NaN and WDM1 < 0 (indicating region1 transform maps into region 1) + //then copy voxel-based gradient value in to region1VoxelBasedMeasureGradientImage + if (warpedDMR2Ptr[n] != warpedDMR2Ptr[n] && warpedDMR1Ptr[n] < 0) + { + region1VBMGPtrX[n] = combinedVBMGPtrX[n]; + region1VBMGPtrY[n] = combinedVBMGPtrY[n]; + if (combinedVBMGPtrZ != NULL) + { + region1VBMGPtrZ[n] = combinedVBMGPtrZ[n]; + } + } + //if both WDMs are not NaN then assign voxel-based gradient value to correct region + //based on WDMs + if (warpedDMR1Ptr[n] == warpedDMR1Ptr[n] && warpedDMR2Ptr[n] == warpedDMR2Ptr[n]) + { + //if sum of WDMs < 0 assign value to region 1, else assign to region 2 + if (warpedDMR1Ptr[n] + warpedDMR2Ptr[n] < 0) + { + region1VBMGPtrX[n] = combinedVBMGPtrX[n]; + region1VBMGPtrY[n] = combinedVBMGPtrY[n]; + if (combinedVBMGPtrZ != NULL) + { + region1VBMGPtrZ[n] = combinedVBMGPtrZ[n]; + } + } + else + { + region2VBMGPtrX[n] = combinedVBMGPtrX[n]; + region2VBMGPtrY[n] = combinedVBMGPtrY[n]; + if (combinedVBMGPtrZ != NULL) + { + region2VBMGPtrZ[n] = combinedVBMGPtrZ[n]; + } + } + } + } + } + + + /*char im_name[100]; + sprintf(im_name, "debug_vox_grad_comb_lev_%i_iter_%i.nii", this->currentLevel, this->optimiser->GetCurrentIterationNumber()); + reg_io_WriteImageFile(this->voxelBasedMeasureGradient, im_name); + sprintf(im_name, "debug_vox_grad_R1_lev_%i_iter_%i.nii", this->currentLevel, this->optimiser->GetCurrentIterationNumber()); + reg_io_WriteImageFile(this->region1VoxelBasedMeasureGradientImage, im_name); + sprintf(im_name, "debug_vox_grad_R2_lev_%i_iter_%i.nii", this->currentLevel, this->optimiser->GetCurrentIterationNumber()); + reg_io_WriteImageFile(this->region2VoxelBasedMeasureGradientImage, im_name);*/ + + + + //convert voxel-based gradienta to CPG gradients for both regions + + //first convolve voxel-based gardient with a spline kernel + int kernel_type = CUBIC_SPLINE_KERNEL; + // Convolution along the x axis + float currentNodeSpacing[3]; + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->controlPointGrid->dx; + bool activeAxis[3] = { 1, 0, 0 }; + reg_tools_kernelConvolution(this->region1VoxelBasedMeasureGradientImage, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + reg_tools_kernelConvolution(this->region2VoxelBasedMeasureGradientImage, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + // Convolution along the y axis + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->controlPointGrid->dy; + activeAxis[0] = 0; + activeAxis[1] = 1; + reg_tools_kernelConvolution(this->region1VoxelBasedMeasureGradientImage, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + reg_tools_kernelConvolution(this->region2VoxelBasedMeasureGradientImage, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + // Convolution along the z axis if required + if (this->region1VoxelBasedMeasureGradientImage->nz > 1) + { + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->controlPointGrid->dz; + activeAxis[1] = 0; + activeAxis[2] = 1; + reg_tools_kernelConvolution(this->region1VoxelBasedMeasureGradientImage, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + reg_tools_kernelConvolution(this->region2VoxelBasedMeasureGradientImage, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + } + //now resample voxel-based gradients at control points to get transformationGradients + //the gradients need to be reorientated to account for the transformation from distance + //map image coordinates to world coordinates + mat44 reorientation; + if (this->currentFloating->sform_code>0) + reorientation = this->currentFloating->sto_ijk; + else reorientation = this->currentFloating->qto_ijk; + reg_voxelCentric2NodeCentric(this->transformationGradient, + this->region1VoxelBasedMeasureGradientImage, + this->similarityWeight, + false, // no update + &reorientation); + reg_voxelCentric2NodeCentric(this->region2TransformationGradient, + this->region2VoxelBasedMeasureGradientImage, + this->similarityWeight, + false, // no update + &reorientation); + + + /*sprintf(im_name, "debug_cp_grad_R1_lev_%i_iter_%i.nii", this->currentLevel, this->optimiser->GetCurrentIterationNumber()); + reg_io_WriteImageFile(this->transformationGradient, im_name); + sprintf(im_name, "debug_cp_grad_R2_lev_%i_iter_%i.nii", this->currentLevel, this->optimiser->GetCurrentIterationNumber()); + reg_io_WriteImageFile(this->region2TransformationGradient, im_name); + reg_exit(0);*/ + + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetSimilarityMeasureGradient()"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +template +void reg_f3d_sli::GetGapOverlapGradient() +{ + //NOTE - this method assumes the deformation fields and the WDMs for each region + //have already been calculated by calling the GetDeformationField() method prior + //to calling this method. + + //first calculate gap-overlap gradient with respect to def field for each region + //then convolve these with a B-spline kernal to get the gap-overlap gradient with + //respect to the transform (i.e. the CPG) for each region + // + //the gap-overlap gradient with respect to the def field for region 1 is: + //dGO/dDF1 = -WDM2*(dWDM1/dDF1) if WDM1*WDM2<0 else 0 + //where dWMD1/dDF1 is the spatial gradient of the distance map warped by the def + //field for region 1 + // + //the gap-overlap gradient with respect to the def field for region 2 is: + //dGO/dDF2 = -WDM1*(dWDM2/dDF2) if WDM1*WDM2<0 else 0 + //where dWMD2/dDF2 is the spatial gradient of the distance map warped by the def + //field for region 2 + + //The gap-overlap gradients WRT the def fields for each region are filled with zeros + reg_tools_multiplyValueToImage(this->gapOverlapGradientWRTDefFieldRegion1, + this->gapOverlapGradientWRTDefFieldRegion1, + 0.f); + reg_tools_multiplyValueToImage(this->gapOverlapGradientWRTDefFieldRegion2, + this->gapOverlapGradientWRTDefFieldRegion2, + 0.f); + + //calculate the spatial gradient of the distance map warped by the def field from + //each region + reg_getImageGradient(this->currentDistanceMap, + this->warpedDistanceMapGradientRegion1, + this->region1DeformationFieldImage, + this->currentMask, + this->interpolation, + this->warpedPaddingValue, + 0);//timepoint 0 + reg_getImageGradient(this->currentDistanceMap, + this->warpedDistanceMapGradientRegion2, + this->region2DeformationFieldImage, + this->currentMask, + this->interpolation, + this->warpedPaddingValue, + 0);//timepoint 0 + + //pointers to warped distance maps + T *warpedDMR1Ptr = static_cast(this->warpedDistanceMapRegion1->data); + T *warpedDMR2Ptr = static_cast(this->warpedDistanceMapRegion2->data); + //pointers to warped spatial gradients + size_t numVox = this->warpedDistanceMapRegion1->nx * + this->warpedDistanceMapRegion1->ny * + this->warpedDistanceMapRegion1->nz; + T *warpedDMGradR1PtrX = static_cast(this->warpedDistanceMapGradientRegion1->data); + T *warpedDMGradR1PtrY = &warpedDMGradR1PtrX[numVox]; + T *warpedDMGradR1PtrZ = NULL; + T *warpedDMGradR2PtrX = static_cast(this->warpedDistanceMapGradientRegion2->data); + T *warpedDMGradR2PtrY = &warpedDMGradR2PtrX[numVox]; + T *warpedDMGradR2PtrZ = NULL; + //pointers to the gap-overlap gradients WRT def field for each region + T *gapOverlapGradR1PtrX = static_cast(this->gapOverlapGradientWRTDefFieldRegion1->data); + T *gapOverlapGradR1PtrY = &gapOverlapGradR1PtrX[numVox]; + T *gapOverlapGradR1PtrZ = NULL; + T *gapOverlapGradR2PtrX = static_cast(this->gapOverlapGradientWRTDefFieldRegion2->data); + T *gapOverlapGradR2PtrY = &gapOverlapGradR2PtrX[numVox]; + T *gapOverlapGradR2PtrZ = NULL; + //check for 3D + if (this->warpedDistanceMapGradientRegion1->nz > 1) + { + warpedDMGradR1PtrZ = &warpedDMGradR1PtrY[numVox]; + warpedDMGradR2PtrZ = &warpedDMGradR2PtrY[numVox]; + gapOverlapGradR1PtrZ = &gapOverlapGradR1PtrY[numVox]; + gapOverlapGradR2PtrZ = &gapOverlapGradR2PtrY[numVox]; + } + + //loop over all voxels and calculate gap-overlap gradient with respect to def field + //for each region + for (size_t n = 0; n < numVox; n++) + { + if (warpedDMR1Ptr[n] * warpedDMR2Ptr[n] < 0) + { + //dGO / dDF1 = -WDM2*(dWDM1 / dDF1) + gapOverlapGradR1PtrX[n] = -warpedDMR2Ptr[n] * warpedDMGradR1PtrX[n]; + gapOverlapGradR1PtrY[n] = -warpedDMR2Ptr[n] * warpedDMGradR1PtrY[n]; + //dGO / dDF2 = -WDM1*(dWDM2 / dDF2) + gapOverlapGradR2PtrX[n] = -warpedDMR1Ptr[n] * warpedDMGradR2PtrX[n]; + gapOverlapGradR2PtrY[n] = -warpedDMR1Ptr[n] * warpedDMGradR2PtrY[n]; + //check for 3D + if (gapOverlapGradR1PtrZ != NULL) + { + gapOverlapGradR1PtrZ[n] = -warpedDMR2Ptr[n] * warpedDMGradR1PtrZ[n]; + gapOverlapGradR2PtrZ[n] = -warpedDMR1Ptr[n] * warpedDMGradR2PtrZ[n]; + } + }//if (warpedDMR1Ptr[n] * warpedDMR2Ptr[n]) + }//for (size_t n = 0; n < numVox; n++) + + //the gap-overlap gradient WRT the def field is convolved with a B-spline kernel + //to calculate the gradient WRT the CPG for each region + int kernel_type = CUBIC_SPLINE_KERNEL; + // Convolution along the x axis + float currentNodeSpacing[3]; + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->controlPointGrid->dx; + bool activeAxis[3] = { 1, 0, 0 }; + reg_tools_kernelConvolution(this->gapOverlapGradientWRTDefFieldRegion1, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + reg_tools_kernelConvolution(this->gapOverlapGradientWRTDefFieldRegion2, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + // Convolution along the y axis + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->controlPointGrid->dy; + activeAxis[0] = 0; + activeAxis[1] = 1; + reg_tools_kernelConvolution(this->gapOverlapGradientWRTDefFieldRegion1, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + reg_tools_kernelConvolution(this->gapOverlapGradientWRTDefFieldRegion2, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + // Convolution along the z axis if required + if (this->gapOverlapGradientWRTDefFieldRegion1->nz > 1) + { + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->controlPointGrid->dz; + activeAxis[1] = 0; + activeAxis[2] = 1; + reg_tools_kernelConvolution(this->gapOverlapGradientWRTDefFieldRegion1, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + reg_tools_kernelConvolution(this->gapOverlapGradientWRTDefFieldRegion2, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); + } + + //the voxel-wise gradients are now resampled at the CPG locations and added to the + //transformation gradients for each region + //the gradients need to be reorientated to account for the transformation from distance + //map image coordinates to world coordinates + mat44 reorientation; + if (this->currentDistanceMap->sform_code>0) + reorientation = this->currentDistanceMap->sto_ijk; + else reorientation = this->currentDistanceMap->qto_ijk; + reg_voxelCentric2NodeCentric(this->transformationGradient, + this->gapOverlapGradientWRTDefFieldRegion1, + this->gapOverlapWeight, + true, // update the transformation gradient + &reorientation); + reg_voxelCentric2NodeCentric(this->region2TransformationGradient, + this->gapOverlapGradientWRTDefFieldRegion2, + this->gapOverlapWeight, + true, // update the transformation gradient + &reorientation); + +} +/* *************************************************************** */ +template +void reg_f3d_sli::GetBendingEnergyGradient() +{ + //check if bending energy used + if (this->bendingEnergyWeight <= 0) return; + + //calculate bending energy gradient for region 1 transform + reg_f3d::GetBendingEnergyGradient(); + + //calculate bending energy gradient for region 2 transform + reg_spline_approxBendingEnergyGradient(this->region2ControlPointGrid, + this->region2TransformationGradient, + this->bendingEnergyWeight); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetBendingEnergyGradient"); +#endif + return; +} +/* *************************************************************** */ +template +void reg_f3d_sli::GetLinearEnergyGradient() +{ + //check if linear energy used + if (this->linearEnergyWeight <= 0) return; + + //calculate linear energy gradient for region 1 transform + reg_f3d::GetLinearEnergyGradient(); + + //calculate linear energy gradient for region 2 transform + reg_spline_approxLinearEnergyGradient(this->region2ControlPointGrid, + this->region2TransformationGradient, + this->linearEnergyWeight); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetLinearEnergyGradient"); +#endif + return; +} +/* *************************************************************** */ +template +void reg_f3d_sli::GetJacobianBasedGradient() +{ + //check if penalty term used, i.e. weight > 0 + if (this->jacobianLogWeight <= 0) return; + + //Jacobian penalty term not currently implemented for sliding region registrations + //so throw error + reg_print_fct_error("reg_f3d_sli::GetJacobianBasedGradient"); + reg_print_msg_error("Jacobian penalty term not currently implemented for sliding region registrations"); + reg_exit(); +} +/* *************************************************************** */ +template +void reg_f3d_sli::GetLandmarkDistanceGradient() +{ + //check if penalty term used, i.e. weight > 0 + if (this->landmarkRegWeight <= 0) return; + + //Landmark penalty term not currently implemented for sliding region registrations + //so throw error + reg_print_fct_error("reg_f3d_sli::GetLandmarkDistanceGradient"); + reg_print_msg_error("Landmark distance penalty term not currently implemented for sliding region registrations"); + reg_exit(); +} +/* *************************************************************** */ +template +void reg_f3d_sli::SetGradientImageToZero() +{ + //call method from reg_f3d to set region 1 gradient image to 0 + reg_f3d::SetGradientImageToZero(); + + //set region 2 gradient image to 0 + T* nodeGradPtr = static_cast(this->region2TransformationGradient->data); + for (size_t i = 0; iregion2TransformationGradient->nvox; ++i) + *nodeGradPtr++ = 0; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::SetGradientImageToZero"); +#endif +} +/* *************************************************************** */ +template +T reg_f3d_sli::NormaliseGradient() +{ + // call method from reg_f3d to calculate max length of region 1 gradient image + // note - this method does not normalise the gradient (as the executable name + // is not "NiftyReg F3D"), it will just return the max length + T region1MaxValue = reg_f3d::NormaliseGradient(); + + // The max length of the region 2 gradient image is calculated + T maxGradValue = 0; + size_t voxNumber = this->region2TransformationGradient->nx * + this->region2TransformationGradient->ny * + this->region2TransformationGradient->nz; + //pointers to gradient data + T *r2PtrX = static_cast(this->region2TransformationGradient->data); + T *r2PtrY = &r2PtrX[voxNumber]; + T *r2PtrZ = NULL; + //check for 3D + if (this->region2TransformationGradient->nz > 1) + r2PtrZ = &r2PtrY[voxNumber]; + //loop over voxels, calculate length of gradient vector (ignoring dimension(s) not + //being optimised), and store value if greater than current max value + for (size_t i = 0; ioptimiseX == true) + valX = *r2PtrX++; + if (this->optimiseY == true) + valY = *r2PtrY++; + if (r2PtrZ != NULL && this->optimiseZ == true) + valZ = *r2PtrZ++; + T length = (T)(sqrt(valX*valX + valY*valY + valZ*valZ)); + maxGradValue = (length > maxGradValue) ? length : maxGradValue; + } + + // The largest value between the region 1 and region 2 gradients is kept + maxGradValue = maxGradValue>region1MaxValue ? maxGradValue : region1MaxValue; +#ifndef NDEBUG + char text[255]; + sprintf(text, "Objective function gradient maximal length: %g", maxGradValue); + reg_print_msg_debug(text); +#endif + + // The region 1 gradient is normalised + T *r1Ptr = static_cast(this->transformationGradient->data); + for (size_t i = 0; i < this->transformationGradient->nvox; ++i) + { + *r1Ptr++ /= maxGradValue; + } + // The region 2 gradient is normalised + T *r2Ptr = static_cast(this->region2TransformationGradient->data); + for (size_t i = 0; iregion2TransformationGradient->nvox; ++i) + { + *r2Ptr++ /= maxGradValue; + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::NormaliseGradient"); +#endif + // Returns the largest gradient distance + return maxGradValue; +} +/* *************************************************************** */ +template +void reg_f3d_sli::SmoothGradient() +{ + //check if gradients require smoothing + if (this->gradientSmoothingSigma != 0) + { + //call method from reg_f3d to smooth gradient for region 1 transform + reg_f3d::SmoothGradient(); + + //smooth the gradient for region 2 transform + float kernel = fabs(this->gradientSmoothingSigma); + reg_tools_kernelConvolution(this->region2TransformationGradient, + &kernel, + GAUSSIAN_KERNEL); + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::SmoothGradient"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::GetApproximatedGradient() +{ + //call method from reg_f3d to approximate gradient for region 1 + reg_f3d::GetApproximatedGradient(); + + // approximate gradient for region 2 using finite differences + // + //pointers to region 2 CPG and gradient + T *r2CPGPtr = static_cast(this->region2ControlPointGrid->data); + T *r2GradPtr = static_cast(this->region2TransformationGradient->data); + //amount to increase/decrease CPG values by + //equal to floating voxel size in x dimension / 1000 + T eps = this->currentFloating->dx / 1000.f; + //loop over CPG values + for (size_t i = 0; iregion2ControlPointGrid->nvox; i++) + { + //get best CPG value from optimiser + T currentValue = this->optimiser->GetBestDOF_b()[i]; + //increase the value by eps and calculate new objective function value + r2CPGPtr[i] = currentValue + eps; + double valPlus = this->GetObjectiveFunctionValue(); + //decrease the value by eps and calculate new objective function value + r2CPGPtr[i] = currentValue - eps; + double valMinus = this->GetObjectiveFunctionValue(); + //reset CPG to best value + r2CPGPtr[i] = currentValue; + //set the value of gradient by approximating using finite differences + r2GradPtr[i] = -(T)((valPlus - valMinus) / (2.0*eps)); + } +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetApproximatedGradient"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::AllocateWarpedGradient() +{ + //clear any previously allocated warped gradient images + this->ClearWarpedGradient(); + + //call method from reg_base to allocate warped (floating) gradient image + reg_base::AllocateWarpedGradient(); + + //allocate warped distance map gradient images using header info from + //warped (floating) gradient image + this->warpedDistanceMapGradientRegion1 = nifti_copy_nim_info(this->warImgGradient); + this->warpedDistanceMapGradientRegion2 = nifti_copy_nim_info(this->warImgGradient); + this->warpedDistanceMapGradientRegion1->data = (void *)calloc(this->warpedDistanceMapGradientRegion1->nvox, + this->warpedDistanceMapGradientRegion1->nbyper); + this->warpedDistanceMapGradientRegion2->data = (void *)calloc(this->warpedDistanceMapGradientRegion2->nvox, + this->warpedDistanceMapGradientRegion2->nbyper); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::AllocateWarpedGradient"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::ClearWarpedGradient() +{ + //call method from reg_base to clear warped (floating) gradient image + reg_base::ClearWarpedGradient(); + + //now clear warped distance map gradient images + if (this->warpedDistanceMapGradientRegion1 != NULL) + { + nifti_image_free(this->warpedDistanceMapGradientRegion1); + this->warpedDistanceMapGradientRegion1 = NULL; + } + if (this->warpedDistanceMapGradientRegion2 != NULL) + { + nifti_image_free(this->warpedDistanceMapGradientRegion2); + this->warpedDistanceMapGradientRegion2 = NULL; + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ClearWarpedGradient"); +#endif + +} +/* *************************************************************** */ +/* *************************************************************** */ +template +void reg_f3d_sli::AllocateVoxelBasedMeasureGradient() +{ + //clear any previously allocated images + this->ClearVoxelBasedMeasureGradient(); + + //call method from reg_base to allocate voxel-based similarity measure gradient image + //for combined transform + reg_base::AllocateVoxelBasedMeasureGradient(); + + //allocate voxel-based similarity measure gradient images for each region + this->region1VoxelBasedMeasureGradientImage = nifti_copy_nim_info(this->voxelBasedMeasureGradient); + this->region2VoxelBasedMeasureGradientImage = nifti_copy_nim_info(this->voxelBasedMeasureGradient); + this->region1VoxelBasedMeasureGradientImage->data = (void *)calloc(this->region1VoxelBasedMeasureGradientImage->nvox, + this->region1VoxelBasedMeasureGradientImage->nbyper); + this->region2VoxelBasedMeasureGradientImage->data = (void *)calloc(this->region2VoxelBasedMeasureGradientImage->nvox, + this->region2VoxelBasedMeasureGradientImage->nbyper); + + //allocate voxel-based gap-overlap peanlty term gradient images + this->gapOverlapGradientWRTDefFieldRegion1 = nifti_copy_nim_info(this->voxelBasedMeasureGradient); + this->gapOverlapGradientWRTDefFieldRegion2 = nifti_copy_nim_info(this->voxelBasedMeasureGradient); + this->gapOverlapGradientWRTDefFieldRegion1->data = (void *)calloc(this->gapOverlapGradientWRTDefFieldRegion1->nvox, + this->gapOverlapGradientWRTDefFieldRegion1->nbyper); + this->gapOverlapGradientWRTDefFieldRegion2->data = (void *)calloc(this->gapOverlapGradientWRTDefFieldRegion2->nvox, + this->gapOverlapGradientWRTDefFieldRegion2->nbyper); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::AllocateVoxelBasedMeasureGradient"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::ClearVoxelBasedMeasureGradient() +{ + //call method from reg_base to clear voxel-based similarity gradient image for combined transform + reg_base::ClearVoxelBasedMeasureGradient(); + + //clear voxel-based similarity gradient images for each region + if (this->region1VoxelBasedMeasureGradientImage != NULL) + { + nifti_image_free(this->region1VoxelBasedMeasureGradientImage); + this->region1VoxelBasedMeasureGradientImage = NULL; + } + if (this->region2VoxelBasedMeasureGradientImage != NULL) + { + nifti_image_free(this->region2VoxelBasedMeasureGradientImage); + this->region2VoxelBasedMeasureGradientImage = NULL; + } + + //clear voxel-based gap-overlap penalty term gradient images + if (this->gapOverlapGradientWRTDefFieldRegion1 != NULL) + { + nifti_image_free(this->gapOverlapGradientWRTDefFieldRegion1); + this->gapOverlapGradientWRTDefFieldRegion1 = NULL; + } + if (this->gapOverlapGradientWRTDefFieldRegion2 != NULL) + { + nifti_image_free(this->gapOverlapGradientWRTDefFieldRegion2); + this->gapOverlapGradientWRTDefFieldRegion2 = NULL; + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ClearVoxelBasedMeasureGradient"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +template +void reg_f3d_sli::AllocateTransformationGradient() +{ + //clear any previously allocated transformation gradients + this->ClearTransformationGradient(); + + //call method from reg_f3d to allocate transformation gradient for region 1 + reg_f3d::AllocateTransformationGradient(); + + //allocate transformation gradient image for region 2 + this->region2TransformationGradient = nifti_copy_nim_info(this->region2ControlPointGrid); + this->region2TransformationGradient->data = (void *)calloc(this->region2TransformationGradient->nvox, + this->region2TransformationGradient->nbyper); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::AllocateTransformationGradient"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::ClearTransformationGradient() +{ + //call method from reg_f3d to clear transformation gradient for region 1 + reg_f3d::ClearTransformationGradient(); + + //clear transformation gradient image for region 2 + if (this->region2TransformationGradient != NULL) + { + nifti_image_free(this->region2TransformationGradient); + this->region2TransformationGradient = NULL; + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ClearTransformationGradient"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +template +T reg_f3d_sli::InitialiseCurrentLevel() +{ + //call method from reg_f3d to calculate max step size for this level and to + //refine gpg for region 1 and modify bending energy weight (and linear energy + //weight?) if required + T maxStepSize = reg_f3d::InitialiseCurrentLevel(); + + // Refine the control point grid for region 2 if required + if (this->gridRefinement && this->currentLevel > 0) + reg_spline_refineControlPointGrid(this->region2ControlPointGrid, this->currentReference); + + //set current distance map + if(this->usePyramid) + this->currentDistanceMap = this->distanceMapPyramid[this->currentLevel]; + else + this->currentDistanceMap = this->distanceMapPyramid[0]; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::InitialiseCurrentLevel"); +#endif + + //return max step size + return maxStepSize; +} +/* *************************************************************** */ +template +void reg_f3d_sli::ClearCurrentInputImage() +{ + //call method from reg_base to clear current reference, floating, and mask image + reg_base::ClearCurrentInputImage(); + + //clear current distance map image + this->currentDistanceMap = NULL; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::ClearCurrentInputImage"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +template +void reg_f3d_sli::SetOptimiser() +{ + //create new optimiser object + //if useConjGradient set then create new conjugate gradient optimiser + if (this->useConjGradient) + this->optimiser = new reg_conjugateGradient(); + //else create standard (gradient ascent) optimiser + else this->optimiser = new reg_optimiser(); + + //initialise optimiser passing pointers to data from transforms and gradients + //from both regions + this->optimiser->Initialise(this->controlPointGrid->nvox,//number of voxels in region 1 CPG + this->controlPointGrid->nz>1 ? 3 : 2, + this->optimiseX, + this->optimiseY, + this->optimiseZ, + this->maxiterationNumber, + 0, // currentIterationNumber + this, //this object, which implements interface for interacting with optimiser + static_cast(this->controlPointGrid->data),//pointer to data from region 1 CPG + static_cast(this->transformationGradient->data),//pointer to data from region 1 gradient + this->region2ControlPointGrid->nvox,//number of voxels in region 2 CPG + static_cast(this->region2ControlPointGrid->data),//pointer to data from region 2 CPG + static_cast(this->region2TransformationGradient->data));//pointer to data from region 2 gradient + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::SetOptimiser"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::UpdateParameters(float scale) +{ + // First update the transformation for region 1 + reg_f3d::UpdateParameters(scale); + + // Create some pointers to the relevant arrays + //noet - call '_b' methods from optimiser to access the region 2 data + T *currentDOFRegion2 = this->optimiser->GetCurrentDOF_b(); + T *bestDOFRegion2 = this->optimiser->GetBestDOF_b(); + T *gradientRegion2 = this->optimiser->GetGradient_b(); + + // update the CPG values for region 2 + size_t voxNumberRegion2 = this->optimiser->GetVoxNumber_b(); + // Update the values for the x-axis displacement + if (this->optimiser->GetOptimiseX() == true) + { + for (size_t i = 0; i < voxNumberRegion2; ++i) + { + currentDOFRegion2[i] = bestDOFRegion2[i] + scale * gradientRegion2[i]; + } + } + // Update the values for the y-axis displacement + if (this->optimiser->GetOptimiseY() == true) + { + //pointers to y-axis displacements + T *currentDOFYRegion2 = ¤tDOFRegion2[voxNumberRegion2]; + T *bestDOFYRegion2 = &bestDOFRegion2[voxNumberRegion2]; + T *gradientYRegion2 = &gradientRegion2[voxNumberRegion2]; + for (size_t i = 0; i < voxNumberRegion2; ++i) + { + currentDOFYRegion2[i] = bestDOFYRegion2[i] + scale * gradientYRegion2[i]; + } + } + // Update the values for the z-axis displacement + if (this->optimiser->GetOptimiseZ() == true && this->optimiser->GetNDim()>2) + { + //pointers to z-axis displacements + T *currentDOFZRegion2 = ¤tDOFRegion2[2 * voxNumberRegion2]; + T *bestDOFZRegion2 = &bestDOFRegion2[2 * voxNumberRegion2]; + T *gradientZRegion2 = &gradientRegion2[2 * voxNumberRegion2]; + for (size_t i = 0; i < voxNumberRegion2; ++i) + { + currentDOFZRegion2[i] = bestDOFZRegion2[i] + scale * gradientZRegion2[i]; + } + } + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::UpdateParameters"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::UpdateBestObjFunctionValue() +{ + // call method from reg_f3d to update all values except gap-overlap term + reg_f3d::UpdateBestObjFunctionValue(); + //now update best gap-overlap value + this->bestWGO = this->currentWGO; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::UpdateBestObjFunctionValue"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::PrintInitialObjFunctionValue() +{ + // if verbose not set don't display anything + if (!this->verbose) return; + + // format text with initial total objective function value + char text[255]; + sprintf(text, "Initial objective function: %g", this->optimiser->GetBestObjFunctionValue()); + // and similarity measure(s) value + sprintf(text + strlen(text), " = (wSIM)%g", this->bestWMeasure); + // and values of penalty terms + if (this->bendingEnergyWeight > 0) + sprintf(text + strlen(text), " - (wBE)%.2e", this->bestWBE); + if (this->linearEnergyWeight > 0) + sprintf(text + strlen(text), " - (wLE)%.2e", this->bestWLE); + //jacobian and landmark penalty terms not currently implemented for sliding region registrations + //if (this->jacobianLogWeight>0) + // sprintf(text + strlen(text), " - (wJAC)%.2e", this->bestWJac); + //if (this->landmarkRegWeight>0) + // sprintf(text + strlen(text), " - (wLAN)%.2e", this->bestWLand); + if (this->gapOverlapWeight > 0) + sprintf(text + strlen(text), " - (wGO)%.2e", this->bestWGO); + + //display text + reg_print_info(this->executableName, text); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::PrintInitialObjFunctionValue"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::PrintCurrentObjFunctionValue(T currentSize) +{ + // if verbose not set don't display anything + if (!this->verbose) return; + + // format text with iteration number and current total objective function value + char text[255]; + sprintf(text, "[%i] Current objective function: %g", + (int)this->optimiser->GetCurrentIterationNumber(), + this->optimiser->GetBestObjFunctionValue()); + // and similarity measure(s) value + sprintf(text + strlen(text), " = (wSIM)%g", this->bestWMeasure); + // and values of penalty terms + if (this->bendingEnergyWeight > 0) + sprintf(text + strlen(text), " - (wBE)%.2e", this->bestWBE); + if (this->linearEnergyWeight > 0) + sprintf(text + strlen(text), " - (wLE)%.2e", this->bestWLE); + //jacobian and landmark penalty terms not currently implemented for sliding region registrations + //if (this->jacobianLogWeight>0) + // sprintf(text + strlen(text), " - (wJAC)%.2e", this->bestWJac); + //if (this->landmarkRegWeight>0) + // sprintf(text + strlen(text), " - (wLAN)%.2e", this->bestWLand); + if (this->gapOverlapWeight > 0) + sprintf(text + strlen(text), " - (wGO)%.2e", this->bestWGO); + //add current step size to text + sprintf(text + strlen(text), " [+ %g mm]", currentSize); + + //display text + reg_print_info(this->executableName, text); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::PrintCurrentObjFunctionValue"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +template +void reg_f3d_sli::SetDistanceMapImage(nifti_image *distanceMapImage) +{ + this->inputDistanceMap = distanceMapImage; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::SetDistanceMapImage"); +#endif +} +/* *************************************************************** */ +template +void reg_f3d_sli::SetGapOverlapWeight(T weight) +{ + this->gapOverlapWeight = weight; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::SetGapOverlapWeight"); +#endif +} +/* *************************************************************** */ +template +nifti_image * reg_f3d_sli::GetRegion2ControlPointPositionImage() +{ + // Create a control point grid nifti image + nifti_image *returnedControlPointGrid = nifti_copy_nim_info(this->region2ControlPointGrid); + + // Allocate the new image data array + returnedControlPointGrid->data = (void *)malloc(returnedControlPointGrid->nvox*returnedControlPointGrid->nbyper); + + // Copy the final region2 control point grid image + memcpy(returnedControlPointGrid->data, this->region2ControlPointGrid->data, + returnedControlPointGrid->nvox * returnedControlPointGrid->nbyper); + + // Return the new control point grid +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetRegion2ControlPointPositionImage"); +#endif + return returnedControlPointGrid; +} +/* *************************************************************** */ +template +void reg_f3d_sli::SetRegion2ControlPointGridImage(nifti_image *controlPointGrid) +{ + this->inputRegion2ControlPointGrid = controlPointGrid; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::SetRegion2ControlPointGridImage"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +template +void reg_f3d_sli::CheckParameters() +{ + //call method from reg_f3d to check standard parameters + reg_f3d::CheckParameters(); + + //check the distance map has been defined + if (this->inputDistanceMap == NULL) + { + reg_print_fct_error("reg_f3d_sli::CheckParameters()"); + reg_print_msg_error("The distance map image is not defined"); + reg_exit(); + } + else + { + //and has the same dimensions as the floating (source) image + if (this->inputDistanceMap->nx != this->inputFloating->nx || + this->inputDistanceMap->ny != this->inputFloating->ny || + this->inputDistanceMap->nz != this->inputFloating->nz) + { + reg_print_fct_error("reg_f3d_sli::CheckParameters()"); + reg_print_msg_error("The distance map has different dimensions to the floating image"); + reg_exit(); + } + } + + //check if an input CPG has only been provided for one region + if ((this->inputControlPointGrid != NULL && this->inputRegion2ControlPointGrid == NULL) || + (this->inputControlPointGrid == NULL && this->inputRegion2ControlPointGrid != NULL)) + { + reg_print_fct_error("reg_f3d_sli::CheckParameters()"); + reg_print_msg_error("An input Control Point Grid has only been provided for one region"); + reg_print_msg_error("You must provide a Control Point Grid for both regions (or none)"); + reg_exit(); + } + + //if input CPGs provided for both regions check they have the same dimensions + if (this->inputControlPointGrid != NULL && this->inputRegion2ControlPointGrid != NULL) + { + if (this->inputControlPointGrid->nx != this->inputRegion2ControlPointGrid->nx || + this->inputControlPointGrid->ny != this->inputRegion2ControlPointGrid->ny || + this->inputControlPointGrid->nz != this->inputRegion2ControlPointGrid->nz) + { + reg_print_fct_error("reg_f3d_sli::CheckParameters()"); + reg_print_msg_error("The input Control Point Grids for the two regions have different dimensions"); + reg_exit(); + } + } + + + //check if jacobian or landmark penalty term weights have been set - if so throw error as + //these terms are not yet implemented for sliding region registrations + if (this->jacobianLogWeight > 0) + { + reg_print_fct_error("reg_f3d_sli::CheckParameters()"); + reg_print_msg_error("Jacobian penalty term weight > 0"); + reg_print_msg_error("Jacobian penalty term has not yet been implemented to work with sliding region registrations"); + reg_exit(); + } + if (this->landmarkRegWeight > 0) + { + reg_print_fct_error("reg_f3d_sli::CheckParameters()"); + reg_print_msg_error("Landmark penalty term weight > 0"); + reg_print_msg_error("Landmark penalty term has not yet been implemented to work with sliding region registrations"); + reg_exit(); + } + + // check if penalty term weights >= 1 + T penaltySum = this->bendingEnergyWeight + this->linearEnergyWeight + this->gapOverlapWeight; + if (penaltySum >= 1) + { + //display a warning saying images will be ignored for the registration + reg_print_fct_warn("reg_f3d_sli::CheckParameters()"); + reg_print_msg_warn("Penalty term weights greater than or equal to 1"); + reg_print_msg_warn("The images will be ignored during the registration "); + this->similarityWeight = 0; + this->bendingEnergyWeight /= penaltySum; + this->linearEnergyWeight /= penaltySum; + this->gapOverlapWeight /= penaltySum; + } + else this->similarityWeight = 1.0 - penaltySum; + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::CheckParameters"); +#endif + return; +} +/* *************************************************************** */ +template +void reg_f3d_sli::Initialise() +{ + //call method from reg_f3d to initialise image pyramids and region 1 control point grid + reg_f3d::Initialise(); + + //initialise control point grid for region 2 + // + //check if an input grid has been provided + if (this->inputRegion2ControlPointGrid != NULL) + { + //if so use input grid to initialise region 2 control point grid + this->region2ControlPointGrid = nifti_copy_nim_info(this->inputRegion2ControlPointGrid); + this->region2ControlPointGrid->data = (void *)malloc(this->region2ControlPointGrid->nvox * this->region2ControlPointGrid->nbyper); + memcpy(this->region2ControlPointGrid->data, this->inputRegion2ControlPointGrid->data, + this->region2ControlPointGrid->nvox * this->region2ControlPointGrid->nbyper); + } + else + { + //if not copy grid from region 1 + this->region2ControlPointGrid = nifti_copy_nim_info(this->controlPointGrid); + this->region2ControlPointGrid->data = (void *)malloc(this->region2ControlPointGrid->nvox * this->region2ControlPointGrid->nbyper); + memcpy(this->region2ControlPointGrid->data, this->controlPointGrid->data, + this->region2ControlPointGrid->nvox * this->region2ControlPointGrid->nbyper); + } + + //check if image pyramids are being used for multi-resolution + if (this->usePyramid) + { + //create image pyramid for distance map, with one image for each resolution level + this->distanceMapPyramid = (nifti_image **)malloc(this->levelToPerform * sizeof(nifti_image *)); + reg_createImagePyramid(this->inputDistanceMap, this->distanceMapPyramid, this->levelNumber, this->levelToPerform); + + //divide each distance map in pyramid by sqrt of number of voxel to normalise gap-overlap penalty term + for (int n = 0; n < this->levelToPerform; n++) + { + float numVox = this->distanceMapPyramid[n]->nx * this->distanceMapPyramid[n]->ny * this->distanceMapPyramid[n]->nz; + reg_tools_divideValueToImage(this->distanceMapPyramid[n], this->distanceMapPyramid[n], sqrt(numVox)); + } + } + else + { + //image pyramids are not used, so create pyramid with just one level (i.e. copy of input image) + this->distanceMapPyramid = (nifti_image **)malloc(sizeof(nifti_image *)); + reg_createImagePyramid(this->inputDistanceMap, this->distanceMapPyramid, 1, 1); + + //divide distance map in pyramid by sqrt of number of voxel to normalise gap-overlap penalty term + float numVox = this->distanceMapPyramid[0]->nx * this->distanceMapPyramid[0]->ny * this->distanceMapPyramid[0]->nz; + reg_tools_divideValueToImage(this->distanceMapPyramid[0], this->distanceMapPyramid[0], sqrt(numVox)); + } + +#ifdef NDEBUG + if (this->verbose) + { +#endif + //print out some info: + std::string text; + //name of distance map image + text = stringFormat("Distance map image used for sliding regions: %s", this->inputDistanceMap->fname); + reg_print_info(this->executableName, text.c_str()); + //weight of gap-overlap penalty term + text = stringFormat("Gap-overlap penalty term weight: %g", this->gapOverlapWeight); + reg_print_info(this->executableName, text.c_str()); + reg_print_info(this->executableName, ""); + +#ifdef NDEBUG + } +#endif +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d::Initialise"); +#endif + +} +/* *************************************************************** */ +template +nifti_image **reg_f3d_sli::GetWarpedImage() +{ + //the input distance map image is used + if (this->inputDistanceMap == NULL) + { + reg_print_fct_error("reg_f3d_sli::GetWarpedImage()"); + reg_print_msg_error("The distance map image has to be defined"); + reg_exit(); + } + this->currentDistanceMap = this->inputDistanceMap; + + //call method from reg_f3d to get image + nifti_image **warpedImage = reg_f3d::GetWarpedImage(); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetWarpedImage"); +#endif + + //and return image + return warpedImage; +} +/* *************************************************************** */ +/* *************************************************************** */ +template class reg_f3d_sli; +#endif diff --git a/reg-lib/_reg_f3d_sli.h b/reg-lib/_reg_f3d_sli.h new file mode 100644 index 00000000..65e955cf --- /dev/null +++ b/reg-lib/_reg_f3d_sli.h @@ -0,0 +1,216 @@ +/* +* @file _reg_f3d_sli.h +* @author Jamie McClelland +* @date 20/08/2017 +* +* Copyright (c) 2017, University College London. All rights reserved. +* Centre for Medical Image Computing (CMIC) +* See the LICENSE.txt file in the nifty_reg root folder +* +*/ + +#ifndef _REG_F3D_SLI_H +#define _REG_F3D_SLI_H + +#include "_reg_f3d.h" + +/// @brief a class for Fast Free Form Deformation registration wiht sliding regions + +/* +FFD registration that allows for two sliding regions, defined in the source image +a distance map (in the source image) must be provided, giving the distance to the +sliding interface for all voxels, with voxels in region1 having negative values, +and voxels in region2 having positive values +The deformation in each region is represented by a separate B-spline transformation +and a penalty term is used to penalise gaps/overlaps between the two regions. +*/ +template +class reg_f3d_sli : public reg_f3d +{ +protected: + //variables for region2 transform + nifti_image *inputRegion2ControlPointGrid; //pointer to external + nifti_image *region2ControlPointGrid; + nifti_image *region2DeformationFieldImage; + nifti_image *region2VoxelBasedMeasureGradientImage; + nifti_image *region2TransformationGradient; + + //variables for region1 transform + nifti_image *region1DeformationFieldImage; + nifti_image *region1VoxelBasedMeasureGradientImage; + + //variables for distance map image + nifti_image *inputDistanceMap; //pointer to external + nifti_image **distanceMapPyramid; + nifti_image *currentDistanceMap; + + //variables for the distance map warped by the transform for each region + nifti_image *warpedDistanceMapRegion1; + nifti_image *warpedDistanceMapRegion2; + + //variables for the spatial gradient of distance map warped by the transform + //for each region + nifti_image *warpedDistanceMapGradientRegion1; + nifti_image *warpedDistanceMapGradientRegion2; + + //variables for gap-overlap penalty term + T gapOverlapWeight; + double currentWGO; + double bestWGO; + + //variables for the gradient of the penalty term with respect to the def field + //for each region + nifti_image *gapOverlapGradientWRTDefFieldRegion1; + nifti_image *gapOverlapGradientWRTDefFieldRegion2; + + + //reimplement method to get deformation field + //combines deformation fields from each region based on warped distance maps + //at each voxel, if sum of warped distance maps < 0 use region 1 def field else + //use region 2 def field. if one warped distance map has a value of NaN (due to + //transform mapping outside distance map image) and other warped distance map + //maps to the same region as used to warp it (i.e. region 1 warped distance map + // < 0 or region 2 warped distance map >= 0) then use def field from non-NaN + //region, else combined def field set to NaN. If both warped distance maps are + //NaN then combined def field set to NaN. + virtual void GetDeformationField(); + //reimplement methods to allocate/clear deformation field + //these methods will allocate/clear the deformation fields for each region as well + //as the combined deformation field. + virtual void AllocateDeformationField(); + virtual void ClearDeformationField(); + //reimplement methods to allocate/clear the warped images so that the warped + //distance maps are also allocated/cleared + virtual void AllocateWarped(); + virtual void ClearWarped(); + + + //methods to calculate objective function + //note - no need to reimplement method to get similarity measure value as ths will + //be calculated using the combined deformation field using the existing methods + // + //reimplement method to calculate objective function value to also include + //value of gap-overlap penalty term + virtual double GetObjectiveFunctionValue(); + //new method to calculate gap-overlap penalty term + //at each voxel, the warped distance maps are multiplied together - if the result + //is less than 0 (indicating the transforms for the two regions map the voxel into + //different regiions, i.e. a gap or overlap is present) then the negative of the + //result is added to the gap-overlap penalty term. + virtual double ComputeGapOverlapPenaltyTerm(); + //reimplement methods to calculate prenalty terms using transforms for both regions + virtual double ComputeBendingEnergyPenaltyTerm(); + virtual double ComputeLinearEnergyPenaltyTerm(); + //Jacobian penalty term not currently implemented to work with sliding region registrations + //this method will throw an error if called + virtual double ComputeJacobianBasedPenaltyTerm(int); + //Landmark distance penalty term not currently implemented to work with sliding region registrations + //this method will throw an error if called + virtual double ComputeLandmarkDistancePenaltyTerm(); + + + //methods to calculate objective function gradient + // + //reimplement method to calculate objective function gradient to include gradient of + //gap-ovlerlap penalty term + virtual void GetObjectiveFunctionGradient(); + //reimplement method to convert voxel-based similarity gradient to CPG based + //gradient(s). splits voxel-based gradient between two regions, based on warped + //distance maps, and then converts voxel-based gradient for each region to CPG + //gradients + virtual void GetSimilarityMeasureGradient(); + //new method to calculate the gap-overlap penalty term gradient + virtual void GetGapOverlapGradient(); + //reimplement methods to calculate penalty term gradients for transforms for both regions + virtual void GetBendingEnergyGradient(); + virtual void GetLinearEnergyGradient(); + //Jacobian penalty term not currently implemented to work with sliding region registrations + //this method will throw an error if called + virtual void GetJacobianBasedGradient(); + //Landmark distance penalty term not currently implemented to work with sliding region registrations + //this method will throw an error if called + virtual void GetLandmarkDistanceGradient(); + //reimplement method to set gradient image to zero to set gradient images for both regions to 0 + virtual void SetGradientImageToZero(); + //reimplement method to normalise gradient so that gradients for both regions are normalised + //using the max value over both gradient images + virtual T NormaliseGradient(); + //reimplement method to smooth gradient so that gradients for both regions are smoothed + virtual void SmoothGradient(); + //remiplement method to approximate gradient so that gradients for both regions are approximated + virtual void GetApproximatedGradient(); + //reimplement methods to allocate/clear warped gradient images so that the warped + //distance map gradients are also allocated/cleared + virtual void AllocateWarpedGradient(); + virtual void ClearWarpedGradient(); + //reimplement methods to allocate/clear 'voxel-based' similarity measure gradient + //image (i.e. the similarity measure gradient WRT the def field) - these methods + //will now also allocate/clear the 'voxel-based' similarity measure gradients for + //each region and the 'voxel-based' gap-overlap penalty term gradient images + virtual void AllocateVoxelBasedMeasureGradient(); + virtual void ClearVoxelBasedMeasureGradient(); + //reimplement methods to allocate/clear transformation gradient images - these methods + //will now also allocate/clear the transformation gradient image for region 2 + virtual void AllocateTransformationGradient(); + virtual void ClearTransformationGradient(); + + //reimplement method to initialise current level to refine CPGs for transforms for both + //regions and to set the current distance map image + virtual T InitialiseCurrentLevel(); + //reimplement method to clear current input images to also clear current distance map + //image + virtual void ClearCurrentInputImage(); + + //reimplement method for setting optimiser so that region 2 transform data and gradient + //data also passed to optimiser. + //note - no modifications to optimiser required as it can already jointly optimise 2 + //transforms for use with symmetric registrations + virtual void SetOptimiser(); + //reimplement method for updating parameters so that region 2 transform is updated as well + virtual void UpdateParameters(float); + //reimplement method for updating best objective function value so that gap-overlap value + //is updated as well + virtual void UpdateBestObjFunctionValue(); + //reimplement methods for printing objective function value so that gap-overlap value is + //also printed + virtual void PrintInitialObjFunctionValue(); + virtual void PrintCurrentObjFunctionValue(T); + +public: + //constructor and destructor methods + reg_f3d_sli(int refTimePoint, int floTimePoint); + ~reg_f3d_sli(); + + + //new method to set distance map image + virtual void SetDistanceMapImage(nifti_image *); + + //new method to set gap-overlap penalty term weight + virtual void SetGapOverlapWeight(T); + + //new methods to get and set transform for region 2 + //note - used similar method names as for methods for region 1 (i.e. standard methods from reg_f3d) + //hence get method called ...Position... and set method called ...Grid... + virtual nifti_image *GetRegion2ControlPointPositionImage(); + virtual void SetRegion2ControlPointGridImage(nifti_image *); + + + //reimplement method to check parameters so that also checks if distance map has been set + //and has same dimensions as floating image. + //Also checks if an input control point grid has been set for one region but not the other, + //and throws an error if so. If input control point grids have been set for both regions + //then checks they have the same dimensions. + //And checks that jacobian and landmark penalty terms have not been set (as not yet + //implemented for sliding region registrations) and normalises penalty term weights + virtual void CheckParameters(); + + //reimplement method to initialise registration so that also initialises CPG for region 2 + //and image pyramid for distance map image + virtual void Initialise(); + + //reimplement method to get warped image so that distance map set to input distance map + //before calling method from ref_f3d + virtual nifti_image **GetWarpedImage(); +}; + +#endif diff --git a/reg-lib/cpu/_reg_lncc.cpp b/reg-lib/cpu/_reg_lncc.cpp index 7b775364..34d362df 100644 --- a/reg-lib/cpu/_reg_lncc.cpp +++ b/reg-lib/cpu/_reg_lncc.cpp @@ -106,8 +106,8 @@ void reg_lncc::UpdateLocalStatImages(nifti_image *refImage, size_t voxelNumber = (size_t)refImage->nx*refImage->ny*refImage->nz; #endif memcpy(combinedMask, refMask, voxelNumber*sizeof(int)); - reg_tools_removeNanFromMask(refImage, combinedMask); - reg_tools_removeNanFromMask(warImage, combinedMask); + reg_tools_removeNanFromMask(refImage, combinedMask,current_timepoint); + reg_tools_removeNanFromMask(warImage, combinedMask,current_timepoint); DTYPE *origRefPtr = static_cast(refImage->data); DTYPE *meanRefPtr = static_cast(meanRefImage->data); diff --git a/reg-lib/cpu/_reg_mind.cpp b/reg-lib/cpu/_reg_mind.cpp index 294d21c3..e3f9866b 100644 --- a/reg-lib/cpu/_reg_mind.cpp +++ b/reg-lib/cpu/_reg_mind.cpp @@ -136,7 +136,7 @@ void GetMINDImageDesciptor_core(nifti_image* inputImage, #if defined (_OPENMP) #pragma omp parallel for default(none) \ shared(samplingNbr, maskPtr, meanImgDataPtr, \ - MINDImgDataPtr) \ + MINDImgDataPtr, voxelNumber) \ private(voxelIndex, meanValue, max_desc, descValue, mindIndex) #endif for(voxelIndex=0;voxelIndexny * referenceImagePointer->nz; int *combinedMask = (int *)malloc(voxelNumber*sizeof(int)); memcpy(combinedMask, this->referenceMaskPointer, voxelNumber*sizeof(int)); - reg_tools_removeNanFromMask(this->referenceImagePointer, combinedMask); - reg_tools_removeNanFromMask(this->warpedFloatingImagePointer, combinedMask); + reg_tools_removeNanFromMask(this->referenceImagePointer, combinedMask, t); + reg_tools_removeNanFromMask(this->warpedFloatingImagePointer, combinedMask, t); if(this->mind_type==MIND_TYPE){ GetMINDImageDesciptor(this->referenceImagePointer, @@ -578,8 +578,8 @@ double reg_mind::GetSimilarityMeasureValue() floatingImagePointer->ny * floatingImagePointer->nz; combinedMask = (int *)malloc(voxelNumber*sizeof(int)); memcpy(combinedMask, this->floatingMaskPointer, voxelNumber*sizeof(int)); - reg_tools_removeNanFromMask(this->floatingImagePointer, combinedMask); - reg_tools_removeNanFromMask(this->warpedReferenceImagePointer, combinedMask); + reg_tools_removeNanFromMask(this->floatingImagePointer, combinedMask,t); + reg_tools_removeNanFromMask(this->warpedReferenceImagePointer, combinedMask,t); if(this->mind_type==MIND_TYPE){ GetMINDImageDesciptor(this->floatingImagePointer, @@ -655,8 +655,8 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int current_timepoint) this->referenceImagePointer->nz; int *combinedMask = (int *)malloc(voxelNumber*sizeof(int)); memcpy(combinedMask, this->referenceMaskPointer, voxelNumber*sizeof(int)); - reg_tools_removeNanFromMask(this->referenceImagePointer, combinedMask); - reg_tools_removeNanFromMask(this->warpedFloatingImagePointer, combinedMask); + reg_tools_removeNanFromMask(this->referenceImagePointer, combinedMask,current_timepoint); + reg_tools_removeNanFromMask(this->warpedFloatingImagePointer, combinedMask,current_timepoint); if(this->mind_type==MIND_TYPE){ // Compute the reference image descriptors @@ -740,8 +740,8 @@ void reg_mind::GetVoxelBasedSimilarityMeasureGradient(int current_timepoint) floatingImagePointer->ny * floatingImagePointer->nz; combinedMask = (int *)malloc(voxelNumber*sizeof(int)); memcpy(combinedMask, this->floatingMaskPointer, voxelNumber*sizeof(int)); - reg_tools_removeNanFromMask(this->floatingImagePointer, combinedMask); - reg_tools_removeNanFromMask(this->warpedReferenceImagePointer, combinedMask); + reg_tools_removeNanFromMask(this->floatingImagePointer, combinedMask,current_timepoint); + reg_tools_removeNanFromMask(this->warpedReferenceImagePointer, combinedMask,current_timepoint); if(this->mind_type==MIND_TYPE){ GetMINDImageDesciptor(this->floatingImagePointer, diff --git a/reg-lib/cpu/_reg_resampling.cpp b/reg-lib/cpu/_reg_resampling.cpp index eb72c89b..73f6fc17 100755 --- a/reg-lib/cpu/_reg_resampling.cpp +++ b/reg-lib/cpu/_reg_resampling.cpp @@ -443,7 +443,10 @@ void ResampleImage3D(nifti_image *floatingImage, intensity=paddingValue; - if((maskPtr[index])>-1) + if((maskPtr[index])>-1 && + deformationFieldPtrX[index] == deformationFieldPtrX[index] && + deformationFieldPtrY[index] == deformationFieldPtrY[index] && + deformationFieldPtrZ[index] == deformationFieldPtrZ[index]) { world[0]=static_cast(deformationFieldPtrX[index]); world[1]=static_cast(deformationFieldPtrY[index]); @@ -644,7 +647,9 @@ void ResampleImage2D(nifti_image *floatingImage, { intensity=paddingValue; - if((maskPtr[index])>-1) + if ((maskPtr[index])>-1 && + deformationFieldPtrX[index] == deformationFieldPtrX[index] && + deformationFieldPtrY[index] == deformationFieldPtrY[index]) { world[0] = static_cast(deformationFieldPtrX[index]); world[1] = static_cast(deformationFieldPtrY[index]); diff --git a/reg-lib/cpu/_reg_ssd.cpp b/reg-lib/cpu/_reg_ssd.cpp index 98d1e206..7b65154a 100755 --- a/reg-lib/cpu/_reg_ssd.cpp +++ b/reg-lib/cpu/_reg_ssd.cpp @@ -823,7 +823,7 @@ void GetDiscretisedValueSSD_core3D_2(nifti_image *controlPointGridImage, #pragma omp parallel for default(none) \ shared(controlPointGridImage, refImage, warImage, grid2img_vox, blockSize, \ padding_value, refBlockValue, mask, refImgPtr, warImgPtr, discretise_radius, \ - discretise_step, discretisedValue) \ + discretise_step, discretisedValue, voxelBlockNumber_t, voxelNumber, voxelBlockNumber, label_nD_number) \ private(cpx, cpy, cpz, x, y, z, a, b, c, t, currentControlPoint, gridVox, imageVox, \ voxIndex, idBlock, blockIndex, definedValueNumber, tid, \ timeV, voxIndex_t, blockIndex_t, discretisedIndex, currentSum, currentValue) diff --git a/reg-lib/cpu/_reg_tools.cpp b/reg-lib/cpu/_reg_tools.cpp index 247c7327..e65fbbdd 100755 --- a/reg-lib/cpu/_reg_tools.cpp +++ b/reg-lib/cpu/_reg_tools.cpp @@ -2386,31 +2386,34 @@ int reg_tools_nanMask_image(nifti_image *image, nifti_image *maskImage, nifti_im /* *************************************************************** */ /* *************************************************************** */ template -int reg_tools_removeNanFromMask_core(nifti_image *image, int *mask) +int reg_tools_removeNanFromMask_core(nifti_image *image, int *mask, int time_point_index) { size_t voxelNumber = (size_t)image->nx*image->ny*image->nz; TYPE *imagePtr = static_cast(image->data); - for(int t=0; tnt; ++t){ - for(size_t i=0; idatatype) { case NIFTI_TYPE_FLOAT32: return reg_tools_removeNanFromMask_core - (image, mask); + (image, mask, time_point_index); case NIFTI_TYPE_FLOAT64: return reg_tools_removeNanFromMask_core - (image, mask); + (image, mask, time_point_index); default: reg_print_fct_error("reg_tools_removeNanFromMask"); reg_print_msg_error("The image data type is not supported"); diff --git a/reg-lib/cpu/_reg_tools.h b/reg-lib/cpu/_reg_tools.h index 692818ca..08da1db7 100755 --- a/reg-lib/cpu/_reg_tools.h +++ b/reg-lib/cpu/_reg_tools.h @@ -296,7 +296,7 @@ int reg_tools_nanMask_image(nifti_image *img, */ extern "C++" int reg_tools_removeNanFromMask(nifti_image *image, - int *mask); + int *mask, int time_point_index); /* *************************************************************** */ /** @brief Get the minimal value of an image * @param img Input image