Commit 834fecf2 authored by Marc Modat's avatar Marc Modat

Over normalisation of the NMI gradient in the v1.1 (not v1.0) has been fixed -...

Over normalisation of the NMI gradient in the v1.1 (not v1.0) has been fixed - thanks Rolf for the feedback
parent c5978143
......@@ -20,8 +20,8 @@
#define MAX_ITERATIONS 20
#define BLOCK_WIDTH 4
#define BLOCK_2D_SIZE 16
#define BLOCK_SIZE 64
#define BLOCK_2D_SIZE 16
#define OVERLAP_SIZE 3
#define STEP_SIZE 1
#define NUM_BLOCKS_TO_COMPARE 216 // We compare in a 6x6x6 neighborhood.
......
......@@ -2112,7 +2112,7 @@ void reg_bspline_bendingEnergyGradient2D( nifti_image *splineControlPoint,
}
metricGradientValue[0] = (PrecisionTYPE)(*gradientXPtr);
metricGradientValue[1] = (PrecisionTYPE)(*gradientYPtr);
// (Marc) I removed the normalisation by the voxel number as each gradient has to be normalised in the same way
// (Marc) I removed the normalisation by the voxel number as each gradient has to be normalised in the same way (NMI, BE, JAC)
*gradientXPtr++ = (SplineTYPE)((1.0-weight)*metricGradientValue[0] + weight*gradientValue[0]);
*gradientYPtr++ = (SplineTYPE)((1.0-weight)*metricGradientValue[1] + weight*gradientValue[1]);
}
......
......@@ -404,10 +404,9 @@ void reg_getVoxelBasedNMIGradientUsingPW2D( nifti_image *targetImage,
} // t
} // 0<r<bin
} // r
// The gradient is computed in a the voxel space - The target orientation has to be taken into account
PrecisionTYPE temp = (PrecisionTYPE)(entropies[2]);
// (Marc) I removed the normalisation by the voxel number as each gradient has to be normalised in the same way
// (Marc) I removed the normalisation by the voxel number as each gradient has to be normalised in the same way (NMI, BE, JAC)
*nmiGradientPtrX = (NMIGradientTYPE)((fixedEntropyDerivative_X + movingEntropyDerivative_X - NMI * jointEntropyDerivative_X) / temp);
*nmiGradientPtrY = (NMIGradientTYPE)((fixedEntropyDerivative_Y + movingEntropyDerivative_Y - NMI * jointEntropyDerivative_Y) / temp);
......@@ -516,9 +515,8 @@ void reg_getVoxelBasedNMIGradientUsingPW3D( nifti_image *targetImage,
} // 0<r<bin
} // r
// The gradient is computed in a the voxel space - The target orientation has to be taken into account
// (Marc) I removed the normalisation by the voxel number as each gradient has to be normalised in the same way (NMI, BE, JAC)
PrecisionTYPE temp = (PrecisionTYPE)(entropies[2]);
// (Marc) I removed the normalisation by the voxel number as each gradient has to be normalised in the same way
*nmiGradientPtrX = (NMIGradientTYPE)((fixedEntropyDerivative_X + movingEntropyDerivative_X - NMI * jointEntropyDerivative_X) / temp);
*nmiGradientPtrY = (NMIGradientTYPE)((fixedEntropyDerivative_Y + movingEntropyDerivative_Y - NMI * jointEntropyDerivative_Y) / temp);
*nmiGradientPtrZ = (NMIGradientTYPE)((fixedEntropyDerivative_Z + movingEntropyDerivative_Z - NMI * jointEntropyDerivative_Z) / temp);
......
......@@ -161,14 +161,14 @@ void reg_smoothImageForCubicSpline1( nifti_image *image,
/* Smoothing along the X axis */
int windowSize = 2*radius[0] + 1;
PrecisionTYPE *window = (PrecisionTYPE *)calloc(windowSize,sizeof(PrecisionTYPE));
PrecisionTYPE coeffSum=0.0;
// PrecisionTYPE coeffSum=0.0;
for(int it=-radius[0]; it<=radius[0]; it++){
PrecisionTYPE coeff = (PrecisionTYPE)(fabs(2.0*(PrecisionTYPE)it/(PrecisionTYPE)radius[0]));
if(coeff<1.0) window[it+radius[0]] = (PrecisionTYPE)(2.0/3.0 - coeff*coeff + 0.5*coeff*coeff*coeff);
else window[it+radius[0]] = (PrecisionTYPE)(-(coeff-2.0)*(coeff-2.0)*(coeff-2.0)/6.0);
coeffSum += window[it+radius[0]];
// coeffSum += window[it+radius[0]];
}
for(int it=0;it<windowSize;it++) window[it] /= coeffSum;
// for(int it=0;it<windowSize;it++) window[it] /= coeffSum;
for(int t=0;t<timePoint;t++){
for(int u=0;u<field;u++){
......@@ -205,14 +205,14 @@ void reg_smoothImageForCubicSpline1( nifti_image *image,
windowSize = 2*radius[1] + 1;
free(window);
window = (PrecisionTYPE *)calloc(windowSize,sizeof(PrecisionTYPE));
coeffSum=0.0;
// coeffSum=0.0;
for(int it=-radius[1]; it<=radius[1]; it++){
PrecisionTYPE coeff = (PrecisionTYPE)(fabs(2.0*(PrecisionTYPE)it/(PrecisionTYPE)radius[1]));
if(coeff<1.0) window[it+radius[1]] = (PrecisionTYPE)(2.0/3.0 - coeff*coeff + 0.5*coeff*coeff*coeff);
else window[it+radius[1]] = (PrecisionTYPE)(-(coeff-2.0)*(coeff-2.0)*(coeff-2.0)/6.0);
coeffSum += window[it+radius[1]];
// coeffSum += window[it+radius[1]];
}
for(int it=0;it<windowSize;it++)window[it] /= coeffSum;
// for(int it=0;it<windowSize;it++)window[it] /= coeffSum;
for(int t=0;t<timePoint;t++){
for(int u=0;u<field;u++){
......@@ -249,14 +249,14 @@ void reg_smoothImageForCubicSpline1( nifti_image *image,
windowSize = 2*radius[2] + 1;
free(window);
window = (PrecisionTYPE *)calloc(windowSize,sizeof(PrecisionTYPE));
coeffSum=0.0;
// coeffSum=0.0;
for(int it=-radius[2]; it<=radius[2]; it++){
PrecisionTYPE coeff = (PrecisionTYPE)(fabs(2.0*(PrecisionTYPE)it/(PrecisionTYPE)radius[2]));
if(coeff<1.0) window[it+radius[2]] = (PrecisionTYPE)(2.0/3.0 - coeff*coeff + 0.5*coeff*coeff*coeff);
else window[it+radius[2]] = (PrecisionTYPE)(-(coeff-2.0)*(coeff-2.0)*(coeff-2.0)/6.0);
coeffSum += window[it+radius[2]];
// coeffSum += window[it+radius[2]];
}
for(int it=0;it<windowSize;it++)window[it] /= coeffSum;
// for(int it=0;it<windowSize;it++)window[it] /= coeffSum;
for(int t=0;t<timePoint;t++){
for(int u=0;u<field;u++){
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment