Skip to content

Commit

Permalink
Refactorisations
Browse files Browse the repository at this point in the history
  • Loading branch information
onurulgen committed Sep 8, 2023
1 parent ef4f55b commit 327d516
Show file tree
Hide file tree
Showing 29 changed files with 132 additions and 274 deletions.
2 changes: 1 addition & 1 deletion niftyreg_build_version.txt
Original file line number Diff line number Diff line change
@@ -1 +1 @@
322
323
16 changes: 8 additions & 8 deletions reg-io/_reg_ReadWriteMatrix.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ void reg_tool_ReadAffineFile(mat44 *mat,
}
affineFile.close();

NR_MAT44(*mat, "Read affine transformation");
NR_MAT44_DEBUG(*mat, "Read affine transformation");

if (flirtFile) {
mat44 absoluteReference;
Expand Down Expand Up @@ -61,11 +61,11 @@ void reg_tool_ReadAffineFile(mat44 *mat,
absoluteReference.m[3][3] = absoluteFloating.m[3][3] = 1.0;

NR_DEBUG("An flirt affine file is assumed and is converted to a real word affine matrix");
NR_MAT44(*mat, "Matrix read from the input file");
NR_MAT44(*referenceMatrix, "Reference Matrix");
NR_MAT44(*floatingMatrix, "Floating Matrix");
NR_MAT44(absoluteReference, "Reference absolute Matrix");
NR_MAT44(absoluteFloating, "Floating absolute Matrix");
NR_MAT44_DEBUG(*mat, "Matrix read from the input file");
NR_MAT44_DEBUG(*referenceMatrix, "Reference Matrix");
NR_MAT44_DEBUG(*floatingMatrix, "Floating Matrix");
NR_MAT44_DEBUG(absoluteReference, "Reference absolute Matrix");
NR_MAT44_DEBUG(absoluteFloating, "Floating absolute Matrix");

absoluteFloating = nifti_mat44_inverse(absoluteFloating);
*mat = nifti_mat44_inverse(*mat);
Expand All @@ -77,7 +77,7 @@ void reg_tool_ReadAffineFile(mat44 *mat,
*mat = reg_mat44_mul(mat, &tmp);
}

NR_MAT44(*mat, "Affine matrix");
NR_MAT44_DEBUG(*mat, "Affine matrix");
}
/* *************************************************************** */
void reg_tool_ReadAffineFile(mat44 *mat, char *fileName) {
Expand Down Expand Up @@ -223,7 +223,7 @@ mat44* reg_tool_ReadMat44File(char *fileName) {
}
matrixFile.close();

NR_MAT44(*mat, "mat44 matrix");
NR_MAT44_DEBUG(*mat, "mat44 matrix");

return mat;
}
Expand Down
9 changes: 2 additions & 7 deletions reg-io/_reg_ReadWriteMatrix.h
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,6 @@
* @param flirtFile If this flag is set to true the matrix is converted
* from a Flirt (FSL) parametrisation to a standard parametrisation
*/
extern "C++"
void reg_tool_ReadAffineFile(mat44 *mat,
nifti_image *referenceImage,
nifti_image *floatingImage,
Expand All @@ -44,7 +43,6 @@ void reg_tool_ReadAffineFile(mat44 *mat,
* @param mat structure that store the affine transformation matrix
* @param filename Filename of the text file that contains the matrix to read
**/
extern "C++"
void reg_tool_ReadAffineFile(mat44 *mat,
char *filename);

Expand All @@ -54,14 +52,12 @@ void reg_tool_ReadAffineFile(mat44 *mat,
* @param filename Filename of the text file that contains the matrix to read
* @return mat44 structure that store the matrix
**/
extern "C++"
mat44* reg_tool_ReadMat44File(char *fileName);

/** @brief This function save a 4-by-4 matrix to the disk as a text file
* @param mat Matrix to be saved on the disk
* @param filename Name of the text file to save on the disk
*/
extern "C++"
void reg_tool_WriteAffineFile(const mat44 *mat,
const char *fileName);

Expand All @@ -70,7 +66,6 @@ void reg_tool_WriteAffineFile(const mat44 *mat,
* @param filename Filename of the text file that contains the matrix to read
* @return pair of values that contains the matrix size
**/
extern "C++"
std::pair<size_t, size_t> reg_tool_sizeInputMatrixFile(char *filename);
/**
* @brief Read a file that contains a m-by-n matrix and store it into
Expand All @@ -80,7 +75,7 @@ std::pair<size_t, size_t> reg_tool_sizeInputMatrixFile(char *filename);
* @param nbColumn number of column of the input matrix
* @return a pointer to a 2D array that points the read matrix
**/
extern "C++" template <class T>
template <class T>
T** reg_tool_ReadMatrixFile(char *filename,
size_t nbLine,
size_t nbColumn);
Expand All @@ -92,7 +87,7 @@ T** reg_tool_ReadMatrixFile(char *filename,
* @param nbLine number of line of the input matrix
* @param nbColumn number of column of the input matrix
**/
extern "C++" template <class T>
template <class T>
void reg_tool_WriteMatrixFile(char *filename,
T **mat,
size_t nbLine,
Expand Down
1 change: 0 additions & 1 deletion reg-lib/Content.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,6 @@ void Content::AllocateDeformationField(size_t bytes) {
deformationField->scl_slope = 1;
deformationField->scl_inter = 0;
deformationField->data = calloc(deformationField->nvox, deformationField->nbyper);
reg_tools_multiplyValueToImage(deformationField, deformationField, 0.f);
// Convert to an identity deformation field
reg_getDeformationFromDisplacement(deformationField);
}
Expand Down
5 changes: 2 additions & 3 deletions reg-lib/Debug.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -77,10 +77,9 @@ inline std::string StripFunctionName(const std::string& funcName) {
#else
#define NR_MAT33(mat, title) reg_mat33_disp(mat, title)
#define NR_MAT44(mat, title) reg_mat44_disp(mat, title)
#define NR_MAT33_DEBUG(mat, title)
#define NR_MAT44_DEBUG(mat, title)
#define NR_MAT33_DEBUG(mat, title)
#define NR_MAT44_DEBUG(mat, title)
#define NR_MAT33_VERBOSE(mat, title) if (this->verbose) NR_MAT33(mat, "[NiftyReg INFO] "s + (title))
#define NR_MAT44_VERBOSE(mat, title) if (this->verbose) NR_MAT44(mat, "[NiftyReg INFO] "s + (title))

#endif
/* *************************************************************** */
10 changes: 5 additions & 5 deletions reg-lib/_reg_aladin.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -329,7 +329,7 @@ template<class T>
void reg_aladin<T>::UpdateTransformationMatrix(int type) {
this->blockMatchingKernel->template castTo<BlockMatchingKernel>()->Calculate();
this->ltsKernel->template castTo<LtsKernel>()->Calculate(type);
NR_MAT44(*this->affineTransformation, "The updated forward matrix");
NR_MAT44_DEBUG(*this->affineTransformation, "The updated forward matrix");
}
/* *************************************************************** */
template<class T>
Expand Down Expand Up @@ -381,11 +381,11 @@ void reg_aladin<T>::Run() {
this->DebugPrintLevelInfoStart();

if (this->con->Content::GetReference()->sform_code > 0)
NR_MAT44(this->con->Content::GetReference()->sto_xyz, "Reference image matrix (sform sto_xyz)");
else NR_MAT44(this->con->Content::GetReference()->qto_xyz, "Reference image matrix (qform qto_xyz)");
NR_MAT44_DEBUG(this->con->Content::GetReference()->sto_xyz, "Reference image matrix (sform sto_xyz)");
else NR_MAT44_DEBUG(this->con->Content::GetReference()->qto_xyz, "Reference image matrix (qform qto_xyz)");
if (this->con->Content::GetFloating()->sform_code > 0)
NR_MAT44(this->con->Content::GetFloating()->sto_xyz, "Floating image matrix (sform sto_xyz)");
else NR_MAT44(this->con->Content::GetFloating()->qto_xyz, "Floating image matrix (qform qto_xyz)");
NR_MAT44_DEBUG(this->con->Content::GetFloating()->sto_xyz, "Floating image matrix (sform sto_xyz)");
else NR_MAT44_DEBUG(this->con->Content::GetFloating()->qto_xyz, "Floating image matrix (qform qto_xyz)");

/* ****************** */
/* Rigid registration */
Expand Down
3 changes: 1 addition & 2 deletions reg-lib/_reg_f3d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -181,9 +181,8 @@ void reg_f3d<T>::Initialise() {
reg_createControlPointGrid<T>(controlPointGrid, this->referencePyramid[0], gridSpacing);

// The control point grid is updated with an identity transformation
if (this->affineTransformation) {
if (this->affineTransformation)
reg_affine_getDeformationField(this->affineTransformation.get(), controlPointGrid);
}
} else {
// The control point grid image is initialised with the provided grid
controlPointGrid = inputControlPointGrid;
Expand Down
2 changes: 0 additions & 2 deletions reg-lib/cpu/_reg_blockMatching.h
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,6 @@ struct _reg_blockMatchingParam
* image to consider for the registration
* @param runningOnGPU Has to be set to true if the registration has to be performed on the GPU
*/
extern "C++"
void initialise_block_matching_method(nifti_image * referenceImage,
_reg_blockMatchingParam *params,
int percentToKeep_block,
Expand All @@ -104,7 +103,6 @@ void initialise_block_matching_method(nifti_image * referenceImage,
* relevant information
* @param mask Mask array where only voxel defined as active are considered
*/
extern "C++"
void block_matching_method(nifti_image * referenceImage,
nifti_image * warpedImage,
_reg_blockMatchingParam *params,
Expand Down
4 changes: 2 additions & 2 deletions reg-lib/cpu/_reg_dti.h
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ class reg_dti: public reg_measure {
* should be considered. If set to nullptr, all voxels are considered
* @return Returns an L2 measure of the distance between the anisotropic components of the diffusion tensors
*/
extern "C++" template <class DataType>
template <class DataType>
double reg_getDtiMeasureValue(const nifti_image *referenceImage,
const nifti_image *warpedImage,
const int *mask,
Expand All @@ -74,7 +74,7 @@ double reg_getDtiMeasureValue(const nifti_image *referenceImage,
* @param mask Array that contains a mask to specify which voxel
* should be considered. If set to nullptr, all voxels are considered
*/
extern "C++" template <class DataType>
template <class DataType>
void reg_getVoxelBasedDtiMeasureGradient(nifti_image *referenceImage,
nifti_image *warpedImage,
nifti_image *warpedGradient,
Expand Down
1 change: 0 additions & 1 deletion reg-lib/cpu/_reg_globalTrans.h
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,6 @@ typedef struct _reg_sorted_point2D _reg_sorted_point2D;
* @param deformationField Image that contains the deformation field
* that is being updated
*/
extern "C++"
void reg_affine_getDeformationField(mat44 *affine,
nifti_image *deformationField,
bool compose=false,
Expand Down
61 changes: 25 additions & 36 deletions reg-lib/cpu/_reg_localTrans.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,15 +35,13 @@ void reg_createControlPointGrid(NiftiImage& controlPointGridImage,
controlPointGridImage->cal_min = 0;
controlPointGridImage->cal_max = 0;
controlPointGridImage->pixdim[0] = 1.0f;
controlPointGridImage->pixdim[1] = controlPointGridImage->dx = spacing[0];
controlPointGridImage->pixdim[2] = controlPointGridImage->dy = spacing[1];
if (referenceImage->nz == 1) {
controlPointGridImage->pixdim[3] = controlPointGridImage->dz = 1.0f;
} else controlPointGridImage->pixdim[3] = controlPointGridImage->dz = spacing[2];
controlPointGridImage->pixdim[4] = controlPointGridImage->dt = 1.0f;
controlPointGridImage->pixdim[5] = controlPointGridImage->du = 1.0f;
controlPointGridImage->pixdim[6] = controlPointGridImage->dv = 1.0f;
controlPointGridImage->pixdim[7] = controlPointGridImage->dw = 1.0f;
controlPointGridImage.setPixDim(NiftiDim::X, spacing[0]);
controlPointGridImage.setPixDim(NiftiDim::Y, spacing[1]);
controlPointGridImage.setPixDim(NiftiDim::Z, referenceImage->nz > 1 ? spacing[2] : 1.0f);
controlPointGridImage.setPixDim(NiftiDim::T, 1.0f);
controlPointGridImage.setPixDim(NiftiDim::U, 1.0f);
controlPointGridImage.setPixDim(NiftiDim::V, 1.0f);
controlPointGridImage.setPixDim(NiftiDim::W, 1.0f);

// Reproduce the orientation of the reference image and add a one voxel shift
if (referenceImage->qform_code + referenceImage->sform_code > 0) {
Expand Down Expand Up @@ -80,7 +78,7 @@ void reg_createControlPointGrid(NiftiImage& controlPointGridImage,
originIndex[1] = -1.0f;
originIndex[2] = 0.0f;
if (referenceImage->nz > 1) originIndex[2] = -1.0f;
reg_mat44_mul(&(controlPointGridImage->qto_xyz), originIndex, originReal);
reg_mat44_mul(&controlPointGridImage->qto_xyz, originIndex, originReal);
controlPointGridImage->qto_xyz.m[0][3] = controlPointGridImage->qoffset_x = originReal[0];
controlPointGridImage->qto_xyz.m[1][3] = controlPointGridImage->qoffset_y = originReal[1];
controlPointGridImage->qto_xyz.m[2][3] = controlPointGridImage->qoffset_z = originReal[2];
Expand Down Expand Up @@ -112,19 +110,17 @@ void reg_createControlPointGrid(NiftiImage& controlPointGridImage,
controlPointGridImage->sto_xyz.m[3][3] = referenceImage->sto_xyz.m[3][3];

// Origin is shifted from 1 control point in the sform
reg_mat44_mul(&(controlPointGridImage->sto_xyz), originIndex, originReal);
reg_mat44_mul(&controlPointGridImage->sto_xyz, originIndex, originReal);
controlPointGridImage->sto_xyz.m[0][3] = originReal[0];
controlPointGridImage->sto_xyz.m[1][3] = originReal[1];
controlPointGridImage->sto_xyz.m[2][3] = originReal[2];
controlPointGridImage->sto_ijk = nifti_mat44_inverse(controlPointGridImage->sto_xyz);
}

// The grid is initialised with an identity transformation
reg_tools_multiplyValueToImage(controlPointGridImage, controlPointGridImage, 0.f);
reg_getDeformationFromDisplacement(controlPointGridImage);
controlPointGridImage->intent_code = NIFTI_INTENT_VECTOR;
memset(controlPointGridImage->intent_name, 0, 16);
strcpy(controlPointGridImage->intent_name, "NREG_TRANS");
controlPointGridImage.setIntentName("NREG_TRANS"s);
controlPointGridImage->intent_p1 = CUB_SPLINE_GRID;
}
template void reg_createControlPointGrid<float>(NiftiImage&, const NiftiImage&, const float*);
Expand All @@ -142,12 +138,12 @@ void reg_createSymmetricControlPointGrids(NiftiImage& forwardGridImage,
mat44 referenceImageSpace = referenceImage->qto_xyz;
if (referenceImage->sform_code > 0)
referenceImageSpace = referenceImage->sto_xyz;
NR_MAT44(referenceImageSpace, "Input reference image orientation");
NR_MAT44_DEBUG(referenceImageSpace, "Input reference image orientation");
// // Get the floating image space
mat44 floatingImageSpace = floatingImage->qto_xyz;
if (floatingImage->sform_code > 0)
floatingImageSpace = floatingImage->sto_xyz;
NR_MAT44(floatingImageSpace, "Input floating image orientation");
NR_MAT44_DEBUG(floatingImageSpace, "Input floating image orientation");
// Check if an affine transformation is specified
mat44 halfForwardAffine, halfBackwardAffine;
if (forwardAffineTrans != nullptr) {
Expand Down Expand Up @@ -290,10 +286,12 @@ void reg_createSymmetricControlPointGrids(NiftiImage& forwardGridImage,
backwardGridImage = NiftiImage(dims, sizeof(DataType) == sizeof(float) ? NIFTI_TYPE_FLOAT32 : NIFTI_TYPE_FLOAT64);

// Set the control point grid spacing
forwardGridImage->pixdim[1] = forwardGridImage->dx = backwardGridImage->pixdim[1] = backwardGridImage->dx = spacing[0];
forwardGridImage->pixdim[2] = forwardGridImage->dy = backwardGridImage->pixdim[2] = backwardGridImage->dy = spacing[1];
if (referenceImage->nz > 1)
forwardGridImage->pixdim[3] = forwardGridImage->dz = backwardGridImage->pixdim[3] = backwardGridImage->dz = spacing[2];
forwardGridImage.setPixDim(NiftiDim::X, spacing[0]);
backwardGridImage.setPixDim(NiftiDim::X, spacing[0]);
forwardGridImage.setPixDim(NiftiDim::Y, spacing[1]);
backwardGridImage.setPixDim(NiftiDim::Y, spacing[1]);
forwardGridImage.setPixDim(NiftiDim::Z, referenceImage->nz > 1 ? spacing[2] : 1.0f);
backwardGridImage.setPixDim(NiftiDim::Z, referenceImage->nz > 1 ? spacing[2] : 1.0f);
// Set the control point grid image orientation
forwardGridImage->qform_code = backwardGridImage->qform_code = 0;
forwardGridImage->sform_code = backwardGridImage->sform_code = 1;
Expand All @@ -313,10 +311,8 @@ void reg_createSymmetricControlPointGrids(NiftiImage& forwardGridImage,
forwardGridImage->sto_ijk = backwardGridImage->sto_ijk = nifti_mat44_inverse(forwardGridImage->sto_xyz);
// Set the intent type
forwardGridImage->intent_code = backwardGridImage->intent_code = NIFTI_INTENT_VECTOR;
memset(forwardGridImage->intent_name, 0, 16);
memset(backwardGridImage->intent_name, 0, 16);
strcpy(forwardGridImage->intent_name, "NREG_TRANS");
strcpy(backwardGridImage->intent_name, "NREG_TRANS");
forwardGridImage.setIntentName("NREG_TRANS"s);
backwardGridImage.setIntentName("NREG_TRANS"s);
forwardGridImage->intent_p1 = backwardGridImage->intent_p1 = CUB_SPLINE_GRID;
// Set the affine matrices
mat44 identity;
Expand All @@ -339,7 +335,7 @@ void reg_createSymmetricControlPointGrids(NiftiImage& forwardGridImage,
forwardGridImage->ext_list[1].edata = (char*)calloc(forwardGridImage->ext_list[1].esize - 8, sizeof(float));
memcpy(forwardGridImage->ext_list[0].edata, &halfForwardAffine, sizeof(mat44));
memcpy(forwardGridImage->ext_list[1].edata, &halfForwardAffine, sizeof(mat44));
NR_MAT44(halfForwardAffine, "Forward transformation half-affine");
NR_MAT44_DEBUG(halfForwardAffine, "Forward transformation half-affine");
// Create extensions to store the affine parametrisations for the backward transformation
backwardGridImage->num_ext = 2;
backwardGridImage->ext_list = (nifti1_extension*)malloc(2 * sizeof(nifti1_extension));
Expand All @@ -351,23 +347,20 @@ void reg_createSymmetricControlPointGrids(NiftiImage& forwardGridImage,
backwardGridImage->ext_list[1].edata = (char*)calloc(backwardGridImage->ext_list[1].esize - 8, sizeof(float));
memcpy(backwardGridImage->ext_list[0].edata, &halfBackwardAffine, sizeof(mat44));
memcpy(backwardGridImage->ext_list[1].edata, &halfBackwardAffine, sizeof(mat44));
NR_MAT44(halfBackwardAffine, "Backward transformation half-affine");
NR_MAT44_DEBUG(halfBackwardAffine, "Backward transformation half-affine");
}
// Initialise the grid with identity transformations
reg_tools_multiplyValueToImage(forwardGridImage, forwardGridImage, 0.f);
reg_tools_multiplyValueToImage(backwardGridImage, backwardGridImage, 0.f);
// Convert the parametrisations into deformation fields
reg_getDeformationFromDisplacement(forwardGridImage);
reg_getDeformationFromDisplacement(backwardGridImage);
}
template void reg_createSymmetricControlPointGrids<float>(NiftiImage&, NiftiImage&, const NiftiImage&, const NiftiImage&, const mat44*, const float*);
template void reg_createSymmetricControlPointGrids<double>(NiftiImage&, NiftiImage&, const NiftiImage&, const NiftiImage&, const mat44*, const float*);
/* *************************************************************** */
extern "C++" template <class DataType>
template <class DataType>
void reg_createDeformationField(NiftiImage& deformationFieldImage,
const nifti_image *referenceImage) {
// The header information from the reference image are copied over
deformationFieldImage = nifti_copy_nim_info(referenceImage);
deformationFieldImage = NiftiImage(const_cast<nifti_image*>(referenceImage), NiftiImage::Copy::ImageInfo);
// The dimension are updated to store the deformation vector along U index
// in a 5D image
deformationFieldImage.setDim(NiftiDim::NDim, 5);
Expand All @@ -390,10 +383,8 @@ void reg_createDeformationField(NiftiImage& deformationFieldImage,
deformationFieldImage->scl_slope = 1;
deformationFieldImage->scl_inter = 0;

// The data is allocated given the new size
// The data is allocated given the new size and filled in with zero to represent an identity displacement field
deformationFieldImage.realloc();
// The image is filled in with zero to represent an identity displacement field
reg_tools_multiplyValueToImage(deformationFieldImage, deformationFieldImage, 0.f);
deformationFieldImage->intent_p1 = DISP_FIELD;
// The displacement field is converted into a deformation field
reg_getDeformationFromDisplacement(deformationFieldImage);
Expand Down Expand Up @@ -1699,7 +1690,6 @@ void reg_voxelCentric2NodeCentric(nifti_image *nodeImage,
} // loop over z
}
/* *************************************************************** */
extern "C++"
void reg_voxelCentric2NodeCentric(nifti_image * nodeImage,
nifti_image * voxelImage,
float weight,
Expand Down Expand Up @@ -2148,7 +2138,6 @@ void reg_spline_refineControlPointGrid3D(nifti_image *splineControlPoint, nifti_
free(oldGrid);
}
/* *************************************************************** */
extern "C++"
void reg_spline_refineControlPointGrid(nifti_image *controlPointGrid,
nifti_image *referenceImage) {
NR_DEBUG("Starting the refine the control point grid");
Expand Down
Loading

0 comments on commit 327d516

Please sign in to comment.