Commit d86a84ed authored by Marc Modat's avatar Marc Modat

Error when using a deformation field with an intersection

parent b08eaf4c
......@@ -138,6 +138,8 @@ int main(int argc, char **argv)
// Create the average image
nifti_image *average_image=nifti_copy_nim_info(tempImage);
average_image->scl_slope=1.f;
average_image->scl_inter=0.f;
nifti_image_free(tempImage);
tempImage=NULL;
average_image->datatype=NIFTI_TYPE_FLOAT32;
......@@ -166,9 +168,9 @@ int main(int argc, char **argv)
reg_print_msg_error(argv[i]);
return EXIT_FAILURE;
}
if(sizeof(PrecisionTYPE)==sizeof(double))
average_norm_intensity<double>(tempImage);
else average_norm_intensity<float>(tempImage);
// if(sizeof(PrecisionTYPE)==sizeof(double))
// average_norm_intensity<double>(tempImage);
// else average_norm_intensity<float>(tempImage);
reg_tools_addImageToImage(average_image,tempImage,average_image);
imageTotalNumber++;
nifti_image_free(tempImage);
......@@ -343,17 +345,11 @@ int main(int argc, char **argv)
averageMatrix = nifti_mat44_inverse(averageMatrix);
averageMatrix = reg_mat44_logm(&averageMatrix);
// Demean all the input affine matrices
float indet=1.f; // HERE
float outdet=1.f; // HERE
for(size_t i=0; i<affineNumber; ++i)
{
indet *= reg_mat44_det(&affineMatrices[i]);// HERE
for(size_t i=0; i<affineNumber; ++i){
affineMatrices[i] = reg_mat44_logm(&affineMatrices[i]);
affineMatrices[i] = averageMatrix + affineMatrices[i];
affineMatrices[i] = reg_mat44_expm(&affineMatrices[i]);
outdet *= reg_mat44_det(&affineMatrices[i]);// HERE
}
printf("Average determinant %g -> %g\n", indet, outdet); // HERE
// Create a deformation field to be used to resample all the floating images
nifti_image *deformationField = nifti_copy_nim_info(referenceImage);
deformationField->dim[0]=deformationField->ndim=5;
......@@ -366,6 +362,8 @@ int main(int argc, char **argv)
deformationField->datatype=NIFTI_TYPE_FLOAT32;
deformationField->nbyper=sizeof(float);
}
deformationField->scl_slope=1.f;
deformationField->scl_inter=0.f;
deformationField->data = (void *)malloc(deformationField->nvox*deformationField->nbyper);
// Create an average image
nifti_image *averageImage = nifti_copy_nim_info(referenceImage);
......@@ -377,6 +375,8 @@ int main(int argc, char **argv)
averageImage->data = (void *)calloc(averageImage->nvox,averageImage->nbyper);
// Create a temporary image
nifti_image *tempImage = nifti_copy_nim_info(averageImage);
tempImage->scl_slope=1.f;
tempImage->scl_inter=0.f;
tempImage->data = (void *)malloc(tempImage->nvox*tempImage->nbyper);
// warp all floating images and sum them up
for(size_t i=5, j=0; i<argc; i+=2,++j)
......@@ -402,9 +402,9 @@ int main(int argc, char **argv)
}
}
reg_resampleImage(floatingImage,tempImage,deformationField,NULL,3,0.f);
if(sizeof(PrecisionTYPE)==sizeof(double))
average_norm_intensity<double>(tempImage);
else average_norm_intensity<float>(tempImage);
// if(sizeof(PrecisionTYPE)==sizeof(double))
// average_norm_intensity<double>(tempImage);
// else average_norm_intensity<float>(tempImage);
reg_tools_addImageToImage(averageImage,tempImage,averageImage);
nifti_image_free(floatingImage);
}
......@@ -435,6 +435,8 @@ int main(int argc, char **argv)
averageField->nbyper=sizeof(float);
}
averageField->data = (void *)calloc(averageField->nvox,averageField->nbyper);
averageField->scl_slope=1.f;
averageField->scl_inter=0.f;
reg_tools_multiplyValueToImage(averageField,averageField,0.f);
// Iterate over all the transformation parametrisations - Note that I don't store them all to save space
#ifndef NDEBUG
......@@ -465,6 +467,8 @@ int main(int argc, char **argv)
nifti_image *deformationField = nifti_copy_nim_info(averageField);
deformationField->data = (void *)malloc(deformationField->nvox*deformationField->nbyper);
reg_tools_multiplyValueToImage(deformationField,deformationField,0.f);
deformationField->scl_slope=1.f;
deformationField->scl_inter=0.f;
deformationField->intent_p1=DISP_FIELD;
reg_getDeformationFromDisplacement(deformationField);
// Generate a deformation field or a flow field depending of the input transformation
......@@ -524,6 +528,8 @@ int main(int argc, char **argv)
// The affine component is substracted
nifti_image *tempField = nifti_copy_nim_info(deformationField);
tempField->data = (void *)malloc(tempField->nvox*tempField->nbyper);
tempField->scl_slope=1.f;
tempField->scl_inter=0.f;
reg_affine_getDeformationField(&affineTransformation,
tempField);
reg_tools_substractImageToImage(deformationField,tempField,deformationField);
......@@ -544,6 +550,8 @@ int main(int argc, char **argv)
averageImage->datatype=NIFTI_TYPE_FLOAT32;
averageImage->nbyper=sizeof(float);
}
averageImage->scl_slope=1.f;
averageImage->scl_inter=0.f;
averageImage->data = (void *)calloc(averageImage->nvox,averageImage->nbyper);
// Create a temporary image
nifti_image *tempImage = nifti_copy_nim_info(averageImage);
......@@ -574,6 +582,8 @@ int main(int argc, char **argv)
reg_tools_multiplyValueToImage(deformationField,deformationField,0.f);
deformationField->intent_p1=DISP_FIELD;
reg_getDeformationFromDisplacement(deformationField);
deformationField->scl_slope=1.f;
deformationField->scl_inter=0.f;
// Generate a deformation field or a flow field depending of the input transformation
switch(static_cast<int>(transformation->intent_p1))
{
......@@ -607,6 +617,8 @@ int main(int argc, char **argv)
nifti_image *tempDef = nifti_copy_nim_info(deformationField);
tempDef->data = (void *)malloc(tempDef->nvox*tempDef->nbyper);
memcpy(tempDef->data,deformationField->data,tempDef->nvox*tempDef->nbyper);
tempDef->scl_slope=1.f;
tempDef->scl_inter=0.f;
reg_defField_getDeformationFieldFromFlowField(tempDef,deformationField,false);
deformationField->intent_p1=DEF_FIELD;
nifti_free_extensions(deformationField);
......
......@@ -226,6 +226,8 @@ int main(int argc, char **argv)
defField->datatype=NIFTI_TYPE_FLOAT32;
defField->nbyper=sizeof(float);
defField->data=(void *)calloc(defField->nvox,defField->nbyper);
defField->scl_slope=1.f;
defField->scl_inter=0.f;
reg_tools_multiplyValueToImage(defField,defField,0.f);
defField->intent_p1=DISP_FIELD;
reg_getDeformationFromDisplacement(defField);
......
......@@ -305,6 +305,8 @@ int main(int argc, char **argv)
deformationFieldImage->nvox =(size_t)deformationFieldImage->nx*
deformationFieldImage->ny*deformationFieldImage->nz*
deformationFieldImage->nt*deformationFieldImage->nu;
deformationFieldImage->scl_slope=1.f;
deformationFieldImage->scl_inter=0.f;
if(inputTransformationImage!=NULL)
{
deformationFieldImage->datatype = inputTransformationImage->datatype;
......
......@@ -331,6 +331,8 @@ int main(int argc, char **argv)
outputTransformationImage->intent_code=NIFTI_INTENT_VECTOR;
memset(outputTransformationImage->intent_name, 0, 16);
strcpy(outputTransformationImage->intent_name,"NREG_TRANS");
outputTransformationImage->scl_slope=1.f;
outputTransformationImage->scl_inter=0.f;
}
else
{
......@@ -619,6 +621,8 @@ int main(int argc, char **argv)
output1TransImage->nvox=(size_t)output1TransImage->nx *
output1TransImage->ny * output1TransImage->nz *
output1TransImage->nt * output1TransImage->nu;
output1TransImage->scl_slope=1.f;
output1TransImage->scl_inter=0.f;
if(referenceImage->datatype!=NIFTI_TYPE_FLOAT32)
{
output1TransImage->nbyper=sizeof(float);
......@@ -744,6 +748,8 @@ int main(int argc, char **argv)
if(referenceImage2!=NULL)
{
output2TransImage=nifti_copy_nim_info(referenceImage2);
output2TransImage->scl_slope=1.f;
output2TransImage->scl_inter=0.f;
printf("[NiftyReg] Transformation 2 is defined in the space of image:\n[NiftyReg] %s\n",
referenceImage2->fname);
}
......@@ -970,6 +976,8 @@ int main(int argc, char **argv)
strcpy(outputTransImage->intent_name,"NREG_TRANS");
outputTransImage->intent_p1=inputTransImage->intent_p1;
outputTransImage->intent_p2=inputTransImage->intent_p2;
outputTransImage->scl_slope=1.f;
outputTransImage->scl_inter=0.f;
outputTransImage->data=(void *)malloc
(outputTransImage->nvox*outputTransImage->nbyper);
// Convert the spline parametrisation into a dense deformation parametrisation
......@@ -1009,6 +1017,8 @@ int main(int argc, char **argv)
tempField->intent_p1=DEF_VEL_FIELD;
tempField->intent_p2=inputTransImage->intent_p2;
}
tempField->scl_slope=1.f;
tempField->scl_inter=0.f;
tempField->data=(void *)calloc(tempField->nvox,tempField->nbyper);
// Compute the dense field
if(inputTransImage->intent_p1==SPLINE_GRID)
......
......@@ -156,12 +156,13 @@ int reg_aladin<T>::Print()
printf("[%s] Reference image name: %s\n", this->ExecutableName, this->InputReference->fname);
printf("[%s] \t%ix%ix%i voxels\n", this->ExecutableName, this->InputReference->nx,this->InputReference->ny,this->InputReference->nz);
printf("[%s] \t%gx%gx%g mm\n", this->ExecutableName, this->InputReference->dx,this->InputReference->dy,this->InputReference->dz);
printf("[%s] floating image name: %s\n", this->ExecutableName, this->InputFloating->fname);
printf("[%s] Floating image name: %s\n", this->ExecutableName, this->InputFloating->fname);
printf("[%s] \t%ix%ix%i voxels\n", this->ExecutableName, this->InputFloating->nx,this->InputFloating->ny,this->InputFloating->nz);
printf("[%s] \t%gx%gx%g mm\n", this->ExecutableName, this->InputFloating->dx,this->InputFloating->dy,this->InputFloating->dz);
printf("[%s] Maximum iteration number: %i", this->ExecutableName, this->MaxIterations);
printf(" (%i during the first level)\n", 2*this->MaxIterations);
printf("[%s] Percentage of blocks: %i %%", this->ExecutableName, this->BlockPercentage);
printf("[%s] Percentage of blocks: %i %%\n", this->ExecutableName, this->BlockPercentage);
printf("* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\n");
#ifdef NDEBUG
}
#endif
......@@ -444,6 +445,8 @@ void reg_aladin<T>::AllocateDeformationField()
fprintf(stderr,"[NiftyReg ERROR] Only float or double are expected for the deformation field. Exit.\n");
reg_exit(1);
}
this->deformationFieldImage->scl_slope=1.f;
this->deformationFieldImage->scl_inter=0.f;
this->deformationFieldImage->data = (void *)calloc(this->deformationFieldImage->nvox, this->deformationFieldImage->nbyper);
return;
}
......
......@@ -549,6 +549,8 @@ void reg_base<T>::AllocateDeformationField()
memset(this->deformationFieldImage->intent_name, 0, 16);
strcpy(this->deformationFieldImage->intent_name,"NREG_TRANS");
this->deformationFieldImage->intent_p1=DEF_FIELD;
this->deformationFieldImage->scl_slope=1.f;
this->deformationFieldImage->scl_inter=0.f;
if(this->measure_dti!=NULL)
this->forwardJacobianMatrix=(mat33 *)malloc(
......
......@@ -253,6 +253,8 @@ void reg_f3d_sym<T>::AllocateDeformationField()
memset(this->backwardDeformationFieldImage->intent_name, 0, 16);
strcpy(this->backwardDeformationFieldImage->intent_name,"NREG_TRANS");
this->backwardDeformationFieldImage->intent_p1=DEF_FIELD;
this->backwardDeformationFieldImage->scl_slope=1.f;
this->backwardDeformationFieldImage->scl_inter=0.f;
if(this->measure_dti!=NULL)
this->backwardJacobianMatrix=(mat33 *)malloc(
......
......@@ -3918,7 +3918,7 @@ void reg_defField_getDeformationFieldFromFlowField(nifti_image *flowFieldImage,
{
// Create a field that contains the affine component only
affineOnly = nifti_copy_nim_info(deformationFieldImage);
affineOnly->data = (void *)malloc(affineOnly->nvox*affineOnly->nbyper);
affineOnly->data = (void *)calloc(affineOnly->nvox,affineOnly->nbyper);
reg_affine_getDeformationField(reinterpret_cast<mat44 *>(flowFieldImage->ext_list[0].edata),
affineOnly,
false);
......@@ -3927,7 +3927,6 @@ void reg_defField_getDeformationFieldFromFlowField(nifti_image *flowFieldImage,
}
else reg_getDisplacementFromDeformation(flowFieldImage);
// Compute the number of scaling value to ensure unfolded transformation
int squaringNumber = 1;
if(updateStepNumber || flowFieldImage->intent_p2==0)
......
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