From fade7570e9ffa4b188f7f1c9972253d70a796967 Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Mon, 21 Aug 2017 15:31:55 +0100 Subject: [PATCH 01/16] initial commit for sliding regions created reg_f3d_sli .h and .cpp files updated CMakeLists to build new class --- niftyreg_build_version.txt | 2 +- reg-lib/CMakeLists.txt | 3 + reg-lib/_reg_f3d_sli.cpp | 127 +++++++++++++++++++++++++++++++++++++ reg-lib/_reg_f3d_sli.h | 67 +++++++++++++++++++ 4 files changed, 198 insertions(+), 1 deletion(-) create mode 100644 reg-lib/_reg_f3d_sli.cpp create mode 100644 reg-lib/_reg_f3d_sli.h diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 920a1396..c739b42c 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -43 +44 diff --git a/reg-lib/CMakeLists.txt b/reg-lib/CMakeLists.txt index 51c49e34..f89c6e0e 100755 --- a/reg-lib/CMakeLists.txt +++ b/reg-lib/CMakeLists.txt @@ -195,6 +195,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 @@ -215,6 +217,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) #----------------------------------------------------------------------------- #----------------------------------------------------------------------------- diff --git a/reg-lib/_reg_f3d_sli.cpp b/reg-lib/_reg_f3d_sli.cpp new file mode 100644 index 00000000..3bb9559d --- /dev/null +++ b/reg-lib/_reg_f3d_sli.cpp @@ -0,0 +1,127 @@ +/* +* _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 SLI"; + + this->region2ControlPointGrid = NULL; + this->region2DeformationFieldImage = NULL; + this->region2VoxelBasedMeasureGradientImage = NULL; + this->region2TransformationGradient = NULL; + + this->region1DeformationFieldImage = NULL; + this->region1VoxelBasedMeasureGradientImage = NULL; + + this->distanceMapImage = NULL; + this->distanceMapPyramid = NULL; + this->currentDistanceMap = NULL; + this->warpedDistanceMapRegion1 = NULL; + this->warpedDistanceMapRegion2 = NULL; + + this->gapOverlapWeight = 0.1; + +#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 + reg_spline_getDeformationField(this->controlPointGrid, + this->region1DeformationFieldImage, + this->currentMask, + false, //composition + true); //bspline + //get deformation field for region2 + reg_spline_getDeformationField(this->region2ControlPointGrid, + this->region2DeformationFieldImage, + this->currentMask, + false, //composition + true); //bspline + + //warp distance map using region1 def field + reg_resampleImage(this->currentDistanceMap, + this->warpedDistanceMapRegion1, + this->region1DeformationFieldImage, + this->currentMask, + 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, + this->currentMask, + 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 + size_t numVox = this->region1DeformationFieldImage->nx * + this->region1DeformationFieldImage->ny * + this->region1DeformationFieldImage->nz; + //need pointers + +} +/* *************************************************************** */ +/* *************************************************************** */ +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..9254708b --- /dev/null +++ b/reg-lib/_reg_f3d_sli.h @@ -0,0 +1,67 @@ +/* +* @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 *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 *distanceMapImage; + nifti_image **distanceMapPyramid; + nifti_image *currentDistanceMap; + nifti_image *warpedDistanceMapRegion1; + nifti_image *warpedDistanceMapRegion2; + + //variables for penalty term + T gapOverlapWeight; + double currentWGO; + double bestWGO; + + + //reimplement function to get deformation field + //combines deformation fields from each region based on warped distance maps + virtual void GetDeformationField(); + + //new functions for Gap-Overlap penalty term + //virtual double GetGapOverlapPenaltyTerm(); + //virtual void GetGapOverlapGradient(); + +public: + reg_f3d_sli(int refTimePoint, int floTimePoint); + ~reg_f3d_sli(); +}; + +#endif From bee7e41419f107be2d4605ee583e56e548b79ddf Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Tue, 22 Aug 2017 18:45:50 +0100 Subject: [PATCH 02/16] implemented GetDeformationField() and GetSimilarityMeasureGradient() methods --- niftyreg_build_version.txt | 2 +- reg-lib/_reg_f3d_sli.cpp | 361 ++++++++++++++++++++++++++++++++++++- reg-lib/_reg_f3d_sli.h | 10 +- 3 files changed, 369 insertions(+), 4 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index c739b42c..ea90ee31 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -44 +45 diff --git a/reg-lib/_reg_f3d_sli.cpp b/reg-lib/_reg_f3d_sli.cpp index 3bb9559d..4064b7b5 100644 --- a/reg-lib/_reg_f3d_sli.cpp +++ b/reg-lib/_reg_f3d_sli.cpp @@ -118,8 +118,367 @@ void reg_f3d_sli::GetDeformationField() size_t numVox = this->region1DeformationFieldImage->nx * this->region1DeformationFieldImage->ny * this->region1DeformationFieldImage->nz; - //need pointers + //pointers to deformation fields + T *region1DFPtrX = static_cast(this->region1DeformationFieldImage->data); + T *region1DFPtrY = ®ion1DFPtrX[numVox]; + T *region2DFPtrX = static_cast(this->region2DeformationFieldImage->data); + T *region2DFPtrY = ®ion2DFPtrX[numVox]; + T *combinedDFPtrX = static_cast(this->deformationFieldImage->data); + T *combinedDFPtrY = &combinedDFPtrX[numVox]; + //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 + T *region1DFPtrZ = ®ion1DFPtrY[numVox]; + T *region2DFPtrZ = ®ion2DFPtrY[numVox]; + T *combinedDFPtrZ = &combinedDFPtrY[numVox]; + + //loop over voxels + for (size_t n = 0; n < numVox; n++) + { + //check in mask + if (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(); + 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(); + combinedDFPtrZ[n] = std::numeric_limits::quiet_NaN(); + } + else + { + //set combined def field to region2 def field + combinedDFPtrX[n] = region2DFPtrX[n]; + combinedDFPtrY[n] = region2DFPtrY[n]; + 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(); + combinedDFPtrZ[n] = std::numeric_limits::quiet_NaN(); + } + else + { + //set combined def field to region1 def field + combinedDFPtrX[n] = region1DFPtrX[n]; + combinedDFPtrY[n] = region1DFPtrY[n]; + 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]; + combinedDFPtrZ[n] = region1DFPtrZ[n]; + } + else + { + combinedDFPtrX[n] = region2DFPtrX[n]; + combinedDFPtrY[n] = region2DFPtrY[n]; + 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 + else + { + combinedDFPtrX[n] = std::numeric_limits::quiet_NaN(); + combinedDFPtrY[n] = std::numeric_limits::quiet_NaN(); + combinedDFPtrZ[n] = std::numeric_limits::quiet_NaN(); + } + }//for (size_t n = 0; n < numVox; n++) + }//if (this->region1DeformationFieldImage->nz > 1) + else + { + //images are 2D + //loop over voxels + for (size_t n = 0; n < numVox; n++) + { + //check in mask + if (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(); + } + 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(); + } + else + { + //set combined def field to region2 def field + combinedDFPtrX[n] = region2DFPtrX[n]; + combinedDFPtrY[n] = region2DFPtrY[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(); + } + else + { + //set combined def field to region1 def field + combinedDFPtrX[n] = region1DFPtrX[n]; + combinedDFPtrY[n] = region1DFPtrY[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]; + } + else + { + combinedDFPtrX[n] = region2DFPtrX[n]; + combinedDFPtrY[n] = region2DFPtrY[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 + else + { + combinedDFPtrX[n] = std::numeric_limits::quiet_NaN(); + combinedDFPtrY[n] = std::numeric_limits::quiet_NaN(); + } + }//for (size_t n = 0; n < numVox; n++) + }//else (this->region1DeformationFieldImage->nz > 1) + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetDeformationField()"); +#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 + T *combinedVBMGPtr = static_cast(this->voxelBasedMeasureGradient->data); + T *region1VBMGPtr = static_cast(this->region1VoxelBasedMeasureGradientImage->data); + T *region2VBMGPtr = static_cast(this->region2VoxelBasedMeasureGradientImage->data); + + //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 + size_t numVox = this->voxelBasedMeasureGradient->nx * + this->voxelBasedMeasureGradient->ny * + this->voxelBasedMeasureGradient->nz; + 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) + { + region2VBMGPtr[n] = combinedVBMGPtr[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) + { + region1VBMGPtr[n] = combinedVBMGPtr[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) + region1VBMGPtr[n] = combinedVBMGPtr[n]; + else + region2VBMGPtr[n] = combinedVBMGPtr[n]; + } + } + } + + + //convert region 1 voxel-based gradient to CPG gradient + + //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 + ); + // 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 + ); + // 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 + ); + } + //now resample voxel-based gradient at control points to get transformationGradient + 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 + ); + + + //convert region 2 voxel=based gradient to CPG gradient + + //first convolve voxel-based gardient with a spline kernel + // Convolution along the x axis + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->region2ControlPointGrid->dx; + activeAxis[0] = 1; + activeAxis[2] = 0; + 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->region2ControlPointGrid->dy; + activeAxis[0] = 0; + activeAxis[1] = 1; + 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->region2VoxelBasedMeasureGradientImage->nz > 1) + { + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->region2ControlPointGrid->dz; + activeAxis[1] = 0; + activeAxis[2] = 1; + reg_tools_kernelConvolution(this->region2VoxelBasedMeasureGradientImage, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis + ); + } + //now resample voxel-based gradient at control points to get transformationGradient + if (this->currentFloating->sform_code>0) + reorientation = this->currentFloating->sto_ijk; + else reorientation = this->currentFloating->qto_ijk; + reg_voxelCentric2NodeCentric(this->region2TransformationGradient, + this->region2VoxelBasedMeasureGradientImage, + this->similarityWeight, + false, // no update + &reorientation + ); + + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetSimilarityMeasureGradient"); +#endif } /* *************************************************************** */ /* *************************************************************** */ diff --git a/reg-lib/_reg_f3d_sli.h b/reg-lib/_reg_f3d_sli.h index 9254708b..f5761a73 100644 --- a/reg-lib/_reg_f3d_sli.h +++ b/reg-lib/_reg_f3d_sli.h @@ -51,11 +51,17 @@ class reg_f3d_sli : public reg_f3d double bestWGO; - //reimplement function to get deformation field + //reimplement method to get deformation field //combines deformation fields from each region based on warped distance maps virtual void GetDeformationField(); - //new functions for Gap-Overlap penalty term + //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 methods for Gap-Overlap penalty term //virtual double GetGapOverlapPenaltyTerm(); //virtual void GetGapOverlapGradient(); From b1f29aed80beb79fda624d7ddc1d28562c439353 Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Mon, 28 Aug 2017 11:51:59 +0100 Subject: [PATCH 03/16] implemented methods for penalty term tidied up GetDeformationField() and GetSimilarityMeasureGradient() methods --- niftyreg_build_version.txt | 2 +- reg-lib/_reg_f3d_sli.cpp | 496 +++++++++++++++++++++++-------------- reg-lib/_reg_f3d_sli.h | 19 +- 3 files changed, 325 insertions(+), 192 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index ea90ee31..9e5feb52 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -45 +46 diff --git a/reg-lib/_reg_f3d_sli.cpp b/reg-lib/_reg_f3d_sli.cpp index 4064b7b5..c900b6da 100644 --- a/reg-lib/_reg_f3d_sli.cpp +++ b/reg-lib/_reg_f3d_sli.cpp @@ -85,16 +85,18 @@ reg_f3d_sli::~reg_f3d_sli() 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, - this->currentMask, + NULL, //no mask false, //composition true); //bspline //get deformation field for region2 reg_spline_getDeformationField(this->region2ControlPointGrid, this->region2DeformationFieldImage, - this->currentMask, + NULL, //no mask false, //composition true); //bspline @@ -102,29 +104,33 @@ void reg_f3d_sli::GetDeformationField() reg_resampleImage(this->currentDistanceMap, this->warpedDistanceMapRegion1, this->region1DeformationFieldImage, - this->currentMask, + 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, - this->currentMask, + 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); @@ -132,179 +138,99 @@ void reg_f3d_sli::GetDeformationField() if (this->region1DeformationFieldImage->nz > 1) { //extra pointers required for 3D - T *region1DFPtrZ = ®ion1DFPtrY[numVox]; - T *region2DFPtrZ = ®ion2DFPtrY[numVox]; - T *combinedDFPtrZ = &combinedDFPtrY[numVox]; + region1DFPtrZ = ®ion1DFPtrY[numVox]; + region2DFPtrZ = ®ion2DFPtrY[numVox]; + combinedDFPtrZ = &combinedDFPtrY[numVox]; + } - //loop over voxels - for (size_t n = 0; n < numVox; n++) + //loop over voxels + for (size_t n = 0; n < numVox; n++) + { + //check in mask + if (this->currentMask[n] > -1) { - //check in mask - if (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]) { - //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]) { - 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(); + //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(); - combinedDFPtrZ[n] = std::numeric_limits::quiet_NaN(); - } - else - { - //set combined def field to region2 def field - combinedDFPtrX[n] = region2DFPtrX[n]; - combinedDFPtrY[n] = region2DFPtrY[n]; - 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]) + //check if region2 transform maps into region1, i.e. if region2 WDM < 0 + if (warpedDMR2Ptr[n] < 0) { - //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(); + //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]; - 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]; - combinedDFPtrZ[n] = region1DFPtrZ[n]; - } - else - { - combinedDFPtrX[n] = region2DFPtrX[n]; - combinedDFPtrY[n] = region2DFPtrY[n]; + //set combined def field to region2 def field + 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 + } + } + }//if (warpedDMR1Ptr[n] != warpedDMR1Ptr[n]) else { - combinedDFPtrX[n] = std::numeric_limits::quiet_NaN(); - combinedDFPtrY[n] = std::numeric_limits::quiet_NaN(); - combinedDFPtrZ[n] = std::numeric_limits::quiet_NaN(); - } - }//for (size_t n = 0; n < numVox; n++) - }//if (this->region1DeformationFieldImage->nz > 1) - else - { - //images are 2D - //loop over voxels - for (size_t n = 0; n < numVox; n++) - { - //check in mask - if (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]) + //region1 WDM is not NaN, but still need to check region2 WDM + if (warpedDMR2Ptr[n] != warpedDMR2Ptr[n]) { - 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) { - //both WDMs are NaN so set combined def field to NaN + //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(); - } - else - { - //set combined def field to region2 def field - combinedDFPtrX[n] = region2DFPtrX[n]; - combinedDFPtrY[n] = region2DFPtrY[n]; - } + //set combined def field to region1 def field + combinedDFPtrX[n] = region1DFPtrX[n]; + combinedDFPtrY[n] = region1DFPtrY[n]; + if (combinedDFPtrZ != NULL) + combinedDFPtrZ[n] = region1DFPtrZ[n]; } - }//if (warpedDMR1Ptr[n] != warpedDMR1Ptr[n]) + } else { - //region1 WDM is not NaN, but still need to check region2 WDM - if (warpedDMR2Ptr[n] != warpedDMR2Ptr[n]) + //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) { - //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(); - } - else - { - //set combined def field to region1 def field - combinedDFPtrX[n] = region1DFPtrX[n]; - combinedDFPtrY[n] = region1DFPtrY[n]; - } + 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]; - } - else - { - combinedDFPtrX[n] = region2DFPtrX[n]; - combinedDFPtrY[n] = region2DFPtrY[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 - else - { - combinedDFPtrX[n] = std::numeric_limits::quiet_NaN(); - combinedDFPtrY[n] = std::numeric_limits::quiet_NaN(); - } - }//for (size_t n = 0; n < numVox; n++) - }//else (this->region1DeformationFieldImage->nz > 1) + 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()"); @@ -372,7 +298,7 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() } - //convert region 1 voxel-based gradient to CPG gradient + //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; @@ -385,8 +311,13 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() kernel_type, NULL, // mask NULL, // all volumes are considered as active - activeAxis - ); + 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; @@ -396,8 +327,13 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() kernel_type, NULL, // mask NULL, // all volumes are considered as active - activeAxis - ); + 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) { @@ -409,10 +345,17 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() kernel_type, NULL, // mask NULL, // all volumes are considered as active - activeAxis - ); + activeAxis); + reg_tools_kernelConvolution(this->region2VoxelBasedMeasureGradientImage, + currentNodeSpacing, + kernel_type, + NULL, // mask + NULL, // all volumes are considered as active + activeAxis); } - //now resample voxel-based gradient at control points to get transformationGradient + //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; @@ -421,64 +364,241 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() this->region1VoxelBasedMeasureGradientImage, this->similarityWeight, false, // no update - &reorientation - ); + &reorientation); + reg_voxelCentric2NodeCentric(this->region2TransformationGradient, + this->region2VoxelBasedMeasureGradientImage, + this->similarityWeight, + false, // no update + &reorientation); + +#ifndef NDEBUG + reg_print_fct_debug("reg_f3d_sli::GetSimilarityMeasureGradient()"); +#endif +} +/* *************************************************************** */ +/* *************************************************************** */ +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. - //convert region 2 voxel=based gradient to CPG gradient + //NOTE2 - the gap-overlap penalty term is calculate 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 - //first convolve voxel-based gardient with a spline kernel + 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 +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 - currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->region2ControlPointGrid->dx; - activeAxis[0] = 1; - activeAxis[2] = 0; - reg_tools_kernelConvolution(this->region2VoxelBasedMeasureGradientImage, + 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 - ); + activeAxis); // Convolution along the y axis - currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->region2ControlPointGrid->dy; + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->controlPointGrid->dy; activeAxis[0] = 0; activeAxis[1] = 1; - reg_tools_kernelConvolution(this->region2VoxelBasedMeasureGradientImage, + reg_tools_kernelConvolution(this->gapOverlapGradientWRTDefFieldRegion1, currentNodeSpacing, kernel_type, NULL, // mask NULL, // all volumes are considered as active - activeAxis - ); + 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->region2VoxelBasedMeasureGradientImage->nz > 1) + if (this->gapOverlapGradientWRTDefFieldRegion1->nz > 1) { - currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->region2ControlPointGrid->dz; + currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = this->controlPointGrid->dz; activeAxis[1] = 0; activeAxis[2] = 1; - reg_tools_kernelConvolution(this->region2VoxelBasedMeasureGradientImage, + 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 - ); + activeAxis); } - //now resample voxel-based gradient at control points to get transformationGradient - if (this->currentFloating->sform_code>0) - reorientation = this->currentFloating->sto_ijk; - else reorientation = this->currentFloating->qto_ijk; + + //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->region2VoxelBasedMeasureGradientImage, - this->similarityWeight, - false, // no update - &reorientation - ); + this->gapOverlapGradientWRTDefFieldRegion2, + this->gapOverlapWeight, + true, // update the transformation gradient + &reorientation); - -#ifndef NDEBUG - reg_print_fct_debug("reg_f3d_sli::GetSimilarityMeasureGradient"); -#endif } /* *************************************************************** */ /* *************************************************************** */ diff --git a/reg-lib/_reg_f3d_sli.h b/reg-lib/_reg_f3d_sli.h index f5761a73..52ea5570 100644 --- a/reg-lib/_reg_f3d_sli.h +++ b/reg-lib/_reg_f3d_sli.h @@ -42,14 +42,26 @@ class reg_f3d_sli : public reg_f3d nifti_image *distanceMapImage; 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 penalty term + //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 @@ -61,9 +73,10 @@ class reg_f3d_sli : public reg_f3d //gradients virtual void GetSimilarityMeasureGradient(); + //new methods for Gap-Overlap penalty term - //virtual double GetGapOverlapPenaltyTerm(); - //virtual void GetGapOverlapGradient(); + virtual double ComputeGapOverlapPenaltyTerm(); + virtual void GetGapOverlapGradient(); public: reg_f3d_sli(int refTimePoint, int floTimePoint); From 8f528e540f46fa0e26c0d0b5cab80ff9ae7a2b64 Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Sun, 3 Sep 2017 23:16:16 +0100 Subject: [PATCH 04/16] added several more methods to reg_f3d_sli, and rearranged order of some already implemented changed reg_base so that methods for allocating warped, def field, and warped gradent are called after calling method to initialise current level as this method assigns currentDistanceMap which is required to allocate the warped distance maps which are allocated in AllocateWarped method --- niftyreg_build_version.txt | 2 +- reg-lib/_reg_base.cpp | 12 +- reg-lib/_reg_f3d_sli.cpp | 732 ++++++++++++++++++++++++++++++++++--- reg-lib/_reg_f3d_sli.h | 96 ++++- 4 files changed, 779 insertions(+), 63 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 9e5feb52..abac1ea7 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -46 +47 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_sli.cpp b/reg-lib/_reg_f3d_sli.cpp index c900b6da..f8b7fb7f 100644 --- a/reg-lib/_reg_f3d_sli.cpp +++ b/reg-lib/_reg_f3d_sli.cpp @@ -238,6 +238,287 @@ void reg_f3d_sli::GetDeformationField() } /* *************************************************************** */ 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->region1DeformationFieldImage = 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 @@ -379,57 +660,6 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() /* *************************************************************** */ /* *************************************************************** */ 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 calculate 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 void reg_f3d_sli::GetGapOverlapGradient() { //NOTE - this method assumes the deformation fields and the WDMs for each region @@ -601,6 +831,406 @@ void reg_f3d_sli::GetGapOverlapGradient() } /* *************************************************************** */ +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 backward 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 class reg_f3d_sli; #endif diff --git a/reg-lib/_reg_f3d_sli.h b/reg-lib/_reg_f3d_sli.h index 52ea5570..29e752ba 100644 --- a/reg-lib/_reg_f3d_sli.h +++ b/reg-lib/_reg_f3d_sli.h @@ -65,20 +65,106 @@ class reg_f3d_sli : public reg_f3d //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 methods for Gap-Overlap penalty term - virtual double ComputeGapOverlapPenaltyTerm(); + //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(); + + + //virtual void SetOptimiser(); public: + //constructor and destructor methods reg_f3d_sli(int refTimePoint, int floTimePoint); ~reg_f3d_sli(); }; From ef8df3be83ac29c89eb065a7e2fbea6dea52962e Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Mon, 18 Sep 2017 09:18:29 +0100 Subject: [PATCH 05/16] work in progress --- niftyreg_build_version.txt | 2 +- reg-lib/_reg_f3d_sli.cpp | 328 ++++++++++++++++++++++++++++++++++++- reg-lib/_reg_f3d_sli.h | 48 +++++- 3 files changed, 369 insertions(+), 9 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index abac1ea7..21e72e8a 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -47 +48 diff --git a/reg-lib/_reg_f3d_sli.cpp b/reg-lib/_reg_f3d_sli.cpp index f8b7fb7f..8e843e78 100644 --- a/reg-lib/_reg_f3d_sli.cpp +++ b/reg-lib/_reg_f3d_sli.cpp @@ -21,8 +21,9 @@ template reg_f3d_sli::reg_f3d_sli(int refTimePoint, int floTimePoint) :reg_f3d::reg_f3d(refTimePoint, floTimePoint) { - this->executableName = (char *)"NiftyReg F3D SLI"; + this->executableName = (char *)"NiftyReg F3D Sliding regions"; + this->inputRegion2ControlPointGrid = NULL; this->region2ControlPointGrid = NULL; this->region2DeformationFieldImage = NULL; this->region2VoxelBasedMeasureGradientImage = NULL; @@ -31,13 +32,21 @@ reg_f3d_sli::reg_f3d_sli(int refTimePoint, int floTimePoint) this->region1DeformationFieldImage = NULL; this->region1VoxelBasedMeasureGradientImage = NULL; - this->distanceMapImage = 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"); @@ -962,7 +971,7 @@ T reg_f3d_sli::NormaliseGradient() { *r1Ptr++ /= maxGradValue; } - // The backward gradient is normalised + // The region 2 gradient is normalised T *r2Ptr = static_cast(this->region2TransformationGradient->data); for (size_t i = 0; iregion2TransformationGradient->nvox; ++i) { @@ -1211,7 +1220,6 @@ T reg_f3d_sli::InitialiseCurrentLevel() return maxStepSize; } /* *************************************************************** */ -/* *************************************************************** */ template void reg_f3d_sli::ClearCurrentInputImage() { @@ -1227,9 +1235,321 @@ void reg_f3d_sli::ClearCurrentInputImage() } /* *************************************************************** */ /* *************************************************************** */ +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 + // + //if no input grid provided just copy region 1 CPG + if () +} /* *************************************************************** */ /* *************************************************************** */ template class reg_f3d_sli; diff --git a/reg-lib/_reg_f3d_sli.h b/reg-lib/_reg_f3d_sli.h index 29e752ba..792bf032 100644 --- a/reg-lib/_reg_f3d_sli.h +++ b/reg-lib/_reg_f3d_sli.h @@ -29,6 +29,7 @@ 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; @@ -39,7 +40,7 @@ class reg_f3d_sli : public reg_f3d nifti_image *region1VoxelBasedMeasureGradientImage; //variables for distance map image - nifti_image *distanceMapImage; + nifti_image *inputDistanceMap; //pointer to external nifti_image **distanceMapPyramid; nifti_image *currentDistanceMap; @@ -160,13 +161,52 @@ class reg_f3d_sli : public reg_f3d //image virtual void ClearCurrentInputImage(); - - //virtual void SetOptimiser(); - + //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(); }; #endif From 70c2ed2a785b07fd0761dae0ad10597e4ab7c2e6 Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Sun, 8 Oct 2017 23:00:03 +0100 Subject: [PATCH 06/16] finished implementing reg_f3d_sli class modified reg_f3d (executable) to parse and set sliding regions options ready to test --- niftyreg_build_version.txt | 2 +- reg-apps/reg_f3d.cpp | 80 ++++++++++++++++++++++++++++++++++++-- reg-lib/_reg_f3d_sli.cpp | 55 +++++++++++++++++++++++++- 3 files changed, 130 insertions(+), 7 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 21e72e8a..95f9650f 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -48 +49 diff --git a/reg-apps/reg_f3d.cpp b/reg-apps/reg_f3d.cpp index f14894db..4f32db5c 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-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) @@ -792,7 +863,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]); diff --git a/reg-lib/_reg_f3d_sli.cpp b/reg-lib/_reg_f3d_sli.cpp index 8e843e78..85465612 100644 --- a/reg-lib/_reg_f3d_sli.cpp +++ b/reg-lib/_reg_f3d_sli.cpp @@ -1547,8 +1547,59 @@ void reg_f3d_sli::Initialise() //initialise control point grid for region 2 // - //if no input grid provided just copy region 1 CPG - if () + //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); + } + 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); + } + +#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 + } /* *************************************************************** */ /* *************************************************************** */ From ea62eae44e51c24460b6727145d53c07d4c587aa Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Mon, 9 Oct 2017 15:13:13 +0100 Subject: [PATCH 07/16] had forgotten to modify reg_f3d (executable) so that saves CPP from both regions now done --- niftyreg_build_version.txt | 2 +- reg-apps/reg_f3d.cpp | 147 ++++++++++++++++++++++++++++--------- 2 files changed, 113 insertions(+), 36 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 95f9650f..e373ee69 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -49 +50 diff --git a/reg-apps/reg_f3d.cpp b/reg-apps/reg_f3d.cpp index 4f32db5c..f16bd4c7 100755 --- a/reg-apps/reg_f3d.cpp +++ b/reg-apps/reg_f3d.cpp @@ -898,44 +898,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"); + } + if (region1Name.find(".nii") != std::string::npos) + { + region1Name.replace(region1Name.find(".nii"), 4, "_region1.nii"); + region2Name.replace(region2Name.find(".nii"), 4, "_region2.nii"); + } + if (region1Name.find(".hdr") != std::string::npos) + { + region1Name.replace(region1Name.find(".hdr"), 4, "_region1.hdr"); + region2Name.replace(region2Name.find(".hdr"), 4, "_region2.hdr"); + } + 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"); + } + if (region1Name.find(".img") != std::string::npos) + { + region1Name.replace(region1Name.find(".img"), 4, "_region1.img"); + region2Name.replace(region2Name.find(".img"), 4, "_region2.img"); + } + if (region1Name.find(".png") != std::string::npos) + { + region1Name.replace(region1Name.find(".png"), 4, "_region1.png"); + region2Name.replace(region2Name.find(".png"), 4, "_region2.png"); + } + 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) @@ -998,6 +1073,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) From c6a78144b67706723af9da6336f93e1b0f80cc70 Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Tue, 10 Oct 2017 17:10:22 +0100 Subject: [PATCH 08/16] fixed copy-paste error made GetWarpedImage method work with sliding registration results code now saves a result, but it was 'blank', so still more bugs to find and fix! --- niftyreg_build_version.txt | 2 +- reg-lib/_reg_f3d.cpp | 8 ++++---- reg-lib/_reg_f3d_sli.cpp | 33 ++++++++++++++++++++++++++++----- reg-lib/_reg_f3d_sli.h | 4 ++++ 4 files changed, 37 insertions(+), 10 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index e373ee69..82cced27 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -50 +51 diff --git a/reg-lib/_reg_f3d.cpp b/reg-lib/_reg_f3d.cpp index 2a9356d6..6a4e7ac2 100644 --- a/reg-lib/_reg_f3d.cpp +++ b/reg-lib/_reg_f3d.cpp @@ -1021,10 +1021,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(); + AllocateWarped(); + AllocateDeformationField(); + WarpFloatingImage(3); // cubic spline interpolation + 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 index 85465612..0850d659 100644 --- a/reg-lib/_reg_f3d_sli.cpp +++ b/reg-lib/_reg_f3d_sli.cpp @@ -155,8 +155,8 @@ void reg_f3d_sli::GetDeformationField() //loop over voxels for (size_t n = 0; n < numVox; n++) { - //check in mask - if (this->currentMask[n] > -1) + //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 @@ -259,7 +259,7 @@ void reg_f3d_sli::AllocateDeformationField() this->region1DeformationFieldImage = nifti_copy_nim_info(this->deformationFieldImage); this->region1DeformationFieldImage->data = (void *)calloc(this->region1DeformationFieldImage->nvox, this->region1DeformationFieldImage->nbyper); - this->region1DeformationFieldImage = nifti_copy_nim_info(this->deformationFieldImage); + this->region2DeformationFieldImage = nifti_copy_nim_info(this->deformationFieldImage); this->region2DeformationFieldImage->data = (void *)calloc(this->region2DeformationFieldImage->nvox, this->region2DeformationFieldImage->nbyper); @@ -278,12 +278,12 @@ void reg_f3d_sli::ClearDeformationField() if (this->region1DeformationFieldImage != NULL) { nifti_image_free(this->region1DeformationFieldImage); - this->region1DeformationFieldImage == NULL; + this->region1DeformationFieldImage = NULL; } if (this->region2DeformationFieldImage != NULL) { nifti_image_free(this->region2DeformationFieldImage); - this->region2DeformationFieldImage == NULL; + this->region2DeformationFieldImage = NULL; } #ifndef NDEBUG @@ -1600,6 +1600,29 @@ void reg_f3d_sli::Initialise() 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; } /* *************************************************************** */ /* *************************************************************** */ diff --git a/reg-lib/_reg_f3d_sli.h b/reg-lib/_reg_f3d_sli.h index 792bf032..65e955cf 100644 --- a/reg-lib/_reg_f3d_sli.h +++ b/reg-lib/_reg_f3d_sli.h @@ -207,6 +207,10 @@ class reg_f3d_sli : public reg_f3d //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 From 27d2ada99135c8dc26dbf6cffe2ddb8c00dc55b4 Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Mon, 20 Nov 2017 17:29:34 +0000 Subject: [PATCH 09/16] corrected tests for different file extensions when appending 'region1' and 'region2' to name of saved CPGs --- niftyreg_build_version.txt | 2 +- reg-apps/reg_f3d.cpp | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 82cced27..0691f67b 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -51 +52 diff --git a/reg-apps/reg_f3d.cpp b/reg-apps/reg_f3d.cpp index f16bd4c7..58d460a6 100755 --- a/reg-apps/reg_f3d.cpp +++ b/reg-apps/reg_f3d.cpp @@ -953,32 +953,32 @@ int main(int argc, char **argv) region1Name.replace(region1Name.find(".nii.gz"), 7, "_region1.nii.gz"); region2Name.replace(region2Name.find(".nii.gz"), 7, "_region2.nii.gz"); } - if (region1Name.find(".nii") != std::string::npos) + else if (region1Name.find(".nii") != std::string::npos) { region1Name.replace(region1Name.find(".nii"), 4, "_region1.nii"); region2Name.replace(region2Name.find(".nii"), 4, "_region2.nii"); } - if (region1Name.find(".hdr") != std::string::npos) + else if (region1Name.find(".hdr") != std::string::npos) { region1Name.replace(region1Name.find(".hdr"), 4, "_region1.hdr"); region2Name.replace(region2Name.find(".hdr"), 4, "_region2.hdr"); } - if (region1Name.find(".img.gz") != std::string::npos) + 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"); } - if (region1Name.find(".img") != std::string::npos) + else if (region1Name.find(".img") != std::string::npos) { region1Name.replace(region1Name.find(".img"), 4, "_region1.img"); region2Name.replace(region2Name.find(".img"), 4, "_region2.img"); } - if (region1Name.find(".png") != std::string::npos) + else if (region1Name.find(".png") != std::string::npos) { region1Name.replace(region1Name.find(".png"), 4, "_region1.png"); region2Name.replace(region2Name.find(".png"), 4, "_region2.png"); } - if (region1Name.find(".nrrd") != std::string::npos) + else if (region1Name.find(".nrrd") != std::string::npos) { region1Name.replace(region1Name.find(".nrrd"), 5, "_region1.nrrd"); region2Name.replace(region2Name.find(".nrrd"), 5, "_region2.nrrd"); From 0f18902c2a67ee9928ecf9c07aab2cbf67ee2886 Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Mon, 20 Nov 2017 17:40:00 +0000 Subject: [PATCH 10/16] corrected GetDeformationField() method in reg_f3d_sli to work when no mask supplied corrected GetSimilarityMeasureGradient() method to work with vector field previously gap-overlap constraint term value was normalised by number of voxels, but gradient was not now distance map image is pre-normalised during initialisation so that both the gap-overlap constraint term and its gradient are normalised note - weighted distance map is divided by square root of number of voxels, as the two warped distance maps are mutliplied together this will normalise the gap-overlap constraint term by the number of voxels also added some debugging comments which are commented out but still in code --- niftyreg_build_version.txt | 2 +- reg-lib/_reg_f3d.cpp | 13 ++++++ reg-lib/_reg_f3d_sli.cpp | 90 ++++++++++++++++++++++++++++++++------ 3 files changed, 91 insertions(+), 14 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 0691f67b..59343b09 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -52 +53 diff --git a/reg-lib/_reg_f3d.cpp b/reg-lib/_reg_f3d.cpp index 6a4e7ac2..5ccd85fe 100644 --- a/reg-lib/_reg_f3d.cpp +++ b/reg-lib/_reg_f3d.cpp @@ -541,6 +541,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 @@ -591,6 +598,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 diff --git a/reg-lib/_reg_f3d_sli.cpp b/reg-lib/_reg_f3d_sli.cpp index 0850d659..c9dc28e5 100644 --- a/reg-lib/_reg_f3d_sli.cpp +++ b/reg-lib/_reg_f3d_sli.cpp @@ -156,7 +156,7 @@ void reg_f3d_sli::GetDeformationField() for (size_t n = 0; n < numVox; n++) { //check mask exists and in mask - if (this->currentMask != NULL && this->currentMask[n] > -1) + 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 @@ -414,7 +414,7 @@ double reg_f3d_sli::ComputeGapOverlapPenaltyTerm() } //normalise by the number of voxels and return weighted value - gapOverlapTotal /= double(numVox); + //gapOverlapTotal /= double(numVox); return double(this->gapOverlapWeight) * gapOverlapTotal; #ifndef NDEBUG @@ -545,17 +545,29 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() T *warpedDMR1Ptr = static_cast(this->warpedDistanceMapRegion1->data); T *warpedDMR2Ptr = static_cast(this->warpedDistanceMapRegion2->data); //pointers to voxel-based similarity gradients - T *combinedVBMGPtr = static_cast(this->voxelBasedMeasureGradient->data); - T *region1VBMGPtr = static_cast(this->region1VoxelBasedMeasureGradientImage->data); - T *region2VBMGPtr = static_cast(this->region2VoxelBasedMeasureGradientImage->data); - + 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 - size_t numVox = this->voxelBasedMeasureGradient->nx * - this->voxelBasedMeasureGradient->ny * - this->voxelBasedMeasureGradient->nz; for (size_t n = 0; n < numVox; n++) { //is in mask? @@ -566,13 +578,23 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() //then copy voxel-based gradient value in to region2VoxelBasedMeasureGradientImage if (warpedDMR1Ptr[n] != warpedDMR1Ptr[n] && warpedDMR2Ptr[n] >= 0) { - region2VBMGPtr[n] = combinedVBMGPtr[n]; + 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) { - region1VBMGPtr[n] = combinedVBMGPtr[n]; + 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 @@ -580,14 +602,38 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() { //if sum of WDMs < 0 assign value to region 1, else assign to region 2 if (warpedDMR1Ptr[n] + warpedDMR2Ptr[n] < 0) - region1VBMGPtr[n] = combinedVBMGPtr[n]; + { + region1VBMGPtrX[n] = combinedVBMGPtrX[n]; + region1VBMGPtrY[n] = combinedVBMGPtrY[n]; + if (combinedVBMGPtrZ != NULL) + { + region1VBMGPtrZ[n] = combinedVBMGPtrZ[n]; + } + } else - region2VBMGPtr[n] = combinedVBMGPtr[n]; + { + 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 @@ -661,6 +707,13 @@ void reg_f3d_sli::GetSimilarityMeasureGradient() 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()"); @@ -1571,12 +1624,23 @@ void reg_f3d_sli::Initialise() //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 From 42d26eebceb061be3f0de272add344424288ef5d Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Mon, 20 Nov 2017 17:40:49 +0000 Subject: [PATCH 11/16] modified resample image methods to work when deformation field contains NaN values --- niftyreg_build_version.txt | 2 +- reg-lib/cpu/_reg_resampling.cpp | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 59343b09..fb1e7bc8 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -53 +54 diff --git a/reg-lib/cpu/_reg_resampling.cpp b/reg-lib/cpu/_reg_resampling.cpp index 26de09eb..1840a63d 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]); From 8cdfb87d2c7a655309cf7f5a7fee4d819cd527c6 Mon Sep 17 00:00:00 2001 From: Jamie McClelland Date: Fri, 24 Nov 2017 11:58:01 +0000 Subject: [PATCH 12/16] added missing minus sign when calculating the gap-overlap gradient sliding region registrations now appear to be working, but sensible default gap-overlap constraint term weight still needs to be determined --- niftyreg_build_version.txt | 2 +- reg-lib/_reg_f3d_sli.cpp | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index c3f407c0..f6b91e0e 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -55 +56 diff --git a/reg-lib/_reg_f3d_sli.cpp b/reg-lib/_reg_f3d_sli.cpp index c9dc28e5..470f2e56 100644 --- a/reg-lib/_reg_f3d_sli.cpp +++ b/reg-lib/_reg_f3d_sli.cpp @@ -803,16 +803,16 @@ void reg_f3d_sli::GetGapOverlapGradient() if (warpedDMR1Ptr[n] * warpedDMR2Ptr[n] < 0) { //dGO / dDF1 = -WDM2*(dWDM1 / dDF1) - gapOverlapGradR1PtrX[n] = warpedDMR2Ptr[n] * warpedDMGradR1PtrX[n]; - gapOverlapGradR1PtrY[n] = warpedDMR2Ptr[n] * warpedDMGradR1PtrY[n]; + 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]; + 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]; + 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++) From 97de7c7482802b35d330f8669021c6de69366867 Mon Sep 17 00:00:00 2001 From: Bjoern Eiben Date: Thu, 27 Sep 2018 18:03:13 +0100 Subject: [PATCH 13/16] Added reference for sliding registration to readme and usage message. --- README.txt | 13 +++++++++++-- niftyreg_build_version.txt | 2 +- reg-apps/reg_f3d.cpp | 3 ++- 3 files changed, 14 insertions(+), 4 deletions(-) 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 f6b91e0e..e1617e84 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -56 +57 diff --git a/reg-apps/reg_f3d.cpp b/reg-apps/reg_f3d.cpp index d38403ed..bb56348c 100755 --- a/reg-apps/reg_f3d.cpp +++ b/reg-apps/reg_f3d.cpp @@ -134,7 +134,8 @@ void Usage(char *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-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"); From 14da6fa9ad4bf4aaddeaad8c1e10eb69488151b9 Mon Sep 17 00:00:00 2001 From: BBillot Date: Thu, 14 Mar 2019 16:01:09 +0000 Subject: [PATCH 14/16] masking different time points with nans --- reg-lib/cpu/_reg_lncc.cpp | 4 ++-- reg-lib/cpu/_reg_mind.cpp | 16 ++++++++-------- reg-lib/cpu/_reg_tools.cpp | 23 +++++++++++++---------- reg-lib/cpu/_reg_tools.h | 2 +- 4 files changed, 24 insertions(+), 21 deletions(-) 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..35ac5a0d 100644 --- a/reg-lib/cpu/_reg_mind.cpp +++ b/reg-lib/cpu/_reg_mind.cpp @@ -512,8 +512,8 @@ double reg_mind::GetSimilarityMeasureValue() referenceImagePointer->ny * 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_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 From 3de63278ab84fd74f29e4a50526faca6729f8280 Mon Sep 17 00:00:00 2001 From: Bjoern Eiben Date: Tue, 30 Apr 2019 11:58:01 +0100 Subject: [PATCH 15/16] Sliding: Fixed compiliation error on Linux (due to polymorphism). --- niftyreg_build_version.txt | 2 +- reg-lib/_reg_f3d.cpp | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index a8fa06e1..4b9026d8 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -62 +63 diff --git a/reg-lib/_reg_f3d.cpp b/reg-lib/_reg_f3d.cpp index 59f47016..29447c1e 100644 --- a/reg-lib/_reg_f3d.cpp +++ b/reg-lib/_reg_f3d.cpp @@ -1036,10 +1036,10 @@ nifti_image **reg_f3d::GetWarpedImage() this->currentFloating = this->inputFloating; this->currentMask=NULL; - AllocateWarped(); - AllocateDeformationField(); - WarpFloatingImage(3); // cubic spline interpolation - 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); From cea262bdffa20fa8efc1bea73232b16acb297331 Mon Sep 17 00:00:00 2001 From: Bjoern Eiben Date: Mon, 24 Aug 2020 13:55:26 +0100 Subject: [PATCH 16/16] Added missing openMP declarations. --- niftyreg_build_version.txt | 2 +- reg-lib/cpu/_reg_mind.cpp | 4 ++-- reg-lib/cpu/_reg_ssd.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index b5489e5e..39f5b693 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -69 +71 diff --git a/reg-lib/cpu/_reg_mind.cpp b/reg-lib/cpu/_reg_mind.cpp index 35ac5a0d..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;voxelIndex