Commit fe5c06a6 authored by Marc Modat's avatar Marc Modat

Minor error correction in the GPU BE gradient computation. extrapolation of...

Minor error correction in the GPU BE gradient computation. extrapolation of the joint histogram [0;64[ -> [0;68[
parent cf56e43d
......@@ -25,6 +25,7 @@ ENDIF(WIN32)
OPTION(BUILD_ALADIN "To build the reg_aladin executable" ON)
OPTION(BUILD_F3D "To build the reg_f3d executable" ON)
OPTION(BUILD_RESAMPLE "To build the reg_resample executable" ON)
OPTION(BUILD_TOOLS "To build the reg_tools executable" ON)
OPTION(USE_VERBOSE "To print out extra information" OFF)
......
#-----------------------------------------------------------------------------
IF(BUILD_RESAMPLE)
ADD_EXECUTABLE(reg_resample reg_resample.cpp)
TARGET_LINK_LIBRARIES(reg_resample _reg_resampling _reg_bspline _reg_tools _reg_affineTransformation reg_nifti ${ZLIB})
INSTALL_TARGETS(/bin reg_resample)
IF(BUILD_TOOLS)
ADD_EXECUTABLE(reg_tools reg_tools.cpp)
TARGET_LINK_LIBRARIES(reg_tools _reg_resampling _reg_tools _reg_affineTransformation reg_nifti ${ZLIB})
INSTALL_TARGETS(/bin reg_tools)
ENDIF(BUILD_TOOLS)
#-----------------------------------------------------------------------------
IF(BUILD_RESAMPLE)
ADD_EXECUTABLE(reg_resample reg_resample.cpp)
TARGET_LINK_LIBRARIES(reg_resample _reg_resampling _reg_bspline _reg_tools _reg_affineTransformation reg_nifti ${ZLIB})
INSTALL_TARGETS(/bin reg_resample)
ENDIF(BUILD_RESAMPLE)
#-----------------------------------------------------------------------------
IF(BUILD_F3D)
......
......@@ -134,7 +134,7 @@ void Usage(char *exec)
printf("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n");
return;
}
#define CONVERGENCE_EPS 0.001
#define CONVERGENCE_EPS 0.00001
bool reg_test_convergence(mat44 *mat)
{
bool convergence=true;
......
......@@ -41,6 +41,9 @@
#endif
#define PrecisionTYPE float
#define JH_TRI 0
#define JH_PARZEN_WIN 1
#define JH_PW_APPROX 2
typedef struct{
char *targetImageName;
......@@ -92,7 +95,6 @@ typedef struct{
bool targetSigmaFlag;
bool sourceSigmaFlag;
bool pyramidFlag;
bool metricPaddingFlag;
#ifdef _USE_CUDA
bool useGPUFlag;
......@@ -152,7 +154,6 @@ void Usage(char *exec)
printf("\t-bgi <int> <int> <int>\tForce the background value during\n\t\t\t\tresampling to have the same value as this voxel in the source image [none]\n");
// printf("\t-ssd\t\t\tTo use the SSD as the similiarity measure [no]\n");
printf("\t-noConj\t\t\tTo not use the conjuage gradient optimisation but a simple gradient ascent/descent\n");
printf("\t-mpad\t\t\tTo include the padding value (0.0) into the metric computation\n");
printf("\t-mem\t\t\tDisplay an approximate memory requierment and exit\n");
#ifdef _USE_CUDA
......@@ -280,9 +281,6 @@ int main(int argc, char **argv)
param->sourceSigmaValue=(float)(atof(argv[++i]));
flag->sourceSigmaFlag=1;
}
else if(strcmp(argv[i], "-mpad") == 0){
flag->metricPaddingFlag=1;
}
else if(strcmp(argv[i], "-bgi") == 0){
param->backgroundIndex[0]=atoi(argv[++i]);
param->backgroundIndex[1]=atoi(argv[++i]);
......@@ -343,6 +341,7 @@ int main(int argc, char **argv)
/* Read the number of bin to use */
if(!flag->binningFlag) param->binning=64;
param->binning += 4; //This is due to the extrapolation of the joint histogram using the Parzen window
if(param->bendingEnergyWeight==0.0f)flag->bendingEnergyFlag=0;
......@@ -567,7 +566,8 @@ int main(int argc, char **argv)
}
else{
for(int i=0; i<targetImage->nvox; i++)
targetMask[i]=activeVoxelNumber++;
targetMask[i]=i;
activeVoxelNumber=targetImage->nvox;
}
}
......@@ -728,57 +728,59 @@ int main(int argc, char **argv)
float currentSize = maxStepSize;
float smallestSize = maxStepSize / 100.0f;
/* the target and source are resampled between 0 and bin-1 */
reg_intensityRescale(targetImage,0.0f,(float)param->binning-1.0f);
reg_intensityRescale(sourceImage,0.0f,(float)param->binning-1.0f);
if(flag->backgroundIndexFlag){
int index[3];
index[0]=param->backgroundIndex[0];
index[1]=param->backgroundIndex[1];
index[2]=param->backgroundIndex[2];
if(flag->pyramidFlag){
for(int l=level; l<param->levelNumber-1; l++){
index[0] /= 2;
index[1] /= 2;
index[2] /= 2;
}
}
param->sourceBGValue = (float)(reg_tool_GetIntensityValue(sourceImage, index));
}
/* 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+1.0f);
reg_intensityRescale(sourceImage,2.0f,(float)param->binning+1.0f);
if(flag->backgroundIndexFlag){
int index[3];
index[0]=param->backgroundIndex[0];
index[1]=param->backgroundIndex[1];
index[2]=param->backgroundIndex[2];
if(flag->pyramidFlag){
for(int l=level; l<param->levelNumber-1; l++){
index[0] /= 2;
index[1] /= 2;
index[2] /= 2;
}
}
param->sourceBGValue = (float)(reg_tool_GetIntensityValue(sourceImage, index));
}
/* the gradient images are allocated */
nifti_image *resultGradientImage;
nifti_image *voxelNMIGradientImage;
nifti_image *nodeNMIGradientImage;
/* Conjugate gradient */
PrecisionTYPE *conjugateG;
PrecisionTYPE *conjugateH;
/* joint histogram related variables */
double *probaJointHistogram = (double *)malloc(param->binning*(param->binning+2)*sizeof(double));
double *logJointHistogram = (double *)malloc(param->binning*(param->binning+2)*sizeof(double));
double *entropies = (double *)malloc(4*sizeof(double));
PrecisionTYPE *bestControlPointPosition;
/* the gradient images are allocated */
nifti_image *resultGradientImage;
nifti_image *voxelNMIGradientImage;
nifti_image *nodeNMIGradientImage;
/* Conjugate gradient */
PrecisionTYPE *conjugateG;
PrecisionTYPE *conjugateH;
/* joint histogram related variables */
double *probaJointHistogram = (double *)malloc(param->binning*(param->binning+2)*sizeof(double));
double *logJointHistogram = (double *)malloc(param->binning*(param->binning+2)*sizeof(double));
double *entropies = (double *)malloc(4*sizeof(double));
PrecisionTYPE *bestControlPointPosition;
#ifdef _USE_CUDA
float *targetImageArray_d;
cudaArray *sourceImageArray_d;
float4 *controlPointImageArray_d;
float *resultImageArray_d;
float4 *positionFieldImageArray_d;
float *targetImageArray_d;
cudaArray *sourceImageArray_d;
float4 *controlPointImageArray_d;
float *resultImageArray_d;
float4 *positionFieldImageArray_d;
int *targetMask_d;
float4 *resultGradientArray_d;
float4 *voxelNMIGradientArray_d;
float4 *nodeNMIGradientArray_d;
float4 *resultGradientArray_d;
float4 *voxelNMIGradientArray_d;
float4 *nodeNMIGradientArray_d;
float4 *conjugateG_d;
float4 *conjugateH_d;
float4 *conjugateG_d;
float4 *conjugateH_d;
float4 *bestControlPointPosition_d;
float4 *bestControlPointPosition_d;
float *logJointHistogram_d;
float *logJointHistogram_d;
if(flag->useGPUFlag){
if(!flag->noConjugateGradient){
......@@ -850,7 +852,7 @@ int main(int argc, char **argv)
smoothingRadius[2] = (int)floor( 2.0*controlPointImage->dz/targetImage->dz );
int iteration=0;
while(iteration<param->maxIteration && currentSize>smallestSize){
while(iteration<param->maxIteration && currentSize>smallestSize){
#ifdef _USE_CUDA
if(flag->useGPUFlag){
/* generate the position field */
......@@ -894,12 +896,11 @@ int main(int argc, char **argv)
double currentValue;
reg_getEntropies<double>( targetImage,
resultImage,
2,
JH_PW_APPROX,
param->binning,
probaJointHistogram,
logJointHistogram,
entropies,
flag->metricPaddingFlag,
targetMask);
currentValue = (1.0-param->bendingEnergyWeight-param->jacobianWeight)*(entropies[0]+entropies[1])/entropies[2];
......@@ -981,19 +982,18 @@ int main(int argc, char **argv)
&targetMask_d,
activeVoxelNumber,
entropies,
param->binning,
flag->metricPaddingFlag);
reg_smoothImageForCubicSpline_gpu( resultImage,
&voxelNMIGradientArray_d,
smoothingRadius);
reg_voxelCentric2NodeCentric_gpu( targetImage,
controlPointImage,
&voxelNMIGradientArray_d,
&nodeNMIGradientArray_d);
param->binning);
reg_smoothImageForCubicSpline_gpu( resultImage,
&voxelNMIGradientArray_d,
smoothingRadius);
reg_voxelCentric2NodeCentric_gpu( targetImage,
controlPointImage,
&voxelNMIGradientArray_d,
&nodeNMIGradientArray_d);
/* The NMI gradient is converted from voxel space to real space */
reg_convertNMIGradientFromVoxelToRealSpace_gpu( sourceMatrix_xyz,
controlPointImage,
&nodeNMIGradientArray_d);
reg_convertNMIGradientFromVoxelToRealSpace_gpu( sourceMatrix_xyz,
controlPointImage,
&nodeNMIGradientArray_d);
/* The other gradients are calculated */
if(flag->beGradFlag && flag->bendingEnergyFlag){
reg_bspline_ApproxBendingEnergyGradient_gpu(controlPointImage,
......@@ -1033,15 +1033,15 @@ int main(int argc, char **argv)
positionFieldImage,
targetMask,
1);
reg_getVoxelBasedNMIGradientUsingPW<double>( targetImage,
resultImage,
resultGradientImage,
param->binning,
logJointHistogram,
entropies,
voxelNMIGradientImage,
targetMask,
flag->metricPaddingFlag);
reg_getVoxelBasedNMIGradientUsingPW<double>(targetImage,
resultImage,
JH_PW_APPROX,
resultGradientImage,
param->binning,
logJointHistogram,
entropies,
voxelNMIGradientImage,
targetMask);
reg_smoothImageForCubicSpline<PrecisionTYPE>(voxelNMIGradientImage,smoothingRadius);
reg_voxelCentric2NodeCentric(nodeNMIGradientImage,voxelNMIGradientImage);
......@@ -1239,12 +1239,11 @@ int main(int argc, char **argv)
/* Computation of the NMI */
reg_getEntropies<double>( targetImage,
resultImage,
2,
JH_PW_APPROX,
param->binning,
probaJointHistogram,
logJointHistogram,
entropies,
flag->metricPaddingFlag,
targetMask);
currentValue = (1.0-param->bendingEnergyWeight-param->jacobianWeight)*(entropies[0]+entropies[1])/entropies[2];
......@@ -1341,20 +1340,28 @@ int main(int argc, char **argv)
free(probaJointHistogram);
free(logJointHistogram);
#ifdef _USE_CUDA
if(flag->useGPUFlag){
if(cudaCommon_transferFromDeviceToNifti(controlPointImage, &controlPointImageArray_d)) return 1;
cudaCommon_free((void **)&resultGradientArray_d);
cudaCommon_free((void **)&voxelNMIGradientArray_d);
cudaCommon_free((void **)&nodeNMIGradientArray_d);
if(!flag->noConjugateGradient){
cudaCommon_free((void **)&conjugateG_d);
cudaCommon_free((void **)&conjugateH_d);
}
cudaCommon_free((void **)&bestControlPointPosition_d);
cudaCommon_free((void **)&logJointHistogram_d);
if(flag->useGPUFlag){
if(cudaCommon_transferFromDeviceToNifti(controlPointImage, &controlPointImageArray_d)) return 1;
cudaCommon_free( (void **)&targetImageArray_d );
cudaCommon_free( &sourceImageArray_d );
cudaCommon_free( (void **)&controlPointImageArray_d );
cudaCommon_free( (void **)&resultImageArray_d );
cudaCommon_free( (void **)&positionFieldImageArray_d );
CUDA_SAFE_CALL(cudaFree(targetMask_d));
}
else{
cudaCommon_free((void **)&resultGradientArray_d);
cudaCommon_free((void **)&voxelNMIGradientArray_d);
cudaCommon_free((void **)&nodeNMIGradientArray_d);
if(!flag->noConjugateGradient){
cudaCommon_free((void **)&conjugateG_d);
cudaCommon_free((void **)&conjugateH_d);
}
cudaCommon_free((void **)&bestControlPointPosition_d);
cudaCommon_free((void **)&logJointHistogram_d);
CUDA_SAFE_CALL(cudaFreeHost(resultImage->data));
resultImage->data = NULL;
}
else{
#endif
nifti_image_free(resultGradientImage);
nifti_image_free(voxelNMIGradientImage);
......@@ -1368,18 +1375,12 @@ int main(int argc, char **argv)
}
#endif
#ifdef _USE_CUDA
if(flag->useGPUFlag){
CUDA_SAFE_CALL(cudaFreeHost(resultImage->data));
resultImage->data = NULL;
}
#endif
nifti_image_free( resultImage );
nifti_image_free( resultImage );
if(level==(param->level2Perform-1)){
/* ****************** */
/* OUTPUT THE RESULTS */
/* ****************** */
if(level==(param->level2Perform-1)){
/* ****************** */
/* OUTPUT THE RESULTS */
/* ****************** */
/* The best result is returned */
nifti_set_filenames(controlPointImage, param->outputCPPName, 0, 0);
......@@ -1434,15 +1435,6 @@ int main(int argc, char **argv)
nifti_image_free( sourceImage );
nifti_image_free( targetImage );
#ifdef _USE_CUDA
if(flag->useGPUFlag){
cudaCommon_free( (void **)&targetImageArray_d );
cudaCommon_free( &sourceImageArray_d );
cudaCommon_free( (void **)&controlPointImageArray_d );
cudaCommon_free( (void **)&resultImageArray_d );
cudaCommon_free( (void **)&positionFieldImageArray_d );
}
#endif
printf("- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n\n");
} // for(int level=0; level<param->levelNumber; level++){
......
......@@ -21,13 +21,15 @@
typedef struct{
char *inputImageName;
char *outputImageName;
char *addImageName;
char *addImageName;
char *rmsImageName;
int smoothValue;
}PARAM;
typedef struct{
bool inputImageFlag;
bool outputImageFlag;
bool addImageFlag;
bool addImageFlag;
bool rmsImageFlag;
bool smoothValueFlag;
bool gradientImageFlag;
}FLAG;
......@@ -44,11 +46,12 @@ void Usage(char *exec)
printf("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n");
printf("Usage:\t%s -in <filename> -out <filename> [OPTIONS].\n",exec);
printf("\t-in <filename>\tFilename of the input image image (mandatory)\n");
printf("\t-out <filename>\t\tFilename out the output image [output.nii]\n");
printf("* * OPTIONS * *\n");
printf("\t-out <filename>\t\tFilename out the output image [output.nii]\n");
printf("\t-grad\t\t\t4D spatial gradient of the input image\n");
printf("\t-add <filename>\t\tThis image is added to the input\n");
printf("\t-smo <int>\t\tThe input image is smoothed using a b-spline curve\n");
printf("\t-smo <int>\t\tThe input image is smoothed using a b-spline curve\n");
printf("\t-rms <filename>\tCompute the mean rms between both image\n");
printf("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n");
return;
}
......@@ -77,10 +80,14 @@ int main(int argc, char **argv)
else if(strcmp(argv[i], "-grad") == 0){
flag->gradientImageFlag=1;
}
else if(strcmp(argv[i], "-add") == 0){
param->addImageName=argv[++i];
flag->addImageFlag=1;
}
else if(strcmp(argv[i], "-add") == 0){
param->addImageName=argv[++i];
flag->addImageFlag=1;
}
else if(strcmp(argv[i], "-rms") == 0){
param->rmsImageName=argv[++i];
flag->rmsImageFlag=1;
}
else if(strcmp(argv[i], "-smo") == 0){
param->smoothValue=atoi(argv[++i]);
flag->smoothValueFlag=1;
......@@ -124,18 +131,19 @@ int main(int argc, char **argv)
affineTransformation->m[1][1]=1.0;
affineTransformation->m[2][2]=1.0;
affineTransformation->m[3][3]=1.0;
nifti_image *fakeDeformationField = nifti_copy_nim_info(spatialGradientImage);
fakeDeformationField->data = (void *)malloc(fakeDeformationField->nvox * fakeDeformationField->nbyper);
reg_affine_deformationField( affineTransformation,
image,
fakeDeformationField);
nifti_image *fakepositionField = nifti_copy_nim_info(spatialGradientImage);
fakepositionField->data = (void *)malloc(fakepositionField->nvox * fakepositionField->nbyper);
reg_affine_positionField( affineTransformation,
image,
fakepositionField);
free(affineTransformation);
reg_getSourceImageGradient<PrecisionTYPE>( image,
image,
spatialGradientImage,
fakeDeformationField,
3); // cubic spline interpolation
nifti_image_free(fakeDeformationField);
image,
spatialGradientImage,
fakepositionField,
NULL,
3); // cubic spline interpolation
nifti_image_free(fakepositionField);
nifti_image_write(spatialGradientImage);
nifti_image_free(spatialGradientImage);
}
......@@ -150,40 +158,63 @@ int main(int argc, char **argv)
printf("%i\n", param->smoothValue);
int radius[3];radius[0]=radius[1]=radius[2]=param->smoothValue;
reg_smoothImageForCubicSpline<PrecisionTYPE>(smoothImg, radius);
// reg_smoothImageForTrilinear<PrecisionTYPE>(smoothImg, radius);
nifti_image_write(smoothImg);
nifti_image_free(smoothImg);
}
if(flag->addImageFlag){
nifti_image *imageToAdd = nifti_image_read(param->addImageName,true);
if(imageToAdd == NULL){
fprintf(stderr,"** ERROR Error when reading the image to add: %s\n",param->addImageName);
return 1;
}
// Check image dimension
if(image->dim[0]!=imageToAdd->dim[0] ||
image->dim[1]!=imageToAdd->dim[1] ||
image->dim[2]!=imageToAdd->dim[2] ||
image->dim[3]!=imageToAdd->dim[3] ||
image->dim[4]!=imageToAdd->dim[4] ||
image->dim[5]!=imageToAdd->dim[5] ||
image->dim[6]!=imageToAdd->dim[6] ||
image->dim[7]!=imageToAdd->dim[7]){
fprintf(stderr,"Both images do not have the same dimension\n");
return 1;
}
nifti_image *sumImage = nifti_copy_nim_info(image);
sumImage->data = (void *)malloc(sumImage->nvox * sumImage->nbyper);
if(flag->outputImageFlag)
nifti_set_filenames(sumImage, param->outputImageName, 0, 0);
else nifti_set_filenames(sumImage, "output.nii", 0, 0);
reg_tools_addImages(image, imageToAdd, sumImage);
nifti_image_write(sumImage);
nifti_image_free(sumImage);
nifti_image_free(imageToAdd);
}
if(flag->addImageFlag){
nifti_image *imageToAdd = nifti_image_read(param->addImageName,true);
if(imageToAdd == NULL){
fprintf(stderr,"** ERROR Error when reading the image to add: %s\n",param->addImageName);
return 1;
}
// Check image dimension
if(image->dim[0]!=imageToAdd->dim[0] ||
image->dim[1]!=imageToAdd->dim[1] ||
image->dim[2]!=imageToAdd->dim[2] ||
image->dim[3]!=imageToAdd->dim[3] ||
image->dim[4]!=imageToAdd->dim[4] ||
image->dim[5]!=imageToAdd->dim[5] ||
image->dim[6]!=imageToAdd->dim[6] ||
image->dim[7]!=imageToAdd->dim[7]){
fprintf(stderr,"Both images do not have the same dimension\n");
return 1;
}
nifti_image *sumImage = nifti_copy_nim_info(image);
sumImage->data = (void *)malloc(sumImage->nvox * sumImage->nbyper);
if(flag->outputImageFlag)
nifti_set_filenames(sumImage, param->outputImageName, 0, 0);
else nifti_set_filenames(sumImage, "output.nii", 0, 0);
reg_tools_addImages(image, imageToAdd, sumImage);
nifti_image_write(sumImage);
nifti_image_free(sumImage);
nifti_image_free(imageToAdd);
}
if(flag->rmsImageFlag){
nifti_image *image2 = nifti_image_read(param->rmsImageName,true);
if(image2 == NULL){
fprintf(stderr,"** ERROR Error when reading the image to add: %s\n",param->rmsImageName);
return 1;
}
// Check image dimension
if(image->dim[0]!=image2->dim[0] ||
image->dim[1]!=image2->dim[1] ||
image->dim[2]!=image2->dim[2] ||
image->dim[3]!=image2->dim[3] ||
image->dim[4]!=image2->dim[4] ||
image->dim[5]!=image2->dim[5] ||
image->dim[6]!=image2->dim[6] ||
image->dim[7]!=image2->dim[7]){
fprintf(stderr,"Both images do not have the same dimension\n");
return 1;
}
double meanRMSerror = reg_tools_getMeanRMS(image, image2);
printf("%g\n", meanRMSerror);
nifti_image_free(image2);
}
nifti_image_free(image);
return 0;
......
......@@ -1318,9 +1318,10 @@ void reg_bspline_bendingEnergyGradient1( nifti_image *splineControlPoint,
metricGradientValue[0] = (PrecisionTYPE)(*gradientXPtr);
metricGradientValue[1] = (PrecisionTYPE)(*gradientYPtr);
metricGradientValue[2] = (PrecisionTYPE)(*gradientZPtr);
*gradientXPtr++ = (SplineTYPE)((1.0-weight)*metricGradientValue[0] + weight*gradientValue[0]/nodeNumber);
*gradientYPtr++ = (SplineTYPE)((1.0-weight)*metricGradientValue[1] + weight*gradientValue[1]/nodeNumber);
*gradientZPtr++ = (SplineTYPE)((1.0-weight)*metricGradientValue[2] + weight*gradientValue[2]/nodeNumber);
// (Marc) I removed the normalisation by the voxel number as each gradient has to be normalised in the same way
*gradientXPtr++ = (SplineTYPE)((1.0-weight)*metricGradientValue[0] + weight*gradientValue[0]);
*gradientYPtr++ = (SplineTYPE)((1.0-weight)*metricGradientValue[1] + weight*gradientValue[1]);
*gradientZPtr++ = (SplineTYPE)((1.0-weight)*metricGradientValue[2] + weight*gradientValue[2]);
}
}
}
......
This diff is collapsed.
......@@ -22,19 +22,16 @@ void reg_getEntropies( nifti_image *targetImage,
PrecisionTYPE *probaJointHistogram,
PrecisionTYPE *logJointHistogram,
PrecisionTYPE *entropies,
bool includePadding,
int *mask
);
int *mask);
extern "C++" template <class PrecisionTYPE>
void reg_getVoxelBasedNMIGradientUsingPW( nifti_image *targetImage,
nifti_image *resultImage,
int type,
nifti_image *resultImageGradient,
int binning,
PrecisionTYPE *logJointHistogram,
PrecisionTYPE *entropies,
nifti_image *nmiGradientImage,
int *mask,
bool includePadding
);
int *mask);
#endif
......@@ -141,20 +141,15 @@ void reg_smoothImageForCubicSpline1( nifti_image *image,
/* Smoothing along the X axis */
int windowSize = 2*radius[0] + 1;
// printf("window size along X: %i\n", windowSize);
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++){
// printf("coeff[%i] = %g -> ", it, window[it]);
//window[it] /= coeffSum;
// printf("%g\n", window[it]);
// }
for(int it=0;it<windowSize;it++) window[it] /= coeffSum;
for(int t=0;t<timePoint;t++){
for(int u=0;u<field;u++){
......@@ -189,17 +184,16 @@ void reg_smoothImageForCubicSpline1( nifti_image *image,