Skip to content

Commit

Permalink
Use NiftiImage in Content classes
Browse files Browse the repository at this point in the history
  • Loading branch information
@onurulgen
onurulgen committed Aug 27, 2024
1 parent 2915900 commit d515493
Show file tree
Hide file tree
Showing 76 changed files with 579 additions and 824 deletions.
2 changes: 1 addition & 1 deletion niftyreg_build_version.txt
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
432
433
2 changes: 1 addition & 1 deletion reg-apps/reg_tools.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1071,7 +1071,7 @@ int main(int argc, char **argv)
if(image->datatype!=NIFTI_TYPE_FLOAT32)
reg_tools_changeDatatype<float>(image);
// Create a temporary mask
const size_t voxelNumber = NiftiImage::calcVoxelNumber(image, 3);
const size_t voxelNumber = image.nVoxelsPerVolume();
int *temp_mask = (int *)malloc(voxelNumber * sizeof(int));
for (size_t i = 0; i < voxelNumber; ++i)
temp_mask[i]=i;
Expand Down
6 changes: 2 additions & 4 deletions reg-apps/reg_transform.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1061,20 +1061,18 @@ int main(int argc, char **argv) {
tempField->ndim = tempField->dim[0] = 5;
tempField->nt = tempField->dim[4] = 1;
tempField->nu = tempField->dim[5] = tempField->nz > 1 ? 3 : 2;
tempField->nvox = NiftiImage::calcVoxelNumber(tempField, tempField->ndim);
tempField->nbyper = inputTransImage->nbyper;
tempField->datatype = inputTransImage->datatype;
tempField.realloc();
tempField->intent_code = NIFTI_INTENT_VECTOR;
memset(tempField->intent_name, 0, 16);
strcpy(tempField->intent_name, "NREG_TRANS");
tempField.setIntentName("NREG_TRANS"s);
tempField->intent_p1 = DEF_FIELD;
if (inputTransImage->intent_p1 == SPLINE_VEL_GRID) {
tempField->intent_p1 = DEF_VEL_FIELD;
tempField->intent_p2 = inputTransImage->intent_p2;
}
tempField->scl_slope = 1.f;
tempField->scl_inter = 0.f;
tempField->data = calloc(tempField->nvox, tempField->nbyper);
// Compute the dense field
if (inputTransImage->intent_p1 == LIN_SPLINE_GRID ||
inputTransImage->intent_p1 == CUB_SPLINE_GRID)
Expand Down
33 changes: 20 additions & 13 deletions reg-io/RNifti/NiftiImage.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1550,17 +1550,6 @@ class NiftiImage
**/
virtual ~NiftiImage () { release(); }

/**
* Disown the wrapped pointer, removing responsibility for freeing it upon destruction
* @return The wrapped pointer
*/
nifti_image* disown ()
{
nifti_image *img = image;
image = nullptr;
return img;
}

/**
* Allows a \c NiftiImage object to be treated as a pointer to a \c const \c nifti_image
**/
Expand Down Expand Up @@ -1813,6 +1802,7 @@ class NiftiImage

/**
* Return the datatype of the image
* @param image A pointer to a NIfTI image
* @return A variant holding a NIfTI datatype
**/
static DataType getDataType (const nifti_image *image)
Expand All @@ -1837,6 +1827,9 @@ class NiftiImage
}
}

// Delete the overload that accepts NiftiImage; use the member function instead
static DataType getDataType(const NiftiImage&) = delete;

/**
* Return the datatype of the image
* @return A variant holding a NIfTI datatype
Expand All @@ -1845,6 +1838,7 @@ class NiftiImage

/**
* Return the datatype of the image, if it is a floating-point type
* @param image A pointer to a NIfTI image
* @return A variant holding a NIfTI datatype
**/
static std::variant<float, double> getFloatingDataType (const nifti_image *image)
Expand All @@ -1861,6 +1855,9 @@ class NiftiImage
}
}

// Delete the overload that accepts NiftiImage; use the member function instead
static std::variant<float, double> getFloatingDataType(const NiftiImage&) = delete;

