Commit 79bc94c7 authored by Marc Modat's avatar Marc Modat

Change in the pyramidal approach (again) to ensure minimal downsampling....

Change in the pyramidal approach (again) to ensure minimal downsampling. Thresholding possible on both input images in F3D
parent 9374592f
......@@ -388,6 +388,8 @@ int main(int argc, char **argv)
/* ********************** */
/* START THE REGISTRATION */
/* ********************** */
bool downsampleAxis[8]={false,false,false,false,false,false,false,false};
for(int level=0; level<param->level2Perform; level++){
/* Read the target and source image */
......@@ -417,10 +419,22 @@ int main(int argc, char **argv)
memcpy( tempMaskImage->data, targetMaskImage->data, tempMaskImage->nvox*tempMaskImage->nbyper);
}
for(int l=level; l<param->levelNumber-1; l++){
reg_downsampleImage<PrecisionTYPE>(targetImage, 1);
reg_downsampleImage<PrecisionTYPE>(sourceImage, 1);
int ratio = (int)pow(2,param->levelNumber-param->levelNumber+l+1);
bool sourceDownsampleAxis[8]={true,true,true,true,true,true,true,true};
if((sourceHeader->nx/ratio) < 32) sourceDownsampleAxis[1]=false;
if((sourceHeader->ny/ratio) < 32) sourceDownsampleAxis[2]=false;
if((sourceHeader->nz/ratio) < 32) sourceDownsampleAxis[3]=false;
reg_downsampleImage<PrecisionTYPE>(sourceImage, 1, sourceDownsampleAxis);
bool targetDownsampleAxis[8]={true,true,true,true,true,true,true,true};
if((targetHeader->nx/ratio) < 32) targetDownsampleAxis[1]=false;
if((targetHeader->ny/ratio) < 32) targetDownsampleAxis[2]=false;
if((targetHeader->nz/ratio) < 32) targetDownsampleAxis[3]=false;
reg_downsampleImage<PrecisionTYPE>(targetImage, 1, targetDownsampleAxis);
if(flag->targetMaskFlag){
reg_downsampleImage<PrecisionTYPE>(tempMaskImage, 0);
reg_downsampleImage<PrecisionTYPE>(tempMaskImage, 0, targetDownsampleAxis);
}
}
......@@ -714,7 +728,7 @@ int main(int argc, char **argv)
cudaCommon_free(&sourceImageArray_d);
cudaCommon_free((void **)&resultImageArray_d);
cudaCommon_free((void **)&positionFieldImageArray_d);
cudaCommon_free((void **)activeBlock_d);
cudaCommon_free((void **)&activeBlock_d);
CUDA_SAFE_CALL(cudaFreeHost(resultImage->data));
resultImage->data=NULL;
}
......
......@@ -26,6 +26,7 @@
#include "_reg_mutualinformation.h"
#include "_reg_ssd.h"
#include "_reg_tools.h"
#include "float.h"
#ifdef _USE_CUDA
#include "_reg_cudaCommon.h"
......@@ -64,6 +65,10 @@ typedef struct{
PrecisionTYPE sourceBGValue;
float targetSigmaValue;
float sourceSigmaValue;
float targetLowThresholdValue;
float targetUpThresholdValue;
float sourceLowThresholdValue;
float sourceUpThresholdValue;
}PARAM;
typedef struct{
bool targetImageFlag;
......@@ -94,7 +99,12 @@ typedef struct{
bool noConjugateGradient;
bool targetSigmaFlag;
bool sourceSigmaFlag;
bool pyramidFlag;
bool pyramidFlag;
bool targetLowThresholdFlag;
bool targetUpThresholdFlag;
bool sourceLowThresholdFlag;
bool sourceUpThresholdFlag;
#ifdef _USE_CUDA
bool useGPUFlag;
......@@ -138,15 +148,20 @@ void Usage(char *exec)
printf("\t-sy <float>\t\tFinal grid spacing along the y axis in mm [sx value]\n");
printf("\t-sz <float>\t\tFinal grid spacing along the z axis in mm [sx value]\n");
printf("\t-bin <int>\t\tNumber of bin to use for the joint histogram [64]\n");
printf("\t-smooT <float>\t\tSmooth the target image using the specified sigma (mm) [0]\n");
printf("\t-smooS <float>\t\tSmooth the source image using the specified sigma (mm) [0]\n");
printf("\t-smooT <float>\t\tSmooth the target image using the specified sigma (mm) [0]\n");
printf("\t-smooS <float>\t\tSmooth the source image using the specified sigma (mm) [0]\n");
printf("\t-tLwTh <float>\t\tLower threshold to apply to the target image intensities [none]\n");
printf("\t-tUpTh <float>\t\tUpper threshold to apply to the target image intensities [none]\n");
printf("\t-sLwTh <float>\t\tLower threshold to apply to the source image intensities [none]\n");
printf("\t-sUpTh <float>\t\tUpper threshold to apply to the source image intensities [none]\n");
printf("\t\t\t\tThe scl_slope and scl_inter from the nifti header are taken into account for the thresholds\n");
printf("\t-ln <int>\t\tNumber of level to perform [3]\n");
printf("\t-lp <int>\t\tOnly perform the first levels [ln]\n");
printf("\t-nopy\t\t\tDo not use a pyramidal approach [no]\n");
printf("\t-be <float>\t\tWeight of the bending energy penalty term [0.01]\n");
printf("\t-noAppBE\t\t\tTo not approximate the BE value only at the control point position\n");
printf("\t-noGradBE\t\t\tTo not use the gradient of the bending energy\n");
printf("\t-noAppBE\t\tTo not approximate the BE value only at the control point position\n");
printf("\t-noGradBE\t\tTo not use the gradient of the bending energy\n");
// printf("\t-jl <float>\t\tWeight of log of the Jacobian determinant penalty term [0.0]\n");
// printf("\t-appJL\t\t\tApproximate the JL value only at the control point position [no]\n");
......@@ -178,6 +193,11 @@ int main(int argc, char **argv)
flag->appBendingEnergyFlag=1;
flag->beGradFlag=1;
param->targetLowThresholdValue=-FLT_MAX;
param->targetUpThresholdValue=FLT_MAX;
param->sourceLowThresholdValue=-FLT_MAX;
param->sourceUpThresholdValue=FLT_MAX;
/* read the input parameter */
for(int i=1;i<argc;i++){
if(strcmp(argv[i], "-help")==0 || strcmp(argv[i], "-Help")==0 ||
......@@ -279,10 +299,26 @@ int main(int argc, char **argv)
param->targetSigmaValue=(float)(atof(argv[++i]));
flag->targetSigmaFlag=1;
}
else if(strcmp(argv[i], "-smooS") == 0){
param->sourceSigmaValue=(float)(atof(argv[++i]));
flag->sourceSigmaFlag=1;
}
else if(strcmp(argv[i], "-smooS") == 0){
param->sourceSigmaValue=(float)(atof(argv[++i]));
flag->sourceSigmaFlag=1;
}
else if(strcmp(argv[i], "-tLwTh") == 0){
param->targetLowThresholdValue=(float)(atof(argv[++i]));
flag->targetLowThresholdFlag=1;
}
else if(strcmp(argv[i], "-tUpTh") == 0){
param->targetUpThresholdValue=(float)(atof(argv[++i]));
flag->targetUpThresholdFlag=1;
}
else if(strcmp(argv[i], "-sLwTh") == 0){
param->sourceLowThresholdValue=(float)(atof(argv[++i]));
flag->sourceLowThresholdFlag=1;
}
else if(strcmp(argv[i], "-sUpTh") == 0){
param->sourceUpThresholdValue=(float)(atof(argv[++i]));
flag->sourceUpThresholdFlag=1;
}
else if(strcmp(argv[i], "-bgi") == 0){
param->backgroundIndex[0]=atoi(argv[++i]);
param->backgroundIndex[1]=atoi(argv[++i]);
......@@ -309,6 +345,7 @@ int main(int argc, char **argv)
}
}
#ifdef _USE_CUDA
if(flag->useGPUFlag){
if(flag->jacobianWeightFlag){
......@@ -469,7 +506,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("Number of bin to used: %i\n",param->binning);
printf("Number of bin to used: %i\n",param->binning-4);
printf("Bending energy weight: %g\n",param->bendingEnergyWeight);
if(flag->appBendingEnergyFlag) printf("Bending energy penalty term evaluated at the control point position only\n");
printf("log of the jacobian determinant weight: %g\n",param->jacobianWeight);
......@@ -554,11 +591,24 @@ int main(int argc, char **argv)
tempMaskImage->data = (void *)malloc(tempMaskImage->nvox * tempMaskImage->nbyper);
memcpy(tempMaskImage->data, targetMaskImage->data, tempMaskImage->nvox*tempMaskImage->nbyper);
}
for(int l=level; l<param->levelNumber-1; l++){
reg_downsampleImage<PrecisionTYPE>(targetImage, 1);
reg_downsampleImage<PrecisionTYPE>(sourceImage, 1);
for(int l=level; l<param->levelNumber-1; l++){
int ratio = (int)pow(2,param->levelNumber-param->levelNumber+l+1);
bool sourceDownsampleAxis[8]={true,true,true,true,true,true,true,true};
if((sourceHeader->nx/ratio) < 32) sourceDownsampleAxis[1]=false;
if((sourceHeader->ny/ratio) < 32) sourceDownsampleAxis[2]=false;
if((sourceHeader->nz/ratio) < 32) sourceDownsampleAxis[3]=false;
reg_downsampleImage<PrecisionTYPE>(sourceImage, 1, sourceDownsampleAxis);
bool targetDownsampleAxis[8]={true,true,true,true,true,true,true,true};
if((targetHeader->nx/ratio) < 32) targetDownsampleAxis[1]=false;
if((targetHeader->ny/ratio) < 32) targetDownsampleAxis[2]=false;
if((targetHeader->nz/ratio) < 32) targetDownsampleAxis[3]=false;
reg_downsampleImage<PrecisionTYPE>(targetImage, 1, targetDownsampleAxis);
if(flag->targetMaskFlag){
reg_downsampleImage<PrecisionTYPE>(tempMaskImage, 0);
reg_downsampleImage<PrecisionTYPE>(tempMaskImage, 0, targetDownsampleAxis);
}
}
targetMask = (int *)malloc(targetImage->nvox*sizeof(int));
......@@ -584,6 +634,12 @@ int main(int argc, char **argv)
reg_gaussianSmoothing<PrecisionTYPE>(sourceImage, param->sourceSigmaValue, smoothAxis);
}
/* the target and source are resampled between 0 and bin-1
* The images are then shifted by two which is the suport of the spline used
* by the parzen window filling of the joint histogram */
reg_intensityRescale(targetImage,2.0f,(float)param->binning-3.0f, param->targetLowThresholdValue, param->targetUpThresholdValue);
reg_intensityRescale(sourceImage,2.0f,(float)param->binning-3.0f, param->sourceLowThresholdValue, param->sourceUpThresholdValue);
if(level==0){
if(!flag->inputCPPFlag){
/* allocate the control point image */
......@@ -728,18 +784,12 @@ int main(int argc, char **argv)
reg_mat44_disp(cppMatrix_xyz, "[VERBOSE] Control point image matrix");
#endif
printf("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n");
float maxStepSize = (targetImage->dx>targetImage->dy)?targetImage->dx:targetImage->dy;
maxStepSize = (targetImage->dz>maxStepSize)?targetImage->dz:maxStepSize;
float currentSize = maxStepSize;
float smallestSize = maxStepSize / 100.0f;
/* the target and source are resampled between 0 and bin-1
* The images are then shifted by two which is the suport of the spline used
* by the parzen window filling of the joint histogram */
reg_intensityRescale(targetImage,2.0f,(float)param->binning-3.0f);
reg_intensityRescale(sourceImage,2.0f,(float)param->binning-3.0f);
if(flag->backgroundIndexFlag){
int index[3];
index[0]=param->backgroundIndex[0];
......
......@@ -27,8 +27,10 @@ int round(PrecisionType x)
template<class DTYPE>
void reg_intensityRescale2( nifti_image *image,
float newMin,
float newMax
float newMin,
float newMax,
float lowThr,
float upThr
)
{
DTYPE *imagePtr = static_cast<DTYPE *>(image->data);
......@@ -41,7 +43,7 @@ void reg_intensityRescale2( nifti_image *image,
break;
case NIFTI_TYPE_INT8:
currentMin=(DTYPE)std::numeric_limits<char>::max();
currentMax=(DTYPE)std::numeric_limits<char>::min();
currentMax=(DTYPE)-std::numeric_limits<char>::max();
break;
case NIFTI_TYPE_UINT16:
currentMin=(DTYPE)std::numeric_limits<unsigned short>::max();
......@@ -49,7 +51,7 @@ void reg_intensityRescale2( nifti_image *image,
break;
case NIFTI_TYPE_INT16:
currentMin=(DTYPE)std::numeric_limits<char>::max();
currentMax=(DTYPE)std::numeric_limits<char>::min();
currentMax=-(DTYPE)std::numeric_limits<char>::max();
break;
case NIFTI_TYPE_UINT32:
currentMin=(DTYPE)std::numeric_limits<unsigned int>::max();
......@@ -57,65 +59,81 @@ void reg_intensityRescale2( nifti_image *image,
break;
case NIFTI_TYPE_INT32:
currentMin=(DTYPE)std::numeric_limits<int>::max();
currentMax=(DTYPE)std::numeric_limits<int>::min();
currentMax=-(DTYPE)std::numeric_limits<int>::max();
break;
case NIFTI_TYPE_FLOAT32:
currentMin=(DTYPE)std::numeric_limits<float>::max();
currentMax=(DTYPE)std::numeric_limits<float>::min();
currentMax=-(DTYPE)std::numeric_limits<float>::max();
break;
case NIFTI_TYPE_FLOAT64:
currentMin=(DTYPE)std::numeric_limits<double>::max();
currentMax=(DTYPE)std::numeric_limits<double>::min();
currentMax=-(DTYPE)std::numeric_limits<double>::max();
break;
}
for(unsigned int index=0; index<image->nvox; index++){
DTYPE value = *imagePtr++;
DTYPE value = (DTYPE)(*imagePtr++ * image->scl_slope + image->scl_inter);
currentMin=(currentMin<value)?currentMin:value;
currentMax=(currentMax>value)?currentMax:value;
}
}
if(currentMin<lowThr) currentMin=(DTYPE)lowThr;
if(currentMax>upThr) currentMax=(DTYPE)upThr;
double currentDiff = (double)(currentMax-currentMin);
double newDiff = (double)(newMax-newMin);
image->cal_min=newMin * image->scl_slope + image->scl_inter;
image->cal_max=newMax * image->scl_slope + image->scl_inter;
imagePtr = static_cast<DTYPE *>(image->data);
for(unsigned int index=0; index<image->nvox; index++){
double value = (double)*imagePtr;
value = (value-(double)currentMin)/currentDiff;
value = value * newDiff + newMin;
double value = (double)*imagePtr * image->scl_slope + image->scl_inter;
if(value<currentMin){
value = newMin;
}
else if(value>currentMax){
value = newMax;
}
else{
value = (value-(double)currentMin)/currentDiff;
value = value * newDiff + newMin;
}
*imagePtr++=(DTYPE)value;
}
}
void reg_intensityRescale( nifti_image *image,
float newMin,
float newMax
float newMin,
float newMax,
float lowThr,
float upThr
)
{
switch(image->datatype){
case NIFTI_TYPE_UINT8:
reg_intensityRescale2<unsigned char>(image, newMin, newMax);
reg_intensityRescale2<unsigned char>(image, newMin, newMax, lowThr, upThr);
break;
case NIFTI_TYPE_INT8:
reg_intensityRescale2<char>(image, newMin, newMax);
reg_intensityRescale2<char>(image, newMin, newMax, lowThr, upThr);
break;
case NIFTI_TYPE_UINT16:
reg_intensityRescale2<unsigned short>(image, newMin, newMax);
reg_intensityRescale2<unsigned short>(image, newMin, newMax, lowThr, upThr);
break;
case NIFTI_TYPE_INT16:
reg_intensityRescale2<short>(image, newMin, newMax);
reg_intensityRescale2<short>(image, newMin, newMax, lowThr, upThr);
break;
case NIFTI_TYPE_UINT32:
reg_intensityRescale2<unsigned int>(image, newMin, newMax);
reg_intensityRescale2<unsigned int>(image, newMin, newMax, lowThr, upThr);
break;
case NIFTI_TYPE_INT32:
reg_intensityRescale2<int>(image, newMin, newMax);
reg_intensityRescale2<int>(image, newMin, newMax, lowThr, upThr);
break;
case NIFTI_TYPE_FLOAT32:
reg_intensityRescale2<float>(image, newMin, newMax);
reg_intensityRescale2<float>(image, newMin, newMax, lowThr, upThr);
break;
case NIFTI_TYPE_FLOAT64:
reg_intensityRescale2<double>(image, newMin, newMax);
reg_intensityRescale2<double>(image, newMin, newMax, lowThr, upThr);
break;
default:
printf("err\treg_intensityRescale\tThe image data type is not supported\n");
......@@ -952,17 +970,11 @@ template void reg_gaussianSmoothing<double>(nifti_image *, float, bool[8]);
/* *************************************************************** */
/* *************************************************************** */
template <class PrecisionTYPE, class ImageTYPE>
void reg_downsampleImage1(nifti_image *image, int type)
void reg_downsampleImage1(nifti_image *image, int type, bool downsampleAxis[8])
{
bool downXYZ[8]={false,false,false,false,false,false,false,false};
for(int i=1;i<=image->dim[0];i++){
if((image->dim[i]/2)>=32) downXYZ[i]=true;
}
if(type==1){
/* the input image is first smooth */
reg_gaussianSmoothing<float>(image, -0.7f, downXYZ);
reg_gaussianSmoothing<float>(image, -0.7f, downsampleAxis);
}
/* the values are copied */
......@@ -974,8 +986,8 @@ void reg_downsampleImage1(nifti_image *image, int type)
int oldDim[4];
for(int i=1; i<4; i++){
oldDim[i]=image->dim[i];
if(image->dim[i]>1 && downXYZ[i]==true) image->dim[i]=(int)(image->dim[i]/2.0);
if(image->pixdim[i]>0 && downXYZ[i]==true) image->pixdim[i]=image->pixdim[i]*2.0f;
if(image->dim[i]>1 && downsampleAxis[i]==true) image->dim[i]=(int)(image->dim[i]/2.0);
if(image->pixdim[i]>0 && downsampleAxis[i]==true) image->pixdim[i]=image->pixdim[i]*2.0f;
}
image->nx=image->dim[1];
image->ny=image->dim[2];
......@@ -991,7 +1003,7 @@ void reg_downsampleImage1(nifti_image *image, int type)
for(int i=0; i<3; i++){
for(int j=0; j<3; j++){
oldMat.m[i][j]=image->qto_ijk.m[i][j];
if(downXYZ[j+1]==true){
if(downsampleAxis[j+1]==true){
image->qto_xyz.m[i][j]=image->qto_xyz.m[i][j]*2.0f;
image->sto_xyz.m[i][j]=image->sto_xyz.m[i][j]*2.0f;
}
......@@ -1099,40 +1111,40 @@ void reg_downsampleImage1(nifti_image *image, int type)
}
/* *************************************************************** */
template <class PrecisionTYPE>
void reg_downsampleImage(nifti_image *image, int type)
void reg_downsampleImage(nifti_image *image, int type, bool downsampleAxis[8])
{
switch(image->datatype){
case NIFTI_TYPE_UINT8:
reg_downsampleImage1<PrecisionTYPE,unsigned char>(image, type);
reg_downsampleImage1<PrecisionTYPE,unsigned char>(image, type, downsampleAxis);
break;
case NIFTI_TYPE_INT8:
reg_downsampleImage1<PrecisionTYPE,char>(image, type);
reg_downsampleImage1<PrecisionTYPE,char>(image, type, downsampleAxis);
break;
case NIFTI_TYPE_UINT16:
reg_downsampleImage1<PrecisionTYPE,unsigned short>(image, type);
reg_downsampleImage1<PrecisionTYPE,unsigned short>(image, type, downsampleAxis);
break;
case NIFTI_TYPE_INT16:
reg_downsampleImage1<PrecisionTYPE,short>(image, type);
reg_downsampleImage1<PrecisionTYPE,short>(image, type, downsampleAxis);
break;
case NIFTI_TYPE_UINT32:
reg_downsampleImage1<PrecisionTYPE,unsigned int>(image, type);
reg_downsampleImage1<PrecisionTYPE,unsigned int>(image, type, downsampleAxis);
break;
case NIFTI_TYPE_INT32:
reg_downsampleImage1<PrecisionTYPE,int>(image, type);
reg_downsampleImage1<PrecisionTYPE,int>(image, type, downsampleAxis);
break;
case NIFTI_TYPE_FLOAT32:
reg_downsampleImage1<PrecisionTYPE,float>(image, type);
reg_downsampleImage1<PrecisionTYPE,float>(image, type, downsampleAxis);
break;
case NIFTI_TYPE_FLOAT64:
reg_downsampleImage1<PrecisionTYPE,double>(image, type);
reg_downsampleImage1<PrecisionTYPE,double>(image, type, downsampleAxis);
break;
default:
printf("err\treg_downsampleImage\tThe image data type is not supported\n");
return;
}
}
template void reg_downsampleImage<float>(nifti_image *, int);
template void reg_downsampleImage<double>(nifti_image *, int);
template void reg_downsampleImage<float>(nifti_image *, int, bool[8]);
template void reg_downsampleImage<double>(nifti_image *, int, bool[8]);
/* *************************************************************** */
/* *************************************************************** */
template <class DTYPE>
......
......@@ -19,7 +19,9 @@
extern "C++"
void reg_intensityRescale( nifti_image *image,
float newMin,
float newMax
float newMax,
float lowThr,
float upThr
);
extern "C++" template <class PrecisionTYPE>
......@@ -38,7 +40,7 @@ void reg_gaussianSmoothing( nifti_image *image,
bool[8]);
extern "C++" template <class PrecisionTYPE>
void reg_downsampleImage(nifti_image *image, int);
void reg_downsampleImage(nifti_image *image, int, bool[8]);
extern "C++" template <class PrecisionTYPE>
PrecisionTYPE reg_getMaximalLength(nifti_image *image);
......
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