Commit ed7f1d00 authored by Marc Modat's avatar Marc Modat

I modified reg_jacobian to deal with any input transformation

parent 50ee490f
This diff is collapsed.
......@@ -337,6 +337,8 @@ int main(int argc, char **argv)
reg_defField_compose(inputTransformationImage,
tempFlowField,
NULL);
tempFlowField->intent_p1=inputTransformationImage->intent_p1;
tempFlowField->intent_p2=inputTransformationImage->intent_p2;
reg_defField_getDeformationFieldFromFlowField(tempFlowField,
deformationFieldImage,
false);
......@@ -344,7 +346,7 @@ int main(int argc, char **argv)
}
break;
case SPLINE_VEL_GRID:
reg_spline_getDeformationFieldFromVelocityGrid(inputTransformationImage,
reg_spline_getDefFieldFromVelocityGrid(inputTransformationImage,
deformationFieldImage,
false);
break;
......
......@@ -460,7 +460,7 @@ int main(int argc, char **argv)
printf("[NiftyReg] The specified transformation is a spline velocity parametrisation:\n[NiftyReg] %s\n",
inputTransformationImage->fname);
// The spline parametrisation is converted into a dense flow and exponentiated
reg_spline_getDeformationFieldFromVelocityGrid(inputTransformationImage,
reg_spline_getDefFieldFromVelocityGrid(inputTransformationImage,
outputTransformationImage,
false // step number is not updated
);
......@@ -672,7 +672,7 @@ int main(int argc, char **argv)
case SPLINE_VEL_GRID:
printf("[NiftyReg] Transformation 1 is a spline velocity field parametrisation:\n[NiftyReg] %s\n",
input1TransImage->fname);
reg_spline_getDeformationFieldFromVelocityGrid(input1TransImage,
reg_spline_getDefFieldFromVelocityGrid(input1TransImage,
output1TransImage,
false // the number of step is not automatically updated
);
......@@ -763,7 +763,7 @@ int main(int argc, char **argv)
(output2TransImage->nvox,output2TransImage->nbyper);
printf("[NiftyReg] Transformation 2 is a spline velocity field parametrisation:\n[NiftyReg] %s\n",
input2TransImage->fname);
reg_spline_getDeformationFieldFromVelocityGrid(input2TransImage,
reg_spline_getDefFieldFromVelocityGrid(input2TransImage,
output2TransImage,
false // the number of step is not automatically updated
);
......
......@@ -99,7 +99,7 @@ void reg_f3d2<T>::GetDeformationField()
printf("[NiftyReg DEBUG] Velocity integration forward. Step number update=%i\n",updateStepNumber);
#endif
// The forward transformation is computed using the scaling-and-squaring approach
reg_spline_getDeformationFieldFromVelocityGrid(this->controlPointGrid,
reg_spline_getDefFieldFromVelocityGrid(this->controlPointGrid,
this->deformationFieldImage,
updateStepNumber
);
......@@ -109,7 +109,7 @@ void reg_f3d2<T>::GetDeformationField()
// The number of step number is copied over from the forward transformation
this->backwardControlPointGrid->intent_p2=this->controlPointGrid->intent_p2;
// The backward transformation is computed using the scaling-and-squaring approach
reg_spline_getDeformationFieldFromVelocityGrid(this->backwardControlPointGrid,
reg_spline_getDefFieldFromVelocityGrid(this->backwardControlPointGrid,
this->backwardDeformationFieldImage,
false
);
......
This diff is collapsed.
......@@ -48,11 +48,11 @@ void reg_createControlPointGrid(nifti_image **controlPointGridImage,
extern "C++" template <class DTYPE>
void reg_createSymmetricControlPointGrids(nifti_image **forwardGridImage,
nifti_image **backwardGridImage,
nifti_image *referenceImage,
nifti_image *floatingImage,
mat44 *forwardAffineTrans,
float *spacing);
nifti_image **backwardGridImage,
nifti_image *referenceImage,
nifti_image *floatingImage,
mat44 *forwardAffineTrans,
float *spacing);
/* *************************************************************** */
/** @brief Compute a dense deformation field in the space of a reference
* image from a grid of control point.
......@@ -73,7 +73,7 @@ void reg_spline_getDeformationField(nifti_image *controlPointGridImage,
int *mask,
bool composition,
bool bspline
);
);
/* *************************************************************** */
/** @brief Compute and return the average bending energy computed using cubic b-spline.
* The value is approximated as the bending energy is computated at
......@@ -97,9 +97,9 @@ double reg_spline_approxBendingEnergy(nifti_image *controlPointGridImage);
*/
extern "C++"
void reg_spline_approxBendingEnergyGradient(nifti_image *controlPointGridImage,
nifti_image *gradientImage,
float weight
);
nifti_image *gradientImage,
float weight
);
/* *************************************************************** */
/** @brief Compute and return the linear elastic energy terms approximated
* at the control point positions only.
......@@ -112,7 +112,7 @@ void reg_spline_approxBendingEnergyGradient(nifti_image *controlPointGridImage,
extern "C++"
void reg_spline_linearEnergy(nifti_image *controlPointGridImage,
double *values
);
);
/* *************************************************************** */
/** @brief Compute the gradient of the linear elastic energy terms
* approximated at the control point positions only.
......@@ -133,7 +133,7 @@ void reg_spline_linearEnergyGradient(nifti_image *controlPointGridImage,
nifti_image *gradientImage,
float weight0,
float weight1
);
);
/* *************************************************************** */
/** @brief Compute and return the L2 norm of the displacement approximated
* at the control point positions only.
......@@ -170,7 +170,7 @@ void reg_spline_L2norm_dispGradient(nifti_image *controlPointGridImage,
extern "C++"
void reg_spline_GetJacobianMap(nifti_image *controlPointGridImage,
nifti_image *jacobianImage
);
);
/* *************************************************************** */
/** @brief Compute the average Jacobian determinant
* @param controlPointGridImage Image that contains the transformation
......@@ -183,10 +183,10 @@ void reg_spline_GetJacobianMap(nifti_image *controlPointGridImage,
*/
extern "C++"
double reg_spline_getJacobianPenaltyTerm(nifti_image *controlPointGridImage,
nifti_image *referenceImage,
bool approx,
bool useHeaderInformation=false
);
nifti_image *referenceImage,
bool approx,
bool useHeaderInformation=false
);
/* *************************************************************** */
/** @brief Compute the gradient at every control point position of the
* Jacobian determinant based penalty term
......@@ -206,12 +206,12 @@ double reg_spline_getJacobianPenaltyTerm(nifti_image *controlPointGridImage,
*/
extern "C++"
void reg_spline_getJacobianPenaltyTermGradient(nifti_image *controlPointGridImage,
nifti_image *referenceImage,
nifti_image *gradientImage,
float weight,
bool approx,
bool useHeaderInformation=false
);
nifti_image *referenceImage,
nifti_image *gradientImage,
float weight,
bool approx,
bool useHeaderInformation=false
);
/* *************************************************************** */
/** @brief Compute the Jacobian matrix at every voxel position
* using a cubic b-spline parametrisation. This function does require
......@@ -227,7 +227,7 @@ extern "C++"
void reg_spline_GetJacobianMatrix(nifti_image *referenceImage,
nifti_image *controlPointGridImage,
mat33 *jacobianImage
);
);
/* *************************************************************** */
/** @brief Correct the folding in the transformation parametrised through
* cubic B-Spline
......@@ -241,7 +241,7 @@ extern "C++"
double reg_spline_correctFolding(nifti_image *controlPointGridImage,
nifti_image *referenceImage,
bool approx
);
);
/* *************************************************************** */
/** @brief Upsample an image from voxel space to node space using
* millimiter correspendences.
......@@ -261,7 +261,7 @@ void reg_voxelCentric2NodeCentric(nifti_image *nodeImage,
float weight,
bool update,
mat44 *voxelToMillimeter = NULL
);
);
/* *************************************************************** */
/** @brief Refine a grid of control points
* @param referenceImage Image that defined the space of the reference
......@@ -272,7 +272,7 @@ void reg_voxelCentric2NodeCentric(nifti_image *nodeImage,
extern "C++"
void reg_spline_refineControlPointGrid(nifti_image *controlPointGridImage,
nifti_image *referenceImage = NULL
);
);
/* *************************************************************** */
/** @brief This function compose the a first control point image with a second one:
* Grid2(x) <= Grid1(Grid2(x)).
......@@ -294,7 +294,7 @@ int reg_spline_cppComposition(nifti_image *grid1,
bool displacement1,
bool displacement2,
bool bspline
);
);
/* *************************************************************** */
......@@ -358,10 +358,10 @@ void reg_defFieldInvert(nifti_image *inputDeformationField,
/* *************************************************************** */
extern "C++"
void reg_defField_getDeformationFieldFromFlowField(nifti_image *flowFieldImage,
nifti_image *deformationFieldImage,
bool updateStepNumber);
nifti_image *deformationFieldImage,
bool updateStepNumber);
/* *********************************************** */
/* **** VELOCITY FIELD BASED FUNCTIONS **** */
/* **** FLOW BASED FUNCTIONS **** */
/* *********************************************** */
/* *************************************************************** */
......@@ -375,42 +375,48 @@ void reg_defField_getDeformationFieldFromFlowField(nifti_image *flowFieldImage,
* the Jacobian matrices of the transformation
*/
extern "C++"
int reg_spline_GetJacobianMatricesFromVelocityField(nifti_image* referenceImage,
nifti_image* velocityFieldImage,
mat33* jacobianMatrices
);
int reg_defField_GetJacobianMatFromFlowField(mat33* jacobianMatrices,
nifti_image* flowFieldImage);
extern "C++"
int reg_spline_GetJacobianMatFromVelocityGrid(mat33* jacobianMatrices,
nifti_image* velocityGridImage,
nifti_image* referenceImage
);
/* *************************************************************** */
/** @brief This function computed a Jacobian determinant map by integrating
* the velocity field
* the velocity grid
* @param jacobianDetImage This image will be filled with the Jacobian
* determinants of every voxel.
* @param velocityFieldImage Image that contains a velocity field
* parametrised using a grid of control points
*/
extern "C++"
int reg_spline_GetJacobianDetFromVelocityField(nifti_image* jacobianDetImage,
nifti_image* velocityFieldImage
);
int reg_defField_GetJacobianDetFromFlowField(nifti_image* jacobianDetImage,
nifti_image* flowFieldImage
);
extern "C++"
int reg_spline_GetJacobianDetFromVelocityGrid(nifti_image* jacobianDetImage,
nifti_image* velocityGridImage);
/* *************************************************************** */
/** @brief The deformation field (img2) is computed by integrating
* a velocity field (img1)
* a velocity Grid (img1)
* @param velocityFieldImage Image that contains a velocity field
* parametrised using a grid of control points
* @param deformationFieldImage Deformation field image that will
* be filled using the exponentiation of the velocity field.
*/
extern "C++"
void reg_spline_getDeformationFieldFromVelocityGrid(nifti_image *velocityFieldGrid,
nifti_image *deformationFieldImage,
bool updateStepNumber);
void reg_spline_getDefFieldFromVelocityGrid(nifti_image *velocityFieldGrid,
nifti_image *deformationFieldImage,
bool updateStepNumber);
/* *************************************************************** */
extern "C++"
void reg_spline_getIntermediateDefFieldFromVelGrid(nifti_image *velocityFieldGrid,
nifti_image **deformationFieldImage);
nifti_image **deformationFieldImage);
/* *************************************************************** */
extern "C++"
void reg_spline_getFlowFieldFromVelocityGrid(nifti_image *velocityFieldGrid,
nifti_image *flowField);
nifti_image *flowField);
/* *************************************************************** */
......
This diff is collapsed.
......@@ -61,9 +61,9 @@ int main(int argc, char **argv)
reg_getDeformationFromDisplacement(controlPointGrid);
size_t controlPointNumber =
controlPointGrid->nx *
controlPointGrid->ny *
controlPointGrid->nz;
controlPointGrid->nx *
controlPointGrid->ny *
controlPointGrid->nz;
// Create a mask
int *mask = (int *)malloc(voxelNumber*sizeof(int));
......@@ -73,14 +73,14 @@ int main(int argc, char **argv)
// Create the displacement field on the device
float4 *field_gpu=NULL;
cudaCommon_allocateArrayToDevice<float4>(&field_gpu,
voxelNumber);
voxelNumber);
// Create the control point grid on the device
float4 *controlPointGrid_gpu=NULL;
cudaCommon_allocateArrayToDevice<float4>(&controlPointGrid_gpu,
controlPointGrid->nx *
controlPointGrid->ny *
controlPointGrid->nz);
controlPointGrid->nx *
controlPointGrid->ny *
controlPointGrid->nz);
// Transfer the ccp coefficients from the host to the device
cudaCommon_transferNiftiToArrayOnDevice(&controlPointGrid_gpu, controlPointGrid);
......@@ -186,14 +186,14 @@ int main(int argc, char **argv)
memset(controlPointGrid->intent_name, 0, 16);
strcpy(controlPointGrid->intent_name,"NREG_VEL_STEP");
// Exponentiate the velocity field using the host
reg_spline_getDeformationFieldFromVelocityGrid(controlPointGrid,
field,
false);
reg_spline_getDefFieldFromVelocityGrid(controlPointGrid,
field,
false);
// Exponentiate the velocity field using the device
reg_getDeformationFieldFromVelocityGrid_gpu(controlPointGrid,
field,
&controlPointGrid_gpu,
&field_gpu);
field,
&controlPointGrid_gpu,
&field_gpu);
// Transfer the device field on the host
cudaCommon_transferFromDeviceToNifti<float4>(field2, &field_gpu);
// Compute the difference between both fields
......@@ -212,7 +212,7 @@ int main(int argc, char **argv)
{
double be_cpu = reg_spline_approxBendingEnergy(controlPointGrid);
double be_gpu = reg_spline_approxBendingEnergy_gpu(controlPointGrid,
&controlPointGrid_gpu);
&controlPointGrid_gpu);
maxDifferenceBE = (be_cpu / be_gpu) - 1.0;
#ifndef NDEBUG
printf("[NiftyReg DEBUG] [dim=%i] Bending energy difference: %g [=(%g/%g)-1]\n",
......@@ -232,7 +232,7 @@ int main(int argc, char **argv)
// Allocate an extra cuda array to store the device gradient
float4 *be_grad_device=NULL;
cudaCommon_allocateArrayToDevice<float4>(&be_grad_device,
controlPointNumber);
controlPointNumber);
// Set the gradients arrays to zero
reg_tools_multiplyValueToImage(be_grad_cpu,be_grad_cpu,0.f);
cudaCommon_transferNiftiToArrayOnDevice<float4>(&be_grad_device,be_grad_cpu);
......@@ -242,9 +242,9 @@ int main(int argc, char **argv)
controlPointNumber); // weight
// Compute the gradient on the device
reg_spline_approxBendingEnergyGradient_gpu(controlPointGrid,
&controlPointGrid_gpu,
&be_grad_device,
controlPointNumber); // weight
&controlPointGrid_gpu,
&be_grad_device,
controlPointNumber); // weight
// Transfer the device field on the host
cudaCommon_transferFromDeviceToNifti<float4>(be_grad_gpu, &be_grad_device);
// Compute the difference between both gradient arrays
......@@ -273,12 +273,12 @@ int main(int argc, char **argv)
bool approximation=false;
if(strcmp(type,"ajac")==0) approximation=true;
double jac_cpu = reg_spline_getJacobianPenaltyTerm(controlPointGrid,
field,
approximation); // aprox
field,
approximation); // aprox
double jac_gpu = reg_spline_getJacobianPenaltyTerm_gpu(field,
controlPointGrid,
&controlPointGrid_gpu,
approximation);
controlPointGrid,
&controlPointGrid_gpu,
approximation);
if(jac_cpu!=jac_cpu)
{
fprintf(stderr,
......@@ -321,23 +321,23 @@ int main(int argc, char **argv)
// Allocate an extra cuda array to store the device gradient
float4 *jac_grad_device=NULL;
cudaCommon_allocateArrayToDevice<float4>(&jac_grad_device,
controlPointNumber);
controlPointNumber);
// Set the gradients arrays to zero
reg_tools_multiplyValueToImage(jac_grad_cpu,jac_grad_cpu,0.f);
cudaCommon_transferNiftiToArrayOnDevice<float4>(&jac_grad_device,jac_grad_cpu);
// Compute the gradient on the host
reg_spline_getJacobianPenaltyTermGradient(controlPointGrid,
field,
jac_grad_cpu,
controlPointNumber,
approximation);
field,
jac_grad_cpu,
controlPointNumber,
approximation);
// Compute the gradient on the device
reg_spline_getJacobianPenaltyTermGradient_gpu(field,
controlPointGrid,
&controlPointGrid_gpu,
&jac_grad_device,
controlPointNumber,
approximation);
controlPointGrid,
&controlPointGrid_gpu,
&jac_grad_device,
controlPointNumber,
approximation);
// Transfer the device field on the host
cudaCommon_transferFromDeviceToNifti<float4>(jac_grad_gpu, &jac_grad_device);
// Compute the difference between both gradient arrays
......
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