/**
* Return the datatype of the image, if it is a floating-point type
* @return A variant holding a NIfTI datatype
Expand Down Expand Up @@ -1888,11 +1885,11 @@ class NiftiImage

/**
* Copy the pixel data from another image
* @param other The image from which to copy the data
* @param source The image from which to copy the data
* @exception runtime_error If the lengths and datatypes of the two images do not match
* @return Self, after copying the data
**/
NiftiImage & copyData (const nifti_image *other);
NiftiImage & copyData (const nifti_image *source);

/**
* Drop the data from the image, retaining only the metadata. This method invalidates any
Expand Down Expand Up @@ -2116,6 +2113,16 @@ class NiftiImage
return voxelNumber;
}

// Delete the overload that accepts NiftiImage; use the member function instead
static size_t calcVoxelNumber(const NiftiImage&, const int) = delete;

/**
* Calculate the number of voxels in the image
* @param dimCount Number of dimensions to consider
* @return The number of voxels in the image
*/
size_t calcVoxelNumber(const int dimCount) const { return calcVoxelNumber(image, dimCount); }

/**
* Recalculate the number of voxels in the image and update the nvox field
*/
Expand Down
18 changes: 9 additions & 9 deletions reg-io/RNifti/NiftiImage_impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1871,23 +1871,23 @@ inline NiftiImage & NiftiImage::replaceData (const NiftiImageData &data)
return *this;
}

inline NiftiImage & NiftiImage::copyData (const nifti_image *other)
inline NiftiImage & NiftiImage::copyData (const nifti_image *source)
{
if (this->isNull())
return *this;
else if (other == nullptr || other->data == nullptr)
else if (source == nullptr || source->data == nullptr)
throw std::runtime_error("Cannot copy data from a null image");
else if (other->nvox != image->nvox)
else if (source->nvox != image->nvox)
throw std::runtime_error("Cannot copy data from an image with a different length");
else if (other->datatype != image->datatype)
else if (source->datatype != image->datatype || source->nbyper != image->nbyper)
throw std::runtime_error("Cannot copy data from an image with a different datatype");

// Copy the data
memcpy(image->data, other->data, totalBytes());
image->scl_slope = other->scl_slope;
image->scl_inter = other->scl_inter;
image->cal_min = other->cal_min;
image->cal_max = other->cal_max;
memcpy(image->data, source->data, totalBytes());
image->scl_slope = source->scl_slope;
image->scl_inter = source->scl_inter;
image->cal_min = source->cal_min;
image->cal_max = source->cal_max;

return *this;
}
Expand Down
4 changes: 2 additions & 2 deletions reg-lib/AladinContent.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@
using namespace std;

