Commit 24582f8e authored by Marc Modat's avatar Marc Modat

I changed F3D2 in order to include an affine transformation within the CPP...

I changed F3D2 in order to include an affine transformation within the CPP file itself, as a nifti extension.
parent e69cc4cf
......@@ -147,8 +147,7 @@ void Usage(char *exec)
printf("\t-ic <float>\t\tWeight of the inverse consistency penalty term [0.01]\n");
printf("\n*** F3D2 options:\n");
printf("\t-vel \t\t\tUse a velocity field integrationto generate the deformation\n");
printf("\t-step <int>\t\tNumber of composition step [6].\n");
printf("\t-vel \t\t\tUse a velocity field integration to generate the deformation\n");
printf("\n*** Other options:\n");
printf("\t-smoothGrad <float>\tTo smooth the metric derivative (in mm) [0]\n");
......@@ -249,11 +248,11 @@ int main(int argc, char **argv)
#endif // _USE_CUDA
reg_f3d<PrecisionTYPE> *REG=NULL;
for(int i=1;i<argc;i++){
if(strcmp(argv[i], "-vel")==0 || strcmp(argv[i], "-step")==0){
if(strcmp(argv[i], "-vel")==0 || strcmp(argv[i], "--vel")==0){
REG=new reg_f3d2<PrecisionTYPE>(referenceImage->nt,floatingImage->nt);
break;
}
if(strcmp(argv[i], "-sym")==0){
if(strcmp(argv[i], "-sym")==0 || strcmp(argv[i], "--sym")==0){
REG=new reg_f3d_sym<PrecisionTYPE>(referenceImage->nt,floatingImage->nt);
break;
}
......@@ -568,9 +567,6 @@ int main(int argc, char **argv)
else if(strcmp(argv[i], "-nogr") ==0){
REG->NoGridRefinement();
}
else if(strcmp(argv[i], "-step")==0 || strcmp(argv[i], "--step")==0){
REG->SetCompositionStepNumber(atoi(argv[++i]));
}
else if(strcmp(argv[i], "-gce")==0 || strcmp(argv[i], "--gce")==0){
REG->UseGradientCumulativeExp();
}
......
......@@ -361,9 +361,8 @@ int main(int argc, char **argv)
#ifndef NDEBUG
printf("[NiftyReg DEBUG] Computation of the deformation field from the affine transformation\n");
#endif
reg_affine_positionField(affineTransformationMatrix,
referenceImage,
deformationFieldImage);
reg_affine_deformationField(affineTransformationMatrix,
deformationFieldImage);
}
}
......
......@@ -207,9 +207,9 @@ int main(int argc, char **argv)
}
else if(strcmp(argv[i],"-makeAff")==0 || strcmp(argv[i],"--makeAff")==0){
flag->makeAffFlag=true;
param->outputTransName=argv[++i];
for(int j=0;j<12;++j)
param->affTransParam[j]=static_cast<float>(atof(argv[++i]));
param->outputTransName=argv[++i];
}
else if(strcmp(argv[i],"-aff2rig")==0 || strcmp(argv[i],"--aff2rig")==0){
flag->aff2rigFlag=true;
......@@ -289,10 +289,8 @@ int main(int argc, char **argv)
outputTransformationImage->nvox=(size_t)outputTransformationImage->nx *
outputTransformationImage->ny * outputTransformationImage->nz *
outputTransformationImage->nt * outputTransformationImage->nu;
if(referenceImage->datatype!=NIFTI_TYPE_FLOAT64){
outputTransformationImage->nbyper=sizeof(float);
outputTransformationImage->datatype=NIFTI_TYPE_FLOAT32;
}
outputTransformationImage->nbyper=sizeof(float);
outputTransformationImage->datatype=NIFTI_TYPE_FLOAT32;
outputTransformationImage->intent_code=NIFTI_INTENT_VECTOR;
memset(outputTransformationImage->intent_name, 0, 16);
strcpy(outputTransformationImage->intent_name,"NREG_TRANS");
......@@ -349,60 +347,62 @@ int main(int argc, char **argv)
// Create a deformation field
else if(flag->outputDefFlag || flag->outputDispFlag){
if(affineTransformation!=NULL){
NULL;
reg_affine_deformationField(affineTransformation,outputTransformationImage);
}
switch(static_cast<int>(reg_round(inputTransformationImage->intent_p1))){
case DEF_FIELD:
// the current in transformation is copied
memcpy(outputTransformationImage->data,inputTransformationImage->data,
outputTransformationImage->nvox*outputTransformationImage->nbyper);
break;
case DISP_FIELD:
// the current in transformation is copied and converted
memcpy(outputTransformationImage->data,inputTransformationImage->data,
outputTransformationImage->nvox*outputTransformationImage->nbyper);
reg_getDeformationFromDisplacement(outputTransformationImage);
break;
case SPLINE_GRID:
// The output field is filled with an identity deformation field
memset(outputTransformationImage->data,
0,
outputTransformationImage->nvox*outputTransformationImage->nbyper);
reg_getDeformationFromDisplacement(outputTransformationImage);
// The spline transformation is composed with the identity field
reg_spline_getDeformationField(inputTransformationImage,
outputTransformationImage,
NULL, // no mask
true, // composition is used,
true // b-spline are used
);
break;
case DEF_VEL_FIELD:
// The flow field is exponentiated
reg_defField_getDeformationFieldFromFlowField(inputTransformationImage,
outputTransformationImage,
false // step number is not updated
);
break;
case DISP_VEL_FIELD:
// The input transformation is converted into a def flow
reg_getDeformationFromDisplacement(outputTransformationImage);
// The flow field is exponentiated
reg_defField_getDeformationFieldFromFlowField(inputTransformationImage,
outputTransformationImage,
false // step number is not updated
);
break;
case SPLINE_VEL_GRID:
// The spline parametrisation is converted into a dense flow and exponentiated
reg_spline_getDeformationFieldFromVelocityGrid(inputTransformationImage,
outputTransformationImage,
false // step number is not updated
);
break;
default:
fprintf(stderr,"[NiftyReg ERROR] Unknown input transformation type\n");
return 1;
else{
switch(static_cast<int>(reg_round(inputTransformationImage->intent_p1))){
case DEF_FIELD:
// the current in transformation is copied
memcpy(outputTransformationImage->data,inputTransformationImage->data,
outputTransformationImage->nvox*outputTransformationImage->nbyper);
break;
case DISP_FIELD:
// the current in transformation is copied and converted
memcpy(outputTransformationImage->data,inputTransformationImage->data,
outputTransformationImage->nvox*outputTransformationImage->nbyper);
reg_getDeformationFromDisplacement(outputTransformationImage);
break;
case SPLINE_GRID:
// The output field is filled with an identity deformation field
memset(outputTransformationImage->data,
0,
outputTransformationImage->nvox*outputTransformationImage->nbyper);
reg_getDeformationFromDisplacement(outputTransformationImage);
// The spline transformation is composed with the identity field
reg_spline_getDeformationField(inputTransformationImage,
outputTransformationImage,
NULL, // no mask
true, // composition is used,
true // b-spline are used
);
break;
case DEF_VEL_FIELD:
// The flow field is exponentiated
reg_defField_getDeformationFieldFromFlowField(inputTransformationImage,
outputTransformationImage,
false // step number is not updated
);
break;
case DISP_VEL_FIELD:
// The input transformation is converted into a def flow
reg_getDeformationFromDisplacement(outputTransformationImage);
// The flow field is exponentiated
reg_defField_getDeformationFieldFromFlowField(inputTransformationImage,
outputTransformationImage,
false // step number is not updated
);
break;
case SPLINE_VEL_GRID:
// The spline parametrisation is converted into a dense flow and exponentiated
reg_spline_getDeformationFieldFromVelocityGrid(inputTransformationImage,
outputTransformationImage,
false // step number is not updated
);
break;
default:
fprintf(stderr,"[NiftyReg ERROR] Unknown input transformation type\n");
return 1;
}
}
outputTransformationImage->intent_p1=DEF_FIELD;
outputTransformationImage->intent_p2=0;
......@@ -521,7 +521,7 @@ int main(int argc, char **argv)
output1TransImage->data=(void *)calloc
(output1TransImage->nvox,output1TransImage->nbyper);
if(affine1Trans!=NULL){
reg_affine_positionField(affine1Trans,referenceImage,output1TransImage);
reg_affine_deformationField(affine1Trans,output1TransImage);
}
else switch(reg_round(input1TransImage->intent_p1)){
case SPLINE_GRID:
......@@ -574,7 +574,7 @@ int main(int argc, char **argv)
output2TransImage->intent_p1=DEF_FIELD;
output2TransImage->data=(void *)calloc
(output2TransImage->nvox,output2TransImage->nbyper);
reg_affine_positionField(affine2Trans,output2TransImage,output2TransImage);
reg_affine_deformationField(affine2Trans,output2TransImage);
reg_defField_compose(output2TransImage,output1TransImage,NULL);
}
else{
......
......@@ -81,8 +81,9 @@ typedef unsigned __int64 airULLong;
#define AIR_LLONG(x) x##i64
#define AIR_ULLONG(x) x##ui64
#else
typedef signed long long airLLong;
typedef unsigned long long airULLong;
#include <stdint.h>
typedef int64_t airLLong;
typedef uint64_t airULLong;
#define AIR_LLONG_FMT "%lld"
#define AIR_ULLONG_FMT "%llu"
#define AIR_LLONG(x) x##ll
......
......@@ -426,9 +426,8 @@ void reg_aladin<T>::InitialiseBlockMatching(int CurrentPercentageOfBlockToUse)
template <class T>
void reg_aladin<T>::GetDeformationField()
{
reg_affine_positionField(this->TransformationMatrix,
this->CurrentReference,
this->deformationFieldImage);
reg_affine_deformationField(this->TransformationMatrix,
this->deformationFieldImage);
}
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
template <class T>
......
......@@ -249,9 +249,8 @@ void reg_aladin_sym<T>::SetCurrentImages()
template <class T>
void reg_aladin_sym<T>::GetBackwardDeformationField()
{
reg_affine_positionField(this->BackwardTransformationMatrix,
this->CurrentFloating,
this->BackwardDeformationFieldImage);
reg_affine_deformationField(this->BackwardTransformationMatrix,
this->BackwardDeformationFieldImage);
}
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
template <class T>
......
......@@ -764,7 +764,7 @@ template <class T>
void reg_base<T>::GetVoxelBasedGradient()
{
// The intensity gradient is first computed
if(measure_dti!=NULL){
if(this->measure_dti!=NULL){
reg_getImageGradient(this->currentFloating,
this->warpedGradientImage,
this->deformationFieldImage,
......
......@@ -237,22 +237,7 @@ void reg_f3d<T>::Initisalise()
// The control point position image is initialised with the affine transformation
if(this->affineTransformation==NULL){
mat44 identityAffine;
identityAffine.m[0][0]=1.f;
identityAffine.m[0][1]=0.f;
identityAffine.m[0][2]=0.f;
identityAffine.m[0][3]=0.f;
identityAffine.m[1][0]=0.f;
identityAffine.m[1][1]=1.f;
identityAffine.m[1][2]=0.f;
identityAffine.m[1][3]=0.f;
identityAffine.m[2][0]=0.f;
identityAffine.m[2][1]=0.f;
identityAffine.m[2][2]=1.f;
identityAffine.m[2][3]=0.f;
identityAffine.m[3][0]=0.f;
identityAffine.m[3][1]=0.f;
identityAffine.m[3][2]=0.f;
identityAffine.m[3][3]=1.f;
reg_mat44_eye(&identityAffine);
if(reg_spline_initialiseControlPointGridWithAffine(&identityAffine, this->controlPointGrid))
reg_exit(1);
}
......@@ -514,17 +499,6 @@ void reg_f3d<T>::GetSimilarityMeasureGradient()
);
}
// float currentNodeSpacing[3]={
// this->controlPointGrid->dx/this->currentReference->dx,
// this->controlPointGrid->dy/this->currentReference->dy,
// this->controlPointGrid->dz/this->currentReference->dz
// };
// reg_tools_CubicSplineKernelConvolution(this->voxelBasedMeasureGradientImage,
// currentNodeSpacing);
// reg_io_WriteImageFile(this->voxelBasedMeasureGradientImage,
// "grad_old.nii");
// reg_exit(1);
// The node based NMI gradient is extracted
reg_voxelCentric2NodeCentric(this->transformationGradient,
this->voxelBasedMeasureGradientImage,
......@@ -533,7 +507,10 @@ void reg_f3d<T>::GetSimilarityMeasureGradient()
/* The similarity measure gradient is converted from voxel space to real space */
mat44 *floatingMatrix_xyz=NULL;
size_t controlPointNumber=(size_t)this->controlPointGrid->nx*this->controlPointGrid->ny*this->controlPointGrid->nz;
size_t controlPointNumber=
(size_t)this->controlPointGrid->nx*
this->controlPointGrid->ny*
this->controlPointGrid->nz;
#ifdef _WIN32
int i;
#else
......@@ -1037,6 +1014,7 @@ void reg_f3d<T>::PrintCurrentObjFunctionValue(T currentSize)
template<class T>
void reg_f3d<T>::GetObjectiveFunctionGradient()
{
if(!this->useApproxGradient){
// Compute the gradient of the similarity measure
if(this->similarityWeight>0){
......
This diff is collapsed.
......@@ -18,14 +18,10 @@ template <class T>
class reg_f3d2 : public reg_f3d_sym<T>
{
protected:
int stepNumber;
bool BCHUpdate;
bool useGradientCumulativeExp;
int BCHUpdateValue;
mat33 *forward2backward_reorient;
mat33 *backward2forward_reorient;
virtual void DefineReorientationMatrices();
virtual void GetDeformationField();
virtual void GetInverseConsistencyErrorField();
virtual void GetInverseConsistencyGradient();
......@@ -36,7 +32,6 @@ class reg_f3d2 : public reg_f3d_sym<T>
virtual void UseGradientCumulativeExp();
public:
virtual void SetCompositionStepNumber(int);
reg_f3d2(int refTimePoint,int floTimePoint);
~reg_f3d2();
virtual void Initisalise();
......
......@@ -338,13 +338,6 @@ void reg_f3d_sym<T>::CheckParameters()
reg_f3d<T>::CheckParameters();
if(this->affineTransformation!=NULL){
fprintf(stderr, "[NiftyReg F3D_SYM ERROR] The inverse consistency parametrisation does not handle affine input\n");
fprintf(stderr, "[NiftyReg F3D_SYM ERROR] Please update your floating image sform using reg_transform\n");
fprintf(stderr, "[NiftyReg F3D_SYM ERROR] and use the updated floating image as an input\n.");
reg_exit(1);
}
// CHECK THE FLOATING MASK DIMENSION IF IT IS DEFINED
if(this->floatingMaskImage!=NULL){
if(this->inputFloating->nx != this->floatingMaskImage->nx ||
......@@ -372,7 +365,7 @@ void reg_f3d_sym<T>::CheckParameters()
this->jacobianLogWeight /= penaltySum;
this->inverseConsistencyWeight /= penaltySum;
}
else this->similarityWeight=1.0 - penaltySum;
else this->similarityWeight = 1.0 - penaltySum;
return;
}
......@@ -411,28 +404,20 @@ void reg_f3d_sym<T>::Initisalise()
this->floatingPyramid[0],
gridSpacing);
// the backward control point is initialised using an affine transformation
mat44 matrixAffine;
matrixAffine.m[0][0]=1.f;
matrixAffine.m[0][1]=0.f;
matrixAffine.m[0][2]=0.f;
matrixAffine.m[0][3]=0.f;
matrixAffine.m[1][0]=0.f;
matrixAffine.m[1][1]=1.f;
matrixAffine.m[1][2]=0.f;
matrixAffine.m[1][3]=0.f;
matrixAffine.m[2][0]=0.f;
matrixAffine.m[2][1]=0.f;
matrixAffine.m[2][2]=1.f;
matrixAffine.m[2][3]=0.f;
matrixAffine.m[3][0]=0.f;
matrixAffine.m[3][1]=0.f;
matrixAffine.m[3][2]=0.f;
matrixAffine.m[3][3]=1.f;
if(reg_spline_initialiseControlPointGridWithAffine(&matrixAffine, this->controlPointGrid))
reg_exit(1);
if(reg_spline_initialiseControlPointGridWithAffine(&matrixAffine, this->backwardControlPointGrid))
reg_exit(1);
// the backward control point grid is initialised using the inverse affine transformation
if(this->affineTransformation!=NULL){
mat44 inverseAffineTransformation=nifti_mat44_inverse(*this->affineTransformation);
if(reg_spline_initialiseControlPointGridWithAffine(&inverseAffineTransformation,
this->backwardControlPointGrid))
reg_exit(1);
}
else{
mat44 identityAffine;
reg_mat44_eye(&identityAffine);
if(reg_spline_initialiseControlPointGridWithAffine(&identityAffine, this->backwardControlPointGrid))
reg_exit(1);
}
// Set the floating mask image pyramid
if(this->usePyramid){
......@@ -661,6 +646,76 @@ double reg_f3d_sym<T>::ComputeL2NormDispPenaltyTerm()
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
template <class T>
void reg_f3d_sym<T>::GetVoxelBasedGradient()
{
// The intensity gradient is first computed
if(this->measure_dti!=NULL){
reg_getImageGradient(this->currentFloating,
this->warpedGradientImage,
this->deformationFieldImage,
this->currentMask,
this->interpolation,
this->warpedPaddingValue,
this->measure_dti->GetActiveTimepoints(),
this->forwardJacobianMatrix,
this->warped);
reg_getImageGradient(this->currentReference,
this->backwardWarpedGradientImage,
this->backwardDeformationFieldImage,
this->currentFloatingMask,
this->interpolation,
this->warpedPaddingValue,
this->measure_dti->GetActiveTimepoints(),
this->backwardJacobianMatrix,
this->backwardWarped);
}
else{
reg_getImageGradient(this->currentFloating,
this->warpedGradientImage,
this->deformationFieldImage,
this->currentMask,
this->interpolation,
this->warpedPaddingValue);
reg_getImageGradient(this->currentReference,
this->backwardWarpedGradientImage,
this->backwardDeformationFieldImage,
this->currentFloatingMask,
this->interpolation,
this->warpedPaddingValue);
}
// The voxel based gradient image is filled with zeros
reg_tools_multiplyValueToImage(this->voxelBasedMeasureGradientImage,
this->voxelBasedMeasureGradientImage,
0.f);
reg_tools_multiplyValueToImage(this->backwardVoxelBasedMeasureGradientImage,
this->backwardVoxelBasedMeasureGradientImage,
0.f);
// The gradient of the various measures of similarity are computed
if(this->measure_nmi!=NULL)
this->measure_nmi->GetVoxelBasedSimilarityMeasureGradient();
if(this->measure_multichannel_nmi!=NULL)
this->measure_multichannel_nmi->GetVoxelBasedSimilarityMeasureGradient();
if(this->measure_ssd!=NULL)
this->measure_ssd->GetVoxelBasedSimilarityMeasureGradient();
if(this->measure_kld!=NULL)
this->measure_kld->GetVoxelBasedSimilarityMeasureGradient();
if(this->measure_lncc!=NULL)
this->measure_lncc->GetVoxelBasedSimilarityMeasureGradient();
if(this->measure_dti!=NULL)
this->measure_dti->GetVoxelBasedSimilarityMeasureGradient();
return;
}
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
/* \/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/ */
template <class T>
void reg_f3d_sym<T>::GetSimilarityMeasureGradient()
{
reg_f3d<T>::GetSimilarityMeasureGradient();
......@@ -705,8 +760,8 @@ void reg_f3d_sym<T>::GetSimilarityMeasureGradient()
/* The gradient is converted from voxel space to real space */
mat44 *referenceMatrix_xyz=NULL;
size_t controlPointNumber=
(size_t)this->backwardControlPointGrid->nx *
this->backwardControlPointGrid->ny *
(size_t)this->backwardControlPointGrid->nx*
this->backwardControlPointGrid->ny*
this->backwardControlPointGrid->nz;
#ifdef _WIN32
int i;
......
......@@ -62,6 +62,7 @@ class reg_f3d_sym : public reg_f3d<T>
virtual double ComputeL2NormDispPenaltyTerm();
virtual void GetDeformationField();
virtual void WarpFloatingImage(int);
virtual void GetVoxelBasedGradient();
virtual void GetSimilarityMeasureGradient();
virtual void GetObjectiveFunctionGradient();
virtual void GetBendingEnergyGradient();
......
......@@ -18,49 +18,47 @@
/* *************************************************************** */
/* *************************************************************** */
template <class FieldTYPE>
void reg_affine_positionField2D(mat44 *affineTransformation,
nifti_image *targetImage,
nifti_image *positionFieldImage)
void reg_affine_deformationField2D(mat44 *affineTransformation,
nifti_image *deformationFieldImage)
{
FieldTYPE *positionFieldPtr = static_cast<FieldTYPE *>(positionFieldImage->data);
FieldTYPE *deformationFieldPtr = static_cast<FieldTYPE *>(deformationFieldImage->data);
size_t positionFieldXIndex=0;
size_t positionFieldYIndex=targetImage->nx*targetImage->ny;
size_t deformationFieldXIndex=0;
size_t deformationFieldYIndex=deformationFieldImage->nx*deformationFieldImage->ny;
mat44 *targetMatrix;
if(targetImage->sform_code>0){
targetMatrix=&(targetImage->sto_xyz);
if(deformationFieldImage->sform_code>0){
targetMatrix=&(deformationFieldImage->sto_xyz);
}
else targetMatrix=&(targetImage->qto_xyz);
else targetMatrix=&(deformationFieldImage->qto_xyz);
mat44 voxelToRealDeformed = reg_mat44_mul(affineTransformation, targetMatrix);
float index[3];
float position[3];
index[2]=0;
for(int y=0; y<targetImage->ny; y++){
for(int y=0; y<deformationFieldImage->ny; y++){
index[1]=(float)y;
for(int x=0; x<targetImage->nx; x++){
for(int x=0; x<deformationFieldImage->nx; x++){
index[0]=(float)x;
reg_mat44_mul(&voxelToRealDeformed, index, position);
/* the deformation field (real coordinates) is stored */
positionFieldPtr[positionFieldXIndex++] = position[0];
positionFieldPtr[positionFieldYIndex++] = position[1];
deformationFieldPtr[deformationFieldXIndex++] = position[0];
deformationFieldPtr[deformationFieldYIndex++] = position[1];
}
}
}
/* *************************************************************** */
template <class FieldTYPE>
void reg_affine_positionField3D(mat44 *affineTransformation,
nifti_image *targetImage,
nifti_image *deformationFieldImage)
void reg_affine_deformationField3D(mat44 *affineTransformation,
nifti_image *deformationFieldImage)
{
size_t voxelNumber=targetImage->nx*targetImage->ny*targetImage->nz;
FieldTYPE *positionFieldPtrX = static_cast<FieldTYPE *>(deformationFieldImage->data);
FieldTYPE *positionFieldPtrY = &positionFieldPtrX[voxelNumber];
FieldTYPE *positionFieldPtrZ = &positionFieldPtrY[voxelNumber];
size_t voxelNumber=deformationFieldImage->nx*deformationFieldImage->ny*deformationFieldImage->nz;
FieldTYPE *deformationFieldPtrX = static_cast<FieldTYPE *>(deformationFieldImage->data);
FieldTYPE *deformationFieldPtrY = &deformationFieldPtrX[voxelNumber];
FieldTYPE *deformationFieldPtrZ = &deformationFieldPtrY[voxelNumber];
mat44 *targetMatrix;
if(deformationFieldImage->sform_code>0){
......@@ -75,8 +73,8 @@ void reg_affine_positionField3D(mat44 *affineTransformation,
size_t index;
#ifdef _OPENMP
#pragma omp parallel for default(none) \
shared(deformationFieldImage, voxelToRealDeformed, positionFieldPtrX, \
positionFieldPtrY, positionFieldPtrZ) \
shared(deformationFieldImage, voxelToRealDeformed, deformationFieldPtrX, \
deformationFieldPtrY, deformationFieldPtrZ) \
private(voxel, position, x, y, z, index)
#endif
for(z=0; z<deformationFieldImage->nz; z++){
......@@ -90,42 +88,41 @@ void reg_affine_positionField3D(mat44 *affineTransformation,
reg_mat44_mul(&voxelToRealDeformed, voxel, position);
/* the deformation field (real coordinates) is stored */
positionFieldPtrX[index] = position[0];
positionFieldPtrY[index] = position[1];
positionFieldPtrZ[index] = position[2];
deformationFieldPtrX[index] = position[0];
deformationFieldPtrY[index] = position[1];
deformationFieldPtrZ[index] = position[2];
index++;
}
}
}
}
/* *************************************************************** */
void reg_affine_positionField(mat44 *affineTransformation,
nifti_image *targetImage,
nifti_image *positionFieldImage)
void reg_affine_deformationField(mat44 *affineTransformation,
nifti_image *deformationFieldImage)
{
if(targetImage->nz==1){
switch(positionFieldImage->datatype){