Commit fe6c0078 authored by Marc Modat's avatar Marc Modat

2D version of various functions have been added

parent b180140d
......@@ -377,7 +377,7 @@ int main(int argc, char **argv)
printf("\t%ix%ix%i voxels\n",sourceHeader->nx,sourceHeader->ny,sourceHeader->nz);
printf("\t%gx%gx%g mm\n",sourceHeader->dx,sourceHeader->dy,sourceHeader->dz);
printf("Maximum iteration number: %i\n",param->maxIteration);
printf("Percentage of blocks: %i\%\n",param->block_percent_to_use);
printf("Percentage of blocks: %i %%\n",param->block_percent_to_use);
#ifdef _USE_CUDA
if(flag->useGPUFlag) printf("The GPU implementation is used\n");
else printf("The CPU implementation is used\n");
......@@ -565,7 +565,7 @@ int main(int argc, char **argv)
else reg_mat44_disp(&sourceImage->qto_xyz, "[VERBOSE] Source image matrix (qform qto_xyz)");
#endif
printf("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n");
reg_mat44_disp(affineTransformation, "Initial affine transformation:");
reg_mat44_disp(affineTransformation, (char *)"Initial affine transformation:");
/* ****************** */
/* Rigid registration */
......@@ -775,7 +775,7 @@ int main(int argc, char **argv)
NULL,
3,
param->sourceBGValue);
if(!flag->outputResultFlag) param->outputResultName="outputResult.nii";
if(!flag->outputResultFlag) param->outputResultName=(char *)"outputResult.nii";
nifti_set_filenames(resultImage, param->outputResultName, 0, 0);
nifti_image_write(resultImage);
......@@ -784,7 +784,7 @@ int main(int argc, char **argv)
nifti_image_free(resultImage);
nifti_image_free(targetImage);
nifti_image_free(sourceImage);
reg_mat44_disp(affineTransformation, "Final affine transformation:");
reg_mat44_disp(affineTransformation, (char *)"Final affine transformation:");
#ifdef _VERBOSE
mat33 tempMat;
for(int i=0; i<3; i++){
......@@ -801,7 +801,7 @@ int main(int argc, char **argv)
/* The affine transformation is saved */
if(flag->outputAffineFlag)
reg_tool_WriteAffineFile(affineTransformation, param->outputAffineName);
else reg_tool_WriteAffineFile(affineTransformation, "outputAffine.txt");
else reg_tool_WriteAffineFile(affineTransformation, (char *)"outputAffine.txt");
free(affineTransformation);
nifti_image_free(targetHeader);
......
This diff is collapsed.
......@@ -210,7 +210,10 @@ int main(int argc, char **argv)
positionFieldImage->dim[2]=positionFieldImage->ny=targetImage->ny;
positionFieldImage->dim[3]=positionFieldImage->nz=targetImage->nz;
positionFieldImage->dim[4]=positionFieldImage->nt=1;positionFieldImage->pixdim[4]=positionFieldImage->dt=1.0;
positionFieldImage->dim[5]=positionFieldImage->nu=3;positionFieldImage->pixdim[5]=positionFieldImage->du=1.0;
if(targetImage->nz>1)
positionFieldImage->dim[5]=positionFieldImage->nu=3;
else positionFieldImage->dim[5]=positionFieldImage->nu=2;
positionFieldImage->pixdim[5]=positionFieldImage->du=1.0;
positionFieldImage->dim[6]=positionFieldImage->nv=1;positionFieldImage->pixdim[6]=positionFieldImage->dv=1.0;
positionFieldImage->dim[7]=positionFieldImage->nw=1;positionFieldImage->pixdim[7]=positionFieldImage->dw=1.0;
positionFieldImage->nvox=positionFieldImage->nx*positionFieldImage->ny*positionFieldImage->nz*positionFieldImage->nt*positionFieldImage->nu;
......@@ -387,21 +390,37 @@ int main(int argc, char **argv)
displacementFieldImage->data = (void *)calloc(displacementFieldImage->nvox, displacementFieldImage->nbyper);
nifti_set_filenames(displacementFieldImage, param->outputDispName, 0, 0);
memcpy(displacementFieldImage->data, positionFieldImage->data, displacementFieldImage->nvox*displacementFieldImage->nbyper);
PrecisionTYPE *fullDefPtrX=static_cast<PrecisionTYPE *>(displacementFieldImage->data);
PrecisionTYPE *fullDefPtrY=&fullDefPtrX[targetImage->nvox];
PrecisionTYPE *fullDefPtrZ=&fullDefPtrY[targetImage->nvox];
PrecisionTYPE position[3];
for(int z=0; z<displacementFieldImage->nz; z++){
if(targetImage->nz>1){
PrecisionTYPE *fullDefPtrX=static_cast<PrecisionTYPE *>(displacementFieldImage->data);
PrecisionTYPE *fullDefPtrY=&fullDefPtrX[targetImage->nvox];
PrecisionTYPE *fullDefPtrZ=&fullDefPtrY[targetImage->nvox];
PrecisionTYPE position[3];
for(int z=0; z<displacementFieldImage->nz; z++){
for(int y=0; y<displacementFieldImage->ny; y++){
for(int x=0; x<displacementFieldImage->nx; x++){
position[0]=x*targetImage->qto_xyz.m[0][0] + y*targetImage->qto_xyz.m[0][1] + z*targetImage->qto_xyz.m[0][2] + targetImage->qto_xyz.m[0][3];
position[1]=x*targetImage->qto_xyz.m[1][0] + y*targetImage->qto_xyz.m[1][1] + z*targetImage->qto_xyz.m[1][2] + targetImage->qto_xyz.m[1][3];
position[2]=x*targetImage->qto_xyz.m[2][0] + y*targetImage->qto_xyz.m[2][1] + z*targetImage->qto_xyz.m[2][2] + targetImage->qto_xyz.m[2][3];
*fullDefPtrX++ -= position[0];
*fullDefPtrY++ -= position[1];
*fullDefPtrZ++ -= position[2];
}
}
}
}
else{
PrecisionTYPE *fullDefPtrX=static_cast<PrecisionTYPE *>(displacementFieldImage->data);
PrecisionTYPE *fullDefPtrY=&fullDefPtrX[targetImage->nvox];
PrecisionTYPE position[3];
for(int y=0; y<displacementFieldImage->ny; y++){
for(int x=0; x<displacementFieldImage->nx; x++){
position[0]=x*targetImage->qto_xyz.m[0][0] + y*targetImage->qto_xyz.m[0][1] + z*targetImage->qto_xyz.m[0][2] + targetImage->qto_xyz.m[0][3];
position[1]=x*targetImage->qto_xyz.m[1][0] + y*targetImage->qto_xyz.m[1][1] + z*targetImage->qto_xyz.m[1][2] + targetImage->qto_xyz.m[1][3];
position[2]=x*targetImage->qto_xyz.m[2][0] + y*targetImage->qto_xyz.m[2][1] + z*targetImage->qto_xyz.m[2][2] + targetImage->qto_xyz.m[2][3];
position[0]=x*targetImage->qto_xyz.m[0][0] + y*targetImage->qto_xyz.m[0][1] + targetImage->qto_xyz.m[0][3];
position[1]=x*targetImage->qto_xyz.m[1][0] + y*targetImage->qto_xyz.m[1][1] + targetImage->qto_xyz.m[1][3];
*fullDefPtrX++ -= position[0];
*fullDefPtrY++ -= position[1];
*fullDefPtrZ++ -= position[2];
}
}
}
nifti_image_write(displacementFieldImage);
nifti_image_free(displacementFieldImage);
......
......@@ -1410,7 +1410,40 @@ template double reg_bspline_jacobian<double>(nifti_image *, nifti_image *, int);
/* *************************************************************** */
/* *************************************************************** */
template<class NodeTYPE, class VoxelTYPE>
void reg_voxelCentric2NodeCentric1( nifti_image *nodeImage,
void reg_voxelCentric2NodeCentric2D( nifti_image *nodeImage,
nifti_image *voxelImage
)
{
NodeTYPE *nodePtrX = static_cast<NodeTYPE *>(nodeImage->data);
NodeTYPE *nodePtrY = &nodePtrX[nodeImage->nx*nodeImage->ny*nodeImage->nz];
VoxelTYPE *voxelPtrX = static_cast<VoxelTYPE *>(voxelImage->data);
VoxelTYPE *voxelPtrY = &voxelPtrX[voxelImage->nx*voxelImage->ny*voxelImage->nz];
float ratio[2];
ratio[0] = nodeImage->dx / voxelImage->dx;
ratio[1] = nodeImage->dy / voxelImage->dy;
for(int y=0;y<nodeImage->ny; y++){
int Y = (int)round((float)(y-1) * ratio[1]);
VoxelTYPE *yvoxelPtrX=&voxelPtrX[Y*voxelImage->nx];
VoxelTYPE *yvoxelPtrY=&voxelPtrY[Y*voxelImage->nx];
for(int x=0;x<nodeImage->nx; x++){
int X = (int)round((float)(x-1) * ratio[0]);
if(Y<voxelImage->ny && -1<X && X<voxelImage->nx){
*nodePtrX++ = (NodeTYPE)(yvoxelPtrX[X]);
*nodePtrY++ = (NodeTYPE)(yvoxelPtrY[X]);
}
else{
*nodePtrX++ = 0.0;
*nodePtrY++ = 0.0;
}
}
}
}
/* *************************************************************** */
template<class NodeTYPE, class VoxelTYPE>
void reg_voxelCentric2NodeCentric3D( nifti_image *nodeImage,
nifti_image *voxelImage
)
{
......@@ -1460,36 +1493,71 @@ void reg_voxelCentric2NodeCentric( nifti_image *nodeImage,
)
{
// it is assumed than node[000] and voxel[000] are aligned.
switch(nodeImage->datatype){
case NIFTI_TYPE_FLOAT32:
switch(voxelImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_voxelCentric2NodeCentric1<float, float>(nodeImage, voxelImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_voxelCentric2NodeCentric1<float, double>(nodeImage, voxelImage);
break;
default:
printf("err\treg_voxelCentric2NodeCentric:v1\tdata type not supported\n");
break;
}
break;
case NIFTI_TYPE_FLOAT64:
switch(voxelImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_voxelCentric2NodeCentric1<double, float>(nodeImage, voxelImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_voxelCentric2NodeCentric1<double, double>(nodeImage, voxelImage);
break;
default:
printf("err\treg_voxelCentric2NodeCentric:v2\tdata type not supported\n");
break;
}
break;
default:
printf("err\treg_voxelCentric2NodeCentric:n\tdata type not supported\n");
break;
if(nodeImage->nz==1){
switch(nodeImage->datatype){
case NIFTI_TYPE_FLOAT32:
switch(voxelImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_voxelCentric2NodeCentric2D<float, float>(nodeImage, voxelImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_voxelCentric2NodeCentric2D<float, double>(nodeImage, voxelImage);
break;
default:
printf("err\treg_voxelCentric2NodeCentric:v1\tdata type not supported\n");
break;
}
break;
case NIFTI_TYPE_FLOAT64:
switch(voxelImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_voxelCentric2NodeCentric2D<double, float>(nodeImage, voxelImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_voxelCentric2NodeCentric2D<double, double>(nodeImage, voxelImage);
break;
default:
printf("err\treg_voxelCentric2NodeCentric:v2\tdata type not supported\n");
break;
}
break;
default:
printf("err\treg_voxelCentric2NodeCentric:n\tdata type not supported\n");
break;
}
}
else{
switch(nodeImage->datatype){
case NIFTI_TYPE_FLOAT32:
switch(voxelImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_voxelCentric2NodeCentric3D<float, float>(nodeImage, voxelImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_voxelCentric2NodeCentric3D<float, double>(nodeImage, voxelImage);
break;
default:
printf("err\treg_voxelCentric2NodeCentric:v1\tdata type not supported\n");
break;
}
break;
case NIFTI_TYPE_FLOAT64:
switch(voxelImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_voxelCentric2NodeCentric3D<double, float>(nodeImage, voxelImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_voxelCentric2NodeCentric3D<double, double>(nodeImage, voxelImage);
break;
default:
printf("err\treg_voxelCentric2NodeCentric:v2\tdata type not supported\n");
break;
}
break;
default:
printf("err\treg_voxelCentric2NodeCentric:n\tdata type not supported\n");
break;
}
}
}
/* *************************************************************** */
......@@ -2578,10 +2646,11 @@ void reg_bspline_refineControlPointGrid2D( nifti_image *targetImage,
int oldDim[4];
oldDim[1]=splineControlPoint->dim[1];
oldDim[2]=splineControlPoint->dim[2];
oldDim[3]=splineControlPoint->dim[3];
oldDim[3]=1;
splineControlPoint->dx = splineControlPoint->pixdim[1] = splineControlPoint->dx / 2.0f;
splineControlPoint->dy = splineControlPoint->pixdim[2] = splineControlPoint->dy / 2.0f;
splineControlPoint->dz = 1.0f;
splineControlPoint->dim[1]=splineControlPoint->nx=(int)floor(targetImage->nx*targetImage->dx/splineControlPoint->dx)+4;
splineControlPoint->dim[2]=splineControlPoint->ny=(int)floor(targetImage->ny*targetImage->dy/splineControlPoint->dy)+4;
......@@ -2610,17 +2679,17 @@ void reg_bspline_refineControlPointGrid2D( nifti_image *targetImage,
+ 6.0f * (GetValue(oldGridPtrX,oldDim,x-1,y,0) + GetValue(oldGridPtrX,oldDim,x+1,y,0) +
GetValue(oldGridPtrX,oldDim,x,y-1,0) + GetValue(oldGridPtrX,oldDim,x,y+1,0) )
+ 36.0f * GetValue(oldGridPtrX,oldDim,x,y,0) ) / 64.0f);
// 1 0
// 1 0
SetValue(gridPtrX, splineControlPoint->dim, X+1, Y, 0,
(GetValue(oldGridPtrX,oldDim,x,y-1,0) + GetValue(oldGridPtrX,oldDim,x+1,y-1,0) +
GetValue(oldGridPtrX,oldDim,x,y+1,0) + GetValue(oldGridPtrX,oldDim,x+1,y+1,0)
+ 6.0f * GetValue(oldGridPtrX,oldDim,x,y,0) + GetValue(oldGridPtrX,oldDim,x+1,y,0) ) / 16.0f);
// 0 1
+ 6.0f * ( GetValue(oldGridPtrX,oldDim,x,y,0) + GetValue(oldGridPtrX,oldDim,x+1,y,0) ) ) / 16.0f);
// 0 1
SetValue(gridPtrX, splineControlPoint->dim, X, Y+1, 0,
(GetValue(oldGridPtrX,oldDim,x-1,y,0) + GetValue(oldGridPtrX,oldDim,x-1,y+1,0) +
GetValue(oldGridPtrX,oldDim,x+1,y,0) + GetValue(oldGridPtrX,oldDim,x+1,y+1,0)
+ 6.0f * GetValue(oldGridPtrX,oldDim,x,y,0) + GetValue(oldGridPtrX,oldDim,x,y+1,0) ) / 16.0f);
// 1 1
+ 6.0f * ( GetValue(oldGridPtrX,oldDim,x,y,0) + GetValue(oldGridPtrX,oldDim,x,y+1,0) ) ) / 16.0f);
// 1 1
SetValue(gridPtrX, splineControlPoint->dim, X+1, Y+1, 0,
(GetValue(oldGridPtrX,oldDim,x,y,0) + GetValue(oldGridPtrX,oldDim,x+1,y,0) +
GetValue(oldGridPtrX,oldDim,x,y+1,0) + GetValue(oldGridPtrX,oldDim,x+1,y+1,0) ) / 4.0f);
......@@ -2633,17 +2702,17 @@ void reg_bspline_refineControlPointGrid2D( nifti_image *targetImage,
+ 6.0f * (GetValue(oldGridPtrY,oldDim,x-1,y,0) + GetValue(oldGridPtrY,oldDim,x+1,y,0) +
GetValue(oldGridPtrY,oldDim,x,y-1,0) + GetValue(oldGridPtrY,oldDim,x,y+1,0) )
+ 36.0f * GetValue(oldGridPtrY,oldDim,x,y,0) ) / 64.0f);
// 1 0
// 1 0
SetValue(gridPtrY, splineControlPoint->dim, X+1, Y, 0,
(GetValue(oldGridPtrY,oldDim,x,y-1,0) + GetValue(oldGridPtrY,oldDim,x+1,y-1,0) +
GetValue(oldGridPtrY,oldDim,x,y+1,0) + GetValue(oldGridPtrY,oldDim,x+1,y+1,0)
+ 6.0f * GetValue(oldGridPtrY,oldDim,x,y,0) + GetValue(oldGridPtrY,oldDim,x+1,y,0) ) / 16.0f);
// 0 1
+ 6.0f * ( GetValue(oldGridPtrY,oldDim,x,y,0) + GetValue(oldGridPtrY,oldDim,x+1,y,0) ) ) / 16.0f);
// 0 1
SetValue(gridPtrY, splineControlPoint->dim, X, Y+1, 0,
(GetValue(oldGridPtrY,oldDim,x-1,y,0) + GetValue(oldGridPtrY,oldDim,x-1,y+1,0) +
GetValue(oldGridPtrY,oldDim,x+1,y,0) + GetValue(oldGridPtrY,oldDim,x+1,y+1,0)
+ 6.0f * GetValue(oldGridPtrY,oldDim,x,y,0) + GetValue(oldGridPtrY,oldDim,x,y+1,0) ) / 16.0f);
// 1 1
+ 6.0f * ( GetValue(oldGridPtrY,oldDim,x,y,0) + GetValue(oldGridPtrY,oldDim,x,y+1,0) ) ) / 16.0f);
// 1 1
SetValue(gridPtrY, splineControlPoint->dim, X+1, Y+1, 0,
(GetValue(oldGridPtrY,oldDim,x,y,0) + GetValue(oldGridPtrY,oldDim,x+1,y,0) +
GetValue(oldGridPtrY,oldDim,x,y+1,0) + GetValue(oldGridPtrY,oldDim,x+1,y+1,0) ) / 4.0f);
......@@ -3455,7 +3524,37 @@ void reg_bspline_GetJacobianMatrix( nifti_image *splineControlPoint,
/* *************************************************************** */
/* *************************************************************** */
template <class DTYPE>
void reg_bspline_initialiseControlPointGridWithAffine1( mat44 *affineTransformation,
void reg_bspline_initialiseControlPointGridWithAffine2D( mat44 *affineTransformation,
nifti_image *controlPointImage)
{
DTYPE *CPPX=static_cast<DTYPE *>(controlPointImage->data);
DTYPE *CPPY=&CPPX[controlPointImage->nx*controlPointImage->ny*controlPointImage->nz];
mat44 *cppMatrix;
if(controlPointImage->sform_code>0)
cppMatrix=&(controlPointImage->sto_xyz);
else cppMatrix=&(controlPointImage->qto_xyz);
mat44 voxelToRealDeformed = reg_mat44_mul(affineTransformation, cppMatrix);
float index[3];
float position[3];
index[2]=0;
for(int y=0; y<controlPointImage->ny; y++){
index[1]=(float)y;
for(int x=0; x<controlPointImage->nx; x++){
index[0]=(float)x;
reg_mat44_mul(&voxelToRealDeformed, index, position);
*CPPX++ = position[0];
*CPPY++ = position[1];
}
}
}
/* *************************************************************** */
template <class DTYPE>
void reg_bspline_initialiseControlPointGridWithAffine3D( mat44 *affineTransformation,
nifti_image *controlPointImage)
{
DTYPE *CPPX=static_cast<DTYPE *>(controlPointImage->data);
......@@ -3491,17 +3590,33 @@ void reg_bspline_initialiseControlPointGridWithAffine1( mat44 *affineTransformat
int reg_bspline_initialiseControlPointGridWithAffine( mat44 *affineTransformation,
nifti_image *controlPointImage)
{
switch(controlPointImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_bspline_initialiseControlPointGridWithAffine1<float>(affineTransformation, controlPointImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_bspline_initialiseControlPointGridWithAffine1<double>(affineTransformation, controlPointImage);
break;
default:
fprintf(stderr,"ERROR:\treg_bspline_initialiseControlPointGridWithAffine\n");
fprintf(stderr,"ERROR:\tOnly single of double precision is implemented for the control point image\n");
return 1;
if(controlPointImage->nz==1){
switch(controlPointImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_bspline_initialiseControlPointGridWithAffine2D<float>(affineTransformation, controlPointImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_bspline_initialiseControlPointGridWithAffine2D<double>(affineTransformation, controlPointImage);
break;
default:
fprintf(stderr,"ERROR:\treg_bspline_initialiseControlPointGridWithAffine\n");
fprintf(stderr,"ERROR:\tOnly single or double precision is implemented for the control point image\n");
return 1;
}
}
else{
switch(controlPointImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_bspline_initialiseControlPointGridWithAffine3D<float>(affineTransformation, controlPointImage);
break;
case NIFTI_TYPE_FLOAT64:
reg_bspline_initialiseControlPointGridWithAffine3D<double>(affineTransformation, controlPointImage);
break;
default:
fprintf(stderr,"ERROR:\treg_bspline_initialiseControlPointGridWithAffine\n");
fprintf(stderr,"ERROR:\tOnly single or double precision is implemented for the control point image\n");
return 1;
}
}
return 0;
}
......
......@@ -327,7 +327,104 @@ template void reg_getEntropies<double>(nifti_image *, nifti_image *, int, int, d
/* *************************************************************** */
/* *************************************************************** */
template<class PrecisionTYPE,class TargetTYPE,class ResultTYPE,class ResultGradientTYPE,class NMIGradientTYPE>
void reg_getVoxelBasedNMIGradientUsingPW4( nifti_image *targetImage,
void reg_getVoxelBasedNMIGradientUsingPW2D( nifti_image *targetImage,
nifti_image *resultImage,
int type,
nifti_image *resultImageGradient,
int binning,
PrecisionTYPE *logJointHistogram,
PrecisionTYPE *entropies,
nifti_image *nmiGradientImage,
int *mask)
{
TargetTYPE *targetPtr = static_cast<TargetTYPE *>(targetImage->data);
ResultTYPE *resultPtr = static_cast<ResultTYPE *>(resultImage->data);
ResultGradientTYPE *resultGradientPtrX = static_cast<ResultGradientTYPE *>(resultImageGradient->data);
ResultGradientTYPE *resultGradientPtrY = &resultGradientPtrX[resultImage->nvox];
NMIGradientTYPE *nmiGradientPtrX = static_cast<NMIGradientTYPE *>(nmiGradientImage->data);
NMIGradientTYPE *nmiGradientPtrY = &nmiGradientPtrX[resultImage->nvox];
int *maskPtr = &mask[0];
// In a first time the NMI gradient is computed for every voxel
memset(nmiGradientPtrX,0,nmiGradientImage->nvox*nmiGradientImage->nbyper);
PrecisionTYPE NMI = (entropies[0]+entropies[1])/entropies[2];
unsigned int binningSquare = binning*binning;
for(int y=0; y<targetImage->ny; y++){
for(int x=0; x<targetImage->nx; x++){
if(*maskPtr++>-1){
TargetTYPE targetValue = *targetPtr;
ResultTYPE resultValue = *resultPtr;
if(targetValue>2.0f && resultValue>2.0f){
// The two is added because the image is resample between 2 and bin +2
// if 64 bins are used the histogram will have 68 bins et the image will be between 2 and 65
if(type!=1){
targetValue = (TargetTYPE)floor((double)targetValue);
resultValue = (ResultTYPE)floor((double)resultValue);
}
PrecisionTYPE resDeriv[2];
resDeriv[0] = (PrecisionTYPE)(*resultGradientPtrX);
resDeriv[1] = (PrecisionTYPE)(*resultGradientPtrY);
PrecisionTYPE jointEntropyDerivative_X = 0.0;
PrecisionTYPE movingEntropyDerivative_X = 0.0;
PrecisionTYPE fixedEntropyDerivative_X = 0.0;
PrecisionTYPE jointEntropyDerivative_Y = 0.0;
PrecisionTYPE movingEntropyDerivative_Y = 0.0;
PrecisionTYPE fixedEntropyDerivative_Y = 0.0;
for(int t=(int)(targetValue-1.0); t<(int)(targetValue+2.0); t++){
if(-1<t && t<binning){
for(int r=(int)(resultValue-1.0); r<(int)(resultValue+2.0); r++){
if(-1<r && r<binning){
PrecisionTYPE commonValue = GetBasisSplineValue<PrecisionTYPE>((PrecisionTYPE)t-(PrecisionTYPE)targetValue) *
GetBasisSplineDerivativeValue<PrecisionTYPE>((PrecisionTYPE)r-(PrecisionTYPE)resultValue);
PrecisionTYPE jointLog = logJointHistogram[t*binning+r];
PrecisionTYPE targetLog = logJointHistogram[binningSquare+t];
PrecisionTYPE resultLog = logJointHistogram[binningSquare+binning+r];
PrecisionTYPE temp = commonValue * resDeriv[0];
jointEntropyDerivative_X -= temp * jointLog;
fixedEntropyDerivative_X -= temp * targetLog;
movingEntropyDerivative_X -= temp * resultLog;
temp = commonValue * resDeriv[1];
jointEntropyDerivative_Y -= temp * jointLog;
fixedEntropyDerivative_Y -= temp * targetLog;
movingEntropyDerivative_Y -= temp * resultLog;
} // O<t<bin
} // 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
*nmiGradientPtrX = (NMIGradientTYPE)((fixedEntropyDerivative_X + movingEntropyDerivative_X - NMI * jointEntropyDerivative_X) / temp);
*nmiGradientPtrY = (NMIGradientTYPE)((fixedEntropyDerivative_Y + movingEntropyDerivative_Y - NMI * jointEntropyDerivative_Y) / temp);
} // value > 0
}// mask > -1
targetPtr++;
resultPtr++;
nmiGradientPtrX++;
nmiGradientPtrY++;
resultGradientPtrX++;
resultGradientPtrY++;
}
}
}
/* *************************************************************** */
template<class PrecisionTYPE,class TargetTYPE,class ResultTYPE,class ResultGradientTYPE,class NMIGradientTYPE>
void reg_getVoxelBasedNMIGradientUsingPW3D( nifti_image *targetImage,
nifti_image *resultImage,
int type,
nifti_image *resultImageGradient,
......@@ -452,18 +549,34 @@ void reg_getVoxelBasedNMIGradientUsingPW3( nifti_image *targetImage,
nifti_image *nmiGradientImage,
int *mask)
{
switch(nmiGradientImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_getVoxelBasedNMIGradientUsingPW4<PrecisionTYPE,TargetTYPE,ResultTYPE,ResultGradientTYPE,float>
(targetImage, resultImage, type, resultImageGradient, binning, logJointHistogram, entropies, nmiGradientImage, mask);
break;
case NIFTI_TYPE_FLOAT64:
reg_getVoxelBasedNMIGradientUsingPW4<PrecisionTYPE,TargetTYPE,ResultTYPE,ResultGradientTYPE,double>
(targetImage, resultImage, type, resultImageGradient, binning, logJointHistogram, entropies, nmiGradientImage, mask);
break;
default:
printf("err\treg_getVoxelBasedNMIGradientUsingPW\tThe result image gradient data type is not supported\n");
return;
if(nmiGradientImage->nz>1){
switch(nmiGradientImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_getVoxelBasedNMIGradientUsingPW3D<PrecisionTYPE,TargetTYPE,ResultTYPE,ResultGradientTYPE,float>
(targetImage, resultImage, type, resultImageGradient, binning, logJointHistogram, entropies, nmiGradientImage, mask);
break;
case NIFTI_TYPE_FLOAT64:
reg_getVoxelBasedNMIGradientUsingPW3D<PrecisionTYPE,TargetTYPE,ResultTYPE,ResultGradientTYPE,double>
(targetImage, resultImage, type, resultImageGradient, binning, logJointHistogram, entropies, nmiGradientImage, mask);
break;
default:
printf("err\treg_getVoxelBasedNMIGradientUsingPW\tThe result image gradient data type is not supported\n");
return;
}
}else{
switch(nmiGradientImage->datatype){
case NIFTI_TYPE_FLOAT32:
reg_getVoxelBasedNMIGradientUsingPW2D<PrecisionTYPE,TargetTYPE,ResultTYPE,ResultGradientTYPE,float>
(targetImage, resultImage, type, resultImageGradient, binning, logJointHistogram, entropies, nmiGradientImage, mask);
break;
case NIFTI_TYPE_FLOAT64:
reg_getVoxelBasedNMIGradientUsingPW2D<PrecisionTYPE,TargetTYPE,ResultTYPE,ResultGradientTYPE,double>
(targetImage, resultImage, type, resultImageGradient, binning, logJointHistogram, entropies, nmiGradientImage, mask);
break;
default:
printf("err\treg_getVoxelBasedNMIGradientUsingPW\tThe result image gradient data type is not supported\n");
return;
}
}
}
/* *************************************************************** */
......
......@@ -244,45 +244,46 @@ void reg_smoothImageForCubicSpline1( nifti_image *image,
}
}
}
/* Smoothing along the Z axis */
windowSize = 2*radius[2] + 1;
free(window);
window = (PrecisionTYPE *)calloc(windowSize,sizeof(PrecisionTYPE));
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]];
}
for(int it=0;it<windowSize;it++)window[it] /= coeffSum;
for(int t=0;t<timePoint;t++){
for(int u=0;u<field;u++){
DTYPE *readingValue=&imageArray[(t+u*timePoint)*image->nx*image->ny*image->nz];
DTYPE *writtingValue=&tempArray[(t+u*timePoint)*image->nx*image->ny*image->nz];
int i=0;
for(int z=0; z<image->nz; z++){
for(int y=0; y<image->ny; y++){
for(int x=0; x<image->nx; x++){
PrecisionTYPE finalValue=0.0;
int index = i - image->nx*image->ny*radius[2];
int Z = z - radius[2];
for(int it=0; it<windowSize; it++){
if(-1<