/* *************************************************************** */
AladinContent::AladinContent(nifti_image *referenceIn,
nifti_image *floatingIn,
AladinContent::AladinContent(NiftiImage& referenceIn,
NiftiImage& floatingIn,
int *referenceMaskIn,
mat44 *transformationMatrixIn,
size_t bytesIn,
Expand Down
6 changes: 3 additions & 3 deletions reg-lib/AladinContent.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@
class AladinContent: public Content {
public:
AladinContent(const AladinContent&) = delete;
AladinContent(nifti_image *referenceIn,
nifti_image *floatingIn,
AladinContent(NiftiImage& referenceIn,
NiftiImage& floatingIn,
int *referenceMaskIn = nullptr,
mat44 *transformationMatrixIn = nullptr,
size_t bytesIn = sizeof(float),
Expand All @@ -28,7 +28,7 @@ class AladinContent: public Content {
virtual _reg_blockMatchingParam* GetBlockMatchingParams() { return blockMatchingParams; }

protected:
_reg_blockMatchingParam* blockMatchingParams;
_reg_blockMatchingParam *blockMatchingParams;
unsigned currentPercentageOfBlockToUse;
unsigned inlierLts;
int stepSizeBlock;
Expand Down
4 changes: 2 additions & 2 deletions reg-lib/AladinContentCreator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@

class AladinContentCreator: public ContentCreator {
public:
virtual AladinContent* Create(nifti_image *reference,
nifti_image *floating,
virtual AladinContent* Create(NiftiImage& reference,
NiftiImage& floating,
int *referenceMask = nullptr,
mat44 *transformationMatrix = nullptr,
size_t bytes = sizeof(float),
Expand Down
81 changes: 33 additions & 48 deletions reg-lib/Compute.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@ void Compute::UpdateControlPointPosition(float *currentDof,
const bool optimiseX,
const bool optimiseY,
const bool optimiseZ) {
const nifti_image *controlPointGrid = dynamic_cast<F3dContent&>(con).F3dContent::GetControlPointGrid();
const NiftiImage& controlPointGrid = dynamic_cast<F3dContent&>(con).F3dContent::GetControlPointGrid();
if (optimiseX && optimiseY && optimiseZ) {
// Update the values for all axis displacement
for (size_t i = 0; i < controlPointGrid->nvox; ++i)
Expand Down Expand Up @@ -139,7 +139,7 @@ void Compute::GetImageGradient(int interpolation, float paddingValue, int active
/* *************************************************************** */
double Compute::GetMaximalLength(bool optimiseX, bool optimiseY, bool optimiseZ) {
if (!optimiseX && !optimiseY && !optimiseZ) return 0;
const nifti_image *transformationGradient = dynamic_cast<F3dContent&>(con).GetTransformationGradient();
const NiftiImage& transformationGradient = dynamic_cast<F3dContent&>(con).GetTransformationGradient();
switch (transformationGradient->datatype) {
case NIFTI_TYPE_FLOAT32:
return reg_getMaximalLength<float>(transformationGradient, optimiseX, optimiseY, optimiseZ);
Expand All @@ -151,7 +151,7 @@ double Compute::GetMaximalLength(bool optimiseX, bool optimiseY, bool optimiseZ)
/* *************************************************************** */
void Compute::NormaliseGradient(double maxGradLength, bool optimiseX, bool optimiseY, bool optimiseZ) {
if (maxGradLength == 0 || (!optimiseX && !optimiseY && !optimiseZ)) return;
NiftiImage transformationGradient = dynamic_cast<F3dContent&>(con).GetTransformationGradient();
NiftiImage& transformationGradient = dynamic_cast<F3dContent&>(con).GetTransformationGradient();
const bool hasZ = transformationGradient->nz > 1;
if (!hasZ) optimiseZ = false;
NiftiImageData ptrX = transformationGradient.data(0);
Expand Down Expand Up @@ -191,8 +191,8 @@ void Compute::SmoothGradient(float sigma) {
/* *************************************************************** */
void Compute::GetApproximatedGradient(InterfaceOptimiser& opt) {
F3dContent& con = dynamic_cast<F3dContent&>(this->con);
nifti_image *controlPointGrid = con.GetControlPointGrid();
nifti_image *transformationGradient = con.GetTransformationGradient();
NiftiImage& controlPointGrid = con.GetControlPointGrid();
NiftiImage& transformationGradient = con.GetTransformationGradient();
std::visit([&](auto&& cppDataType) {
using Type = std::decay_t<decltype(cppDataType)>;

Expand All @@ -217,7 +217,7 @@ void Compute::GetApproximatedGradient(InterfaceOptimiser& opt) {
// Update the changes for GPU
con.UpdateControlPointGrid();
con.UpdateTransformationGradient();
}, NiftiImage::getFloatingDataType(controlPointGrid));
}, controlPointGrid.getFloatingDataType());
}
/* *************************************************************** */
void Compute::GetDefFieldFromVelocityGrid(const bool updateStepNumber) {
Expand All @@ -227,8 +227,8 @@ void Compute::GetDefFieldFromVelocityGrid(const bool updateStepNumber) {
updateStepNumber);
}
/* *************************************************************** */
void Compute::ConvolveImage(nifti_image *image) {
const nifti_image *controlPointGrid = dynamic_cast<F3dContent&>(con).F3dContent::GetControlPointGrid();
void Compute::ConvolveImage(NiftiImage& image) {
const NiftiImage& controlPointGrid = dynamic_cast<F3dContent&>(con).F3dContent::GetControlPointGrid();
constexpr ConvKernelType kernelType = ConvKernelType::Cubic;
float currentNodeSpacing[3];
currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = controlPointGrid->dx;
Expand Down Expand Up @@ -283,27 +283,27 @@ void Compute::ConvolveVoxelBasedMeasureGradient(float weight) {
void Compute::ExponentiateGradient(Content& conBwIn) {
F3dContent& con = dynamic_cast<F3dContent&>(this->con);
F3dContent& conBw = dynamic_cast<F3dContent&>(conBwIn);
const nifti_image *deformationField = con.Content::GetDeformationField();
nifti_image *voxelBasedMeasureGradient = con.GetVoxelBasedMeasureGradient();
nifti_image *controlPointGridBw = conBw.GetControlPointGrid();
const NiftiImage& deformationField = con.Content::GetDeformationField();
NiftiImage& voxelBasedMeasureGradient = con.GetVoxelBasedMeasureGradient();
NiftiImage& controlPointGridBw = conBw.GetControlPointGrid();
mat44 *affineTransformationBw = conBw.GetTransformationMatrix();
const size_t compNum = size_t(fabs(controlPointGridBw->intent_p2)); // The number of composition

/* Allocate a temporary gradient image to store the backward gradient */
nifti_image *tempGrad = nifti_dup(*voxelBasedMeasureGradient, false);
NiftiImage tempGrad(voxelBasedMeasureGradient, NiftiImage::Copy::ImageInfoAndAllocData);

// Create all deformation field images needed for resampling
nifti_image **tempDef = (nifti_image**)malloc((compNum + 1) * sizeof(nifti_image*));
unique_ptr<NiftiImage[]> tempDef(new NiftiImage[compNum + 1]);
for (size_t i = 0; i <= compNum; ++i)
tempDef[i] = nifti_dup(*deformationField, false);
tempDef[i] = NiftiImage(deformationField, NiftiImage::Copy::ImageInfoAndAllocData);

// Generate all intermediate deformation fields
reg_spline_getIntermediateDefFieldFromVelGrid(controlPointGridBw, tempDef);
reg_spline_getIntermediateDefFieldFromVelGrid(controlPointGridBw, tempDef.get());

// Remove the affine component
nifti_image *affineDisp = nullptr;
NiftiImage affineDisp;
if (affineTransformationBw) {
affineDisp = nifti_dup(*deformationField, false);
affineDisp = NiftiImage(deformationField, NiftiImage::Copy::ImageInfoAndAllocData);
reg_affine_getDeformationField(affineTransformationBw, affineDisp);
reg_getDisplacementFromDeformation(affineDisp);
}
Expand All @@ -325,25 +325,18 @@ void Compute::ExponentiateGradient(Content& conBwIn) {
reg_tools_divideValueToImage(voxelBasedMeasureGradient, // in
voxelBasedMeasureGradient, // out
pow(2, compNum)); // value

for (size_t i = 0; i <= compNum; ++i)
nifti_image_free(tempDef[i]);
free(tempDef);
nifti_image_free(tempGrad);
if (affineDisp)
nifti_image_free(affineDisp);
}
/* *************************************************************** */
nifti_image* Compute::ScaleGradient(const nifti_image& transformationGradient, float scale) {
nifti_image *scaledGradient = nifti_dup(transformationGradient, false);
reg_tools_multiplyValueToImage(&transformationGradient, scaledGradient, scale);
NiftiImage Compute::ScaleGradient(const NiftiImage& transformationGradient, float scale) {
NiftiImage scaledGradient(transformationGradient, NiftiImage::Copy::ImageInfoAndAllocData);
reg_tools_multiplyValueToImage(transformationGradient, scaledGradient, scale);
return scaledGradient;
}
/* *************************************************************** */
void Compute::UpdateVelocityField(float scale, bool optimiseX, bool optimiseY, bool optimiseZ) {
F3dContent& con = dynamic_cast<F3dContent&>(this->con);
nifti_image *scaledGradient = ScaleGradient(*con.GetTransformationGradient(), scale);
nifti_image *controlPointGrid = con.GetControlPointGrid();
NiftiImage scaledGradient = ScaleGradient(con.GetTransformationGradient(), scale);
NiftiImage& controlPointGrid = con.GetControlPointGrid();

// Reset the gradient along the axes if appropriate
reg_setGradientToZero(scaledGradient, !optimiseX, !optimiseY, !optimiseZ);
Expand All @@ -352,36 +345,31 @@ void Compute::UpdateVelocityField(float scale, bool optimiseX, bool optimiseY, b
reg_tools_addImageToImage(controlPointGrid, // in
scaledGradient, // in
controlPointGrid); // out

nifti_image_free(scaledGradient);
}
/* *************************************************************** */
void Compute::BchUpdate(float scale, int bchUpdateValue) {
F3dContent& con = dynamic_cast<F3dContent&>(this->con);
nifti_image *scaledGradient = ScaleGradient(*con.GetTransformationGradient(), scale);
nifti_image *controlPointGrid = con.GetControlPointGrid();

NiftiImage scaledGradient = ScaleGradient(con.GetTransformationGradient(), scale);
NiftiImage& controlPointGrid = con.GetControlPointGrid();
compute_BCH_update(controlPointGrid, scaledGradient, bchUpdateValue);

nifti_image_free(scaledGradient);
}
/* *************************************************************** */
void Compute::SymmetriseVelocityFields(Content& conBwIn) {
nifti_image *controlPointGrid = dynamic_cast<F3dContent&>(this->con).GetControlPointGrid();
nifti_image *controlPointGridBw = dynamic_cast<F3dContent&>(conBwIn).GetControlPointGrid();
NiftiImage& controlPointGrid = dynamic_cast<F3dContent&>(this->con).GetControlPointGrid();
NiftiImage& controlPointGridBw = dynamic_cast<F3dContent&>(conBwIn).GetControlPointGrid();

// In order to ensure symmetry, the forward and backward velocity fields
// are averaged in both image spaces: reference and floating
nifti_image *warpedTrans = nifti_dup(*controlPointGridBw, false);
nifti_image *warpedTransBw = nifti_dup(*controlPointGrid, false);
NiftiImage warpedTrans(controlPointGridBw, NiftiImage::Copy::ImageInfoAndAllocData);
NiftiImage warpedTransBw(controlPointGrid, NiftiImage::Copy::ImageInfoAndAllocData);

// Both parametrisations are converted into displacement
reg_getDisplacementFromDeformation(controlPointGrid);
reg_getDisplacementFromDeformation(controlPointGridBw);

// Both parametrisations are copied over
memcpy(warpedTransBw->data, controlPointGridBw->data, warpedTransBw->nvox * warpedTransBw->nbyper);
memcpy(warpedTrans->data, controlPointGrid->data, warpedTrans->nvox * warpedTrans->nbyper);
warpedTrans.copyData(controlPointGrid);
warpedTransBw.copyData(controlPointGridBw);

// and subtracted (sum and negation)
reg_tools_subtractImageFromImage(controlPointGridBw, // displacement
Expand All @@ -402,19 +390,16 @@ void Compute::SymmetriseVelocityFields(Content& conBwIn) {
// Convert the velocity field from displacement to deformation
reg_getDeformationFromDisplacement(controlPointGrid);
reg_getDeformationFromDisplacement(controlPointGridBw);

nifti_image_free(warpedTrans);
nifti_image_free(warpedTransBw);
}
/* *************************************************************** */
void Compute::DefFieldCompose(const nifti_image *defField) {
void Compute::DefFieldCompose(const NiftiImage& defField) {
reg_defField_compose(defField, con.GetDeformationField(), nullptr);
}
/* *************************************************************** */
NiftiImage Compute::ResampleGradient(int interpolation, float padding) {
DefContent& con = dynamic_cast<DefContent&>(this->con);
nifti_image *voxelBasedMeasureGradient = con.GetVoxelBasedMeasureGradient();
NiftiImage warpedImage = NiftiImage(voxelBasedMeasureGradient, NiftiImage::Copy::ImageInfoAndAllocData);
NiftiImage& voxelBasedMeasureGradient = con.GetVoxelBasedMeasureGradient();
NiftiImage warpedImage(voxelBasedMeasureGradient, NiftiImage::Copy::ImageInfoAndAllocData);
reg_resampleGradient(voxelBasedMeasureGradient, warpedImage, con.GetDeformationField(), interpolation, padding);
return warpedImage;
}
Expand Down
Loading

0 comments on commit d515493

Please sign in to comment.