Commit e4eab3bf authored by Marc Modat's avatar Marc Modat

Error in the LNCC computation when a mask is used (it was due to the density...

Error in the LNCC computation when a mask is used (it was due to the density function of the convolution function). I also corrected the gradient of the log space resampling for DTI images
parent a3162e13
......@@ -548,7 +548,7 @@ int main(int argc, char **argv)
fprintf(stderr, "Error when reading the floating mask image: %s\n",argv[i-1]);
return 1;
}
REG->SetReferenceMask(floatingMaskImage);
REG->SetFloatingMask(floatingMaskImage);
}
else if(strcmp(argv[i], "-ic")==0 || strcmp(argv[i], "--ic")==0){
REG->SetInverseConsistencyWeight(atof(argv[++i]));
......
......@@ -51,7 +51,7 @@ typedef struct{
void PetitUsage(char *exec)
{
fprintf(stderr,"Usage:\t%s -target <referenceImageName> -source <floatingImageName> [OPTIONS].\n",exec);
fprintf(stderr,"Usage:\t%s -ref <referenceImageName> -flo <floatingImageName> [OPTIONS].\n",exec);
fprintf(stderr,"\tSee the help for more details (-h).\n");
return;
}
......@@ -72,13 +72,12 @@ void Usage(char *exec)
printf("\t-cpp <filename>\t\tFilename of the control point grid image (from reg_f3d)\n");
printf("\t-def <filename>\t\tFilename of the deformation field image (from reg_transform)\n");
printf("\t*\tThere are no limit for the required output number from the following\n");
printf("\t-res <filename> \tFilename of the resampled image [none]\n");
printf("\t-blank <filename> \tFilename of the resampled blank grid [none]\n");
printf("\t*\tOthers\n");
printf("\t-NN \t\t\tUse a Nearest Neighbor interpolation for the source resampling (cubic spline by default)\n");
printf("\t-LIN \t\t\tUse a Linear interpolation for the source resampling (cubic spline by default)\n");
printf("\t*\tOthers\n");
printf("\t-inter <int> \t\tInterpolation order (0,1,3)[3]\n");
printf("\t-pad <int> \t\tInterpolation padding value [0]\n");
printf("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n");
return;
}
......@@ -211,19 +210,19 @@ int main(int argc, char **argv)
return 1;
}
/* Read the target image */
/* Read the reference image */
nifti_image *referenceImage = reg_io_ReadImageHeader(param->referenceImageName);
if(referenceImage == NULL){
fprintf(stderr,"[NiftyReg ERROR] Error when reading the target image: %s\n",
fprintf(stderr,"[NiftyReg ERROR] Error when reading the reference image: %s\n",
param->referenceImageName);
return 1;
}
reg_checkAndCorrectDimension(referenceImage);
/* Read the source image */
/* Read the floating image */
nifti_image *floatingImage = reg_io_ReadImageFile(param->floatingImageName);
if(floatingImage == NULL){
fprintf(stderr,"[NiftyReg ERROR] Error when reading the source image: %s\n",
fprintf(stderr,"[NiftyReg ERROR] Error when reading the floating image: %s\n",
param->floatingImageName);
return 1;
}
......@@ -243,10 +242,10 @@ int main(int argc, char **argv)
for(int i=0;i<argc;i++) printf(" %s", argv[i]);
printf("\n\n");
printf("Parameters\n");
printf("Target image name: %s\n",referenceImage->fname);
printf("Reference image name: %s\n",referenceImage->fname);
printf("\t%ix%ix%i voxels, %i volumes\n",referenceImage->nx,referenceImage->ny,referenceImage->nz,referenceImage->nt);
printf("\t%gx%gx%g mm\n",referenceImage->dx,referenceImage->dy,referenceImage->dz);
printf("Source image name: %s\n",floatingImage->fname);
printf("Floating image name: %s\n",floatingImage->fname);
printf("\t%ix%ix%i voxels, %i volumes\n",floatingImage->nx,floatingImage->ny,floatingImage->nz,floatingImage->nt);
printf("\t%gx%gx%g mm\n",floatingImage->dx,floatingImage->dy,floatingImage->dz);
printf("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n\n");
......@@ -372,7 +371,7 @@ int main(int argc, char **argv)
progressXML(2, "Deformation field ready...");
/* ************************* */
/* RESAMPLE THE SOURCE IMAGE */
/* WARP THE FLOATING IMAGE */
/* ************************* */
if(flag->outputResultFlag){
switch(param->interpolation){
......@@ -386,21 +385,21 @@ int main(int argc, char **argv)
param->interpolation=3;
break;
}
nifti_image *resultImage = nifti_copy_nim_info(referenceImage);
resultImage->dim[0]=resultImage->ndim=floatingImage->dim[0];
resultImage->dim[4]=resultImage->nt=floatingImage->dim[4];
resultImage->cal_min=floatingImage->cal_min;
resultImage->cal_max=floatingImage->cal_max;
resultImage->scl_slope=floatingImage->scl_slope;
resultImage->scl_inter=floatingImage->scl_inter;
resultImage->datatype = floatingImage->datatype;
resultImage->nbyper = floatingImage->nbyper;
resultImage->nvox = (size_t)resultImage->dim[1] * (size_t)resultImage->dim[2] *
(size_t)resultImage->dim[3] * (size_t)resultImage->dim[4];
resultImage->data = (void *)calloc(resultImage->nvox, resultImage->nbyper);
nifti_image *warpedImage = nifti_copy_nim_info(referenceImage);
warpedImage->dim[0]=warpedImage->ndim=floatingImage->dim[0];
warpedImage->dim[4]=warpedImage->nt=floatingImage->dim[4];
warpedImage->cal_min=floatingImage->cal_min;
warpedImage->cal_max=floatingImage->cal_max;
warpedImage->scl_slope=floatingImage->scl_slope;
warpedImage->scl_inter=floatingImage->scl_inter;
warpedImage->datatype = floatingImage->datatype;
warpedImage->nbyper = floatingImage->nbyper;
warpedImage->nvox = (size_t)warpedImage->dim[1] * (size_t)warpedImage->dim[2] *
(size_t)warpedImage->dim[3] * (size_t)warpedImage->dim[4];
warpedImage->data = (void *)calloc(warpedImage->nvox, warpedImage->nbyper);
if(floatingImage->dim[4]==6 || floatingImage->dim[4]==7){
#ifndef _NDEBUG
#ifndef NDEBUG
printf("[NiftyReg DEBUG] DTI-based resampling\n");
#endif
// Compute first the Jacobian matrices
......@@ -414,7 +413,7 @@ int main(int argc, char **argv)
bool timepoints[7]; for(int i=0;i<7;++i) timepoints[i]=true;
if(floatingImage->dim[4]==7) timepoints[0]=false;
reg_resampleImage(floatingImage,
resultImage,
warpedImage,
deformationFieldImage,
NULL,
param->interpolation,
......@@ -425,17 +424,17 @@ int main(int argc, char **argv)
}
else{
reg_resampleImage(floatingImage,
resultImage,
warpedImage,
deformationFieldImage,
NULL,
param->interpolation,
param->paddingValue);
}
memset(resultImage->descrip, 0, 80);
strcpy (resultImage->descrip,"Warped image using NiftyReg (reg_resample)");
reg_io_WriteImageFile(resultImage,param->outputResultName);
memset(warpedImage->descrip, 0, 80);
strcpy (warpedImage->descrip,"Warped image using NiftyReg (reg_resample)");
reg_io_WriteImageFile(warpedImage,param->outputResultName);
printf("[NiftyReg] Resampled image has been saved: %s\n", param->outputResultName);
nifti_image_free(resultImage);
nifti_image_free(warpedImage);
}
/* *********************** */
......@@ -482,26 +481,26 @@ int main(int argc, char **argv)
}
}
nifti_image *resultImage = nifti_copy_nim_info(referenceImage);
resultImage->dim[0]=resultImage->ndim=3;
resultImage->dim[4]=resultImage->nt=1;
resultImage->cal_min=floatingImage->cal_min;
resultImage->cal_max=floatingImage->cal_max;
resultImage->scl_slope=floatingImage->scl_slope;
resultImage->scl_inter=floatingImage->scl_inter;
resultImage->datatype =NIFTI_TYPE_UINT8;
resultImage->nbyper = sizeof(unsigned char);
resultImage->data = (void *)calloc(resultImage->nvox, resultImage->nbyper);
nifti_image *warpedImage = nifti_copy_nim_info(referenceImage);
warpedImage->dim[0]=warpedImage->ndim=3;
warpedImage->dim[4]=warpedImage->nt=1;
warpedImage->cal_min=floatingImage->cal_min;
warpedImage->cal_max=floatingImage->cal_max;
warpedImage->scl_slope=floatingImage->scl_slope;
warpedImage->scl_inter=floatingImage->scl_inter;
warpedImage->datatype =NIFTI_TYPE_UINT8;
warpedImage->nbyper = sizeof(unsigned char);
warpedImage->data = (void *)calloc(warpedImage->nvox, warpedImage->nbyper);
reg_resampleImage(gridImage,
resultImage,
warpedImage,
deformationFieldImage,
NULL,
1, // linear interpolation
0);
memset(resultImage->descrip, 0, 80);
strcpy (resultImage->descrip,"Warped regular grid using NiftyReg (reg_resample)");
reg_io_WriteImageFile(resultImage,param->outputBlankName);
nifti_image_free(resultImage);
memset(warpedImage->descrip, 0, 80);
strcpy (warpedImage->descrip,"Warped regular grid using NiftyReg (reg_resample)");
reg_io_WriteImageFile(warpedImage,param->outputBlankName);
nifti_image_free(warpedImage);
nifti_image_free(gridImage);
printf("[NiftyReg] Resampled grid has been saved: %s\n", param->outputBlankName);
}
......
......@@ -252,24 +252,24 @@ int main(int argc, char **argv)
bool boolX[3]={1,0,0};
for(int i=0; i<smoothImg->nt*smoothImg->nu;++i) kernelSize[i]=param->smoothValueX;
if(flag->smoothMeanFlag)
reg_tools_kernelConvolution(smoothImg,kernelSize,2,timePoint,boolX);
reg_tools_kernelConvolution(smoothImg,kernelSize,2,NULL,timePoint,boolX);
else if(flag->smoothSplineFlag)
reg_tools_kernelConvolution(smoothImg,kernelSize,1,timePoint,boolX);
else reg_tools_kernelConvolution(smoothImg,kernelSize,0,timePoint,boolX);
reg_tools_kernelConvolution(smoothImg,kernelSize,1,NULL,timePoint,boolX);
else reg_tools_kernelConvolution(smoothImg,kernelSize,0,NULL,timePoint,boolX);
bool boolY[3]={0,1,0};
for(int i=0; i<smoothImg->nt*smoothImg->nu;++i) kernelSize[i]=param->smoothValueY;
if(flag->smoothMeanFlag)
reg_tools_kernelConvolution(smoothImg,kernelSize,2,timePoint,boolY);
reg_tools_kernelConvolution(smoothImg,kernelSize,2,NULL,timePoint,boolY);
else if(flag->smoothSplineFlag)
reg_tools_kernelConvolution(smoothImg,kernelSize,1,timePoint,boolY);
else reg_tools_kernelConvolution(smoothImg,kernelSize,0,timePoint,boolY);
reg_tools_kernelConvolution(smoothImg,kernelSize,1,NULL,timePoint,boolY);
else reg_tools_kernelConvolution(smoothImg,kernelSize,0,NULL,timePoint,boolY);
bool boolZ[3]={0,0,1};
for(int i=0; i<smoothImg->nt*smoothImg->nu;++i) kernelSize[i]=param->smoothValueZ;
if(flag->smoothMeanFlag)
reg_tools_kernelConvolution(smoothImg,kernelSize,2,timePoint,boolZ);
reg_tools_kernelConvolution(smoothImg,kernelSize,2,NULL,timePoint,boolZ);
else if(flag->smoothSplineFlag)
reg_tools_kernelConvolution(smoothImg,kernelSize,1,timePoint,boolZ);
else reg_tools_kernelConvolution(smoothImg,kernelSize,0,timePoint,boolZ);
reg_tools_kernelConvolution(smoothImg,kernelSize,1,NULL,timePoint,boolZ);
else reg_tools_kernelConvolution(smoothImg,kernelSize,0,NULL,timePoint,boolZ);
delete []kernelSize;
delete []timePoint;
if(flag->outputImageFlag)
......
......@@ -42,9 +42,10 @@ foreach(NAME ${LIST})
RUNTIME DESTINATION bin COMPONENT Development
LIBRARY DESTINATION lib COMPONENT Development
ARCHIVE DESTINATION lib COMPONENT Development
)
install(FILES ${NAME}.h DESTINATION include COMPONENT Development)
)
install(FILES ${NAME}.h DESTINATION include COMPONENT Development)
endforeach(NAME)
install(FILES _reg_measure.h DESTINATION include COMPONENT Development)
#-----------------------------------------------------------------------------
set(NAME _reg_optimiser)
if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
......
......@@ -237,7 +237,7 @@ void reg_aladin<T>::InitialiseRegistration()
for(int i=1;i<this->ReferencePyramid[l]->nt;++i)
active[i]=false;
sigma[0]=this->ReferenceSigma;
reg_tools_kernelConvolution(this->ReferencePyramid[l], sigma, 0, active);
reg_tools_kernelConvolution(this->ReferencePyramid[l], sigma, 0, NULL, active);
delete []active;
delete []sigma;
}
......@@ -249,7 +249,7 @@ void reg_aladin<T>::InitialiseRegistration()
for(int i=1;i<this->FloatingPyramid[l]->nt;++i)
active[i]=false;
sigma[0]=this->FloatingSigma;
reg_tools_kernelConvolution(this->FloatingPyramid[l], sigma, 0, active);
reg_tools_kernelConvolution(this->FloatingPyramid[l], sigma, 0, NULL, active);
delete []active;
delete []sigma;
}
......
......@@ -538,9 +538,12 @@ void reg_base<T>::CheckParameters()
// CHECK THE MASK DIMENSION IF IT IS DEFINED
if(this->maskImage!=NULL){
if(this->inputReference->nx != maskImage->nx ||
this->inputReference->ny != maskImage->ny ||
this->inputReference->nz != maskImage->nz){
if(this->inputReference->nx != this->maskImage->nx ||
this->inputReference->ny != this->maskImage->ny ||
this->inputReference->nz != this->maskImage->nz ){
printf("x: %i %i\n",this->inputReference->nx, this->maskImage->nx);
printf("y: %i %i\n",this->inputReference->ny, this->maskImage->ny);
printf("z: %i %i\n",this->inputReference->nz, this->maskImage->nz);
fprintf(stderr,"[NiftyReg ERROR] The mask image has different x, y or z dimension than the reference image.\n");
reg_exit(1);
}
......@@ -690,7 +693,7 @@ void reg_base<T>::Initisalise()
for(int i=1;i<this->referencePyramid[l]->nt;++i)
active[i]=false;
sigma[0]=this->referenceSmoothingSigma;
reg_tools_kernelConvolution(this->referencePyramid[l], sigma, 0, active);
reg_tools_kernelConvolution(this->referencePyramid[l], sigma, 0, NULL, active);
delete []active;
delete []sigma;
}
......@@ -702,7 +705,7 @@ void reg_base<T>::Initisalise()
for(int i=1;i<this->floatingPyramid[l]->nt;++i)
active[i]=false;
sigma[0]=this->floatingSmoothingSigma;
reg_tools_kernelConvolution(this->floatingPyramid[l], sigma, 0, active);
reg_tools_kernelConvolution(this->floatingPyramid[l], sigma, 0, NULL, active);
delete []active;
delete []sigma;
}
......@@ -769,7 +772,8 @@ void reg_base<T>::GetVoxelBasedGradient()
this->interpolation,
this->warpedPaddingValue,
this->measure_dti->GetActiveTimepoints(),
this->forwardJacobianMatrix);
this->forwardJacobianMatrix,
this->warped);
}
else{
reg_getImageGradient(this->currentFloating,
......@@ -790,7 +794,7 @@ void reg_base<T>::GetVoxelBasedGradient()
if(this->measure_multichannel_nmi!=NULL)
this->measure_multichannel_nmi->GetVoxelBasedSimilarityMeasureGradient();
if(this->measure_ssd!=NULL)
if(this->measure_ssd!=NULL)
this->measure_ssd->GetVoxelBasedSimilarityMeasureGradient();
if(this->measure_kld!=NULL)
......@@ -799,31 +803,9 @@ void reg_base<T>::GetVoxelBasedGradient()
if(this->measure_lncc!=NULL)
this->measure_lncc->GetVoxelBasedSimilarityMeasureGradient();
if(this->measure_dti!=NULL)
if(this->measure_dti!=NULL)
this->measure_dti->GetVoxelBasedSimilarityMeasureGradient();
// nifti_image *temp=nifti_copy_nim_info(this->currentReference);
// temp->dim[0]=temp->ndim=3;
// temp->dim[4] = temp->nt=1;
// temp->dim[5] = temp->nu=1;
// temp->nvox=(size_t)temp->nx*temp->ny*temp->nz;
// temp->data=(void *)malloc(temp->nvox*temp->nbyper);
// T *tempPtr = static_cast<T *>(temp->data);
// T *gradPtr = static_cast<T *>(this->warpedGradientImage->data);
// for(int i=0; i<this->warpedGradientImage->nt*this->warpedGradientImage->nu; ++i){
// char name[255];
// sprintf(name, "grad_%i.nii", i);
// memcpy(tempPtr,&gradPtr[i*temp->nvox],temp->nvox*temp->nbyper);
// reg_io_WriteImageFile(temp,name);
// }
// nifti_image_free(temp);
// reg_io_WriteImageFile(this->warpedGradientImage,
// "spatialGradient.nii");
// reg_io_WriteImageFile(this->voxelBasedMeasureGradientImage,
// "measureGradient.nii");
// reg_exit(1);
return;
}
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
......@@ -928,7 +910,7 @@ void reg_base<T>::WarpFloatingImage(int inter)
// Compute the deformation field
this->GetDeformationField();
if(measure_dti==NULL){
if(this->measure_dti==NULL){
// Resample the floating image
reg_resampleImage(this->currentFloating,
this->warped,
......
......@@ -37,6 +37,8 @@ reg_f3d_sym<T>::reg_f3d_sym(int refTimePoint,int floTimePoint)
this->floatingMaskPyramid=NULL;
this->backwardActiveVoxelNumber=NULL;
this->backwardJacobianMatrix=NULL;
this->inverseConsistencyWeight=0.1;
#ifndef NDEBUG
......@@ -157,11 +159,11 @@ void reg_f3d_sym<T>::AllocateWarped()
template <class T>
void reg_f3d_sym<T>::ClearWarped()
{
reg_f3d<T>::ClearWarped();
if(this->backwardWarped!=NULL){
nifti_image_free(this->backwardWarped);
this->backwardWarped=NULL;
}
reg_f3d<T>::ClearWarped();
if(this->backwardWarped!=NULL){
nifti_image_free(this->backwardWarped);
this->backwardWarped=NULL;
}
return;
}
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
......@@ -210,6 +212,13 @@ void reg_f3d_sym<T>::AllocateDeformationField()
strcpy(this->backwardDeformationFieldImage->intent_name,"NREG_TRANS");
this->backwardDeformationFieldImage->intent_p1=DEF_FIELD;
if(this->measure_dti!=NULL)
this->backwardJacobianMatrix=(mat33 *)malloc(
this->backwardDeformationFieldImage->nx *
this->backwardDeformationFieldImage->ny *
this->backwardDeformationFieldImage->nz *
sizeof(mat33));
return;
}
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
......@@ -221,6 +230,10 @@ void reg_f3d_sym<T>::ClearDeformationField()
nifti_image_free(this->backwardDeformationFieldImage);
this->backwardDeformationFieldImage=NULL;
}
if(this->backwardJacobianMatrix!=NULL){
free(this->backwardJacobianMatrix);
this->backwardJacobianMatrix=NULL;
}
return;
}
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
......@@ -493,20 +506,49 @@ void reg_f3d_sym<T>::WarpFloatingImage(int inter)
// Compute the deformation fields
this->GetDeformationField();
// Resample the floating image
reg_resampleImage(this->currentFloating, // input image
this->warped, // warped input image
this->deformationFieldImage, // deformation field
this->currentMask, // mask
inter, // interpolation
this->warpedPaddingValue); // padding value
// Resample the reference image
reg_resampleImage(this->currentReference, // input image
this->backwardWarped, // warped input image
this->backwardDeformationFieldImage, // deformation field
this->currentFloatingMask, // mask
inter, // interpolation type
this->warpedPaddingValue); // padding value
// Resample the floating image
if(this->measure_dti==NULL){
reg_resampleImage(this->currentFloating,
this->warped,
this->deformationFieldImage,
this->currentMask,
inter,
this->warpedPaddingValue);
}
else{
reg_defField_getJacobianMatrix(this->deformationFieldImage,
this->forwardJacobianMatrix);
reg_resampleImage(this->currentFloating,
this->warped,
this->deformationFieldImage,
this->currentMask,
inter,
this->warpedPaddingValue,
this->measure_dti->GetActiveTimepoints(),
this->forwardJacobianMatrix);
}
// Resample the reference image
if(this->measure_dti==NULL){
reg_resampleImage(this->currentReference, // input image
this->backwardWarped, // warped input image
this->backwardDeformationFieldImage, // deformation field
this->currentFloatingMask, // mask
inter, // interpolation type
this->warpedPaddingValue); // padding value
}
else{
reg_defField_getJacobianMatrix(this->backwardDeformationFieldImage,
this->backwardJacobianMatrix);
reg_resampleImage(this->currentReference, // input image
this->backwardWarped, // warped input image
this->backwardDeformationFieldImage, // deformation field
this->currentFloatingMask, // mask
inter, // interpolation type
this->warpedPaddingValue, // padding value
this->measure_dti->GetActiveTimepoints(),
this->backwardJacobianMatrix);
}
return;
}
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
......
......@@ -37,6 +37,8 @@ class reg_f3d_sym : public reg_f3d<T>
double *backwardLogJointHistogram;
double backwardEntropies[4];
mat33 *backwardJacobianMatrix;
T inverseConsistencyWeight;
double currentIC;
double bestIC;
......
......@@ -30,6 +30,10 @@ reg_lncc::reg_lncc()
this->floatingSdevImage=NULL;
this->warpedReferenceMeanImage=NULL;
this->warpedReferenceSdevImage=NULL;
for(int i=0;i<255;++i)
kernelStandardDeviation[i]=-5.f;
#ifndef NDEBUG
printf("[NiftyReg DEBUG] reg_lncc constructor called\n");
#endif
......@@ -75,7 +79,8 @@ reg_lncc::~reg_lncc()
template <class DTYPE>
void reg_lncc::UpdateLocalStatImages(nifti_image *originalImage,
nifti_image *meanImage,
nifti_image *stdDevImage)
nifti_image *stdDevImage,
int *mask)
{
DTYPE *origPtr = static_cast<DTYPE *>(originalImage->data);
DTYPE *meanPtr = static_cast<DTYPE *>(meanImage->data);
......@@ -83,8 +88,8 @@ void reg_lncc::UpdateLocalStatImages(nifti_image *originalImage,
memcpy(meanPtr, origPtr, originalImage->nvox*originalImage->nbyper);
memcpy(sdevPtr, origPtr, originalImage->nvox*originalImage->nbyper);
reg_tools_multiplyImageToImage(stdDevImage, stdDevImage, stdDevImage);
reg_tools_kernelConvolution(meanImage, this->kernelStandardDeviation, KERNEL_TYPE, this->activeTimePoint);
reg_tools_kernelConvolution(stdDevImage, this->kernelStandardDeviation, KERNEL_TYPE, this->activeTimePoint);
reg_tools_kernelConvolution(meanImage, this->kernelStandardDeviation, KERNEL_TYPE, mask, this->activeTimePoint);
reg_tools_kernelConvolution(stdDevImage, this->kernelStandardDeviation, KERNEL_TYPE, mask, this->activeTimePoint);
for(size_t voxel=0;voxel<originalImage->nvox;++voxel){
// G*(I^2) - (G*I)^2
sdevPtr[voxel] = sqrt(sdevPtr[voxel] - reg_pow2(meanPtr[voxel]));
......@@ -185,12 +190,14 @@ void reg_lncc::InitialiseMeasure(nifti_image *refImgPtr,
case NIFTI_TYPE_FLOAT32:
this->UpdateLocalStatImages<float>(this->referenceImagePointer,
this->referenceMeanImage,
this->referenceSdevImage);
this->referenceSdevImage,
this->referenceMaskPointer);
break;
case NIFTI_TYPE_FLOAT64:
this->UpdateLocalStatImages<double>(this->referenceImagePointer,
this->referenceMeanImage,
this->referenceSdevImage);
this->referenceSdevImage,
this->referenceMaskPointer);
break;
}
if(this->isSymmetric){
......@@ -218,12 +225,14 @@ void reg_lncc::InitialiseMeasure(nifti_image *refImgPtr,
case NIFTI_TYPE_FLOAT32:
this->UpdateLocalStatImages<float>(this->floatingImagePointer,
this->floatingMeanImage,
this->floatingSdevImage);
this->floatingSdevImage,
this->floatingMaskPointer);
break;
case NIFTI_TYPE_FLOAT64:
this->UpdateLocalStatImages<double>(this->floatingImagePointer,
this->floatingMeanImage,
this->floatingSdevImage);
this->floatingSdevImage,
this->floatingMaskPointer);
break;
}
}
......@@ -251,7 +260,7 @@ double reg_getLNCCValue(nifti_image *referenceImage,
{
// Compute the local correlation
reg_tools_multiplyImageToImage(referenceImage, warpedImage, correlationImage);
reg_tools_kernelConvolution(correlationImage, kernelStandardDeviation, KERNEL_TYPE, activeTimePoint);
reg_tools_kernelConvolution(correlationImage, kernelStandardDeviation, KERNEL_TYPE, refMask, activeTimePoint);
double lncc_value_sum = 0., lncc_value;
double activeVoxel_num = 0.;
......@@ -263,6 +272,7 @@ double reg_getLNCCValue(nifti_image *referenceImage,
DTYPE *correlaPtr=static_cast<DTYPE *>(correlationImage->data);
size_t voxel, voxelNumber = (size_t)referenceImage->nx*
referenceImage->ny*referenceImage->nz;
// Iteration over all time points
for(int t=0; t<referenceImage->nt; ++t){
if(activeTimePoint[t]==true){
......@@ -289,10 +299,10 @@ double reg_getLNCCValue(nifti_image *referenceImage,
) /
(refSdevPtr0[voxel]*warSdevPtr0[voxel]);
if(lncc_value==lncc_value && isinf(lncc_value)==0){
if(lncc_value==lncc_value && fabs(lncc_value)<1.01 && isinf(lncc_value)==0){
lncc_value_sum += fabs(lncc_value);
++activeVoxel_num;
}
}
}
}
}
......@@ -309,12 +319,14 @@ double reg_lncc::GetSimilarityMeasureValue()
case NIFTI_TYPE_FLOAT32:
this->UpdateLocalStatImages<float>(this->warpedFloatingImagePointer,
this->warpedFloatingMeanImage,
this->warpedFloatingSdevImage);
this->warpedFloatingSdevImage,
this->referenceMaskPointer);
break;
case NIFTI_TYPE_FLOAT64:
this->UpdateLocalStatImages<double>(this->warpedFloatingImagePointer,
this->warpedFloatingMeanImage,
this->warpedFloatingSdevImage);
this->warpedFloatingSdevImage,
this->referenceMaskPointer);
break;
}
// Compute the LNCC - Forward
......@@ -350,12 +362,14 @@ double reg_lncc::GetSimilarityMeasureValue()
case NIFTI_TYPE_FLOAT32:
this->UpdateLocalStatImages<float>(this->warpedReferenceImagePointer,
this->warpedReferenceMeanImage,
this->warpedReferenceSdevImage);
this->warpedReferenceSdevImage,
this->floatingMaskPointer);
break;
case NIFTI_TYPE_FLOAT64:
this->UpdateLocalStatImages<double>(this->warpedReferenceImagePointer,
this->warpedReferenceMeanImage,
this->warpedReferenceSdevImage);
this->warpedReferenceSdevImage,
this->floatingMaskPointer);
break;
}
// Compute the LNCC - Backward
......@@ -406,7 +420,7 @@ void reg_getVoxelBasedLNCCGradient(nifti_image *referenceImage,
{
// Compute the local correlation