diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 84796bf1..21fbd2eb 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -432 +433 diff --git a/reg-apps/reg_tools.cpp b/reg-apps/reg_tools.cpp index 451def0d..033077e5 100755 --- a/reg-apps/reg_tools.cpp +++ b/reg-apps/reg_tools.cpp @@ -1071,7 +1071,7 @@ int main(int argc, char **argv) if(image->datatype!=NIFTI_TYPE_FLOAT32) reg_tools_changeDatatype(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; diff --git a/reg-apps/reg_transform.cpp b/reg-apps/reg_transform.cpp index 4ab60681..cb19fd66 100755 --- a/reg-apps/reg_transform.cpp +++ b/reg-apps/reg_transform.cpp @@ -1061,12 +1061,11 @@ 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; @@ -1074,7 +1073,6 @@ int main(int argc, char **argv) { } 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) diff --git a/reg-io/RNifti/NiftiImage.h b/reg-io/RNifti/NiftiImage.h index 16461ddd..20646c77 100644 --- a/reg-io/RNifti/NiftiImage.h +++ b/reg-io/RNifti/NiftiImage.h @@ -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 **/ @@ -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) @@ -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 @@ -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 getFloatingDataType (const nifti_image *image) @@ -1861,6 +1855,9 @@ class NiftiImage } } + // Delete the overload that accepts NiftiImage; use the member function instead + static std::variant getFloatingDataType(const NiftiImage&) = delete; + /** * Return the datatype of the image, if it is a floating-point type * @return A variant holding a NIfTI datatype @@ -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 @@ -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 */ diff --git a/reg-io/RNifti/NiftiImage_impl.h b/reg-io/RNifti/NiftiImage_impl.h index 7672f407..15e503ec 100644 --- a/reg-io/RNifti/NiftiImage_impl.h +++ b/reg-io/RNifti/NiftiImage_impl.h @@ -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; } diff --git a/reg-lib/AladinContent.cpp b/reg-lib/AladinContent.cpp index ab1a07af..f8299cdb 100755 --- a/reg-lib/AladinContent.cpp +++ b/reg-lib/AladinContent.cpp @@ -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, diff --git a/reg-lib/AladinContent.h b/reg-lib/AladinContent.h index 19cf8c28..8c8ba85b 100755 --- a/reg-lib/AladinContent.h +++ b/reg-lib/AladinContent.h @@ -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), @@ -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; diff --git a/reg-lib/AladinContentCreator.h b/reg-lib/AladinContentCreator.h index 91d03be8..939fa524 100644 --- a/reg-lib/AladinContentCreator.h +++ b/reg-lib/AladinContentCreator.h @@ -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), diff --git a/reg-lib/Compute.cpp b/reg-lib/Compute.cpp index 329e48fd..5b23e2da 100644 --- a/reg-lib/Compute.cpp +++ b/reg-lib/Compute.cpp @@ -95,7 +95,7 @@ void Compute::UpdateControlPointPosition(float *currentDof, const bool optimiseX, const bool optimiseY, const bool optimiseZ) { - const nifti_image *controlPointGrid = dynamic_cast(con).F3dContent::GetControlPointGrid(); + const NiftiImage& controlPointGrid = dynamic_cast(con).F3dContent::GetControlPointGrid(); if (optimiseX && optimiseY && optimiseZ) { // Update the values for all axis displacement for (size_t i = 0; i < controlPointGrid->nvox; ++i) @@ -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(con).GetTransformationGradient(); + const NiftiImage& transformationGradient = dynamic_cast(con).GetTransformationGradient(); switch (transformationGradient->datatype) { case NIFTI_TYPE_FLOAT32: return reg_getMaximalLength(transformationGradient, optimiseX, optimiseY, optimiseZ); @@ -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(con).GetTransformationGradient(); + NiftiImage& transformationGradient = dynamic_cast(con).GetTransformationGradient(); const bool hasZ = transformationGradient->nz > 1; if (!hasZ) optimiseZ = false; NiftiImageData ptrX = transformationGradient.data(0); @@ -191,8 +191,8 @@ void Compute::SmoothGradient(float sigma) { /* *************************************************************** */ void Compute::GetApproximatedGradient(InterfaceOptimiser& opt) { F3dContent& con = dynamic_cast(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; @@ -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) { @@ -227,8 +227,8 @@ void Compute::GetDefFieldFromVelocityGrid(const bool updateStepNumber) { updateStepNumber); } /* *************************************************************** */ -void Compute::ConvolveImage(nifti_image *image) { - const nifti_image *controlPointGrid = dynamic_cast(con).F3dContent::GetControlPointGrid(); +void Compute::ConvolveImage(NiftiImage& image) { + const NiftiImage& controlPointGrid = dynamic_cast(con).F3dContent::GetControlPointGrid(); constexpr ConvKernelType kernelType = ConvKernelType::Cubic; float currentNodeSpacing[3]; currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = controlPointGrid->dx; @@ -283,27 +283,27 @@ void Compute::ConvolveVoxelBasedMeasureGradient(float weight) { void Compute::ExponentiateGradient(Content& conBwIn) { F3dContent& con = dynamic_cast(this->con); F3dContent& conBw = dynamic_cast(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 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); } @@ -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(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); @@ -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(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(this->con).GetControlPointGrid(); - nifti_image *controlPointGridBw = dynamic_cast(conBwIn).GetControlPointGrid(); + NiftiImage& controlPointGrid = dynamic_cast(this->con).GetControlPointGrid(); + NiftiImage& controlPointGridBw = dynamic_cast(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 @@ -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(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; } diff --git a/reg-lib/Compute.h b/reg-lib/Compute.h index fdf3e673..a8a8852d 100644 --- a/reg-lib/Compute.h +++ b/reg-lib/Compute.h @@ -39,11 +39,11 @@ class Compute { #ifdef NR_TESTING public: #endif - virtual void DefFieldCompose(const nifti_image *defField); + virtual void DefFieldCompose(const NiftiImage& defField); virtual NiftiImage ResampleGradient(int interpolation, float padding); virtual void VoxelCentricToNodeCentric(float weight); private: - void ConvolveImage(nifti_image*); - nifti_image* ScaleGradient(const nifti_image&, float); + void ConvolveImage(NiftiImage&); + NiftiImage ScaleGradient(const NiftiImage&, float); }; diff --git a/reg-lib/Content.cpp b/reg-lib/Content.cpp index b64a48b8..2b8fbcbd 100644 --- a/reg-lib/Content.cpp +++ b/reg-lib/Content.cpp @@ -2,63 +2,49 @@ #include "_reg_tools.h" /* *************************************************************** */ -Content::Content(nifti_image *referenceIn, - nifti_image *floatingIn, +Content::Content(NiftiImage& referenceIn, + NiftiImage& floatingIn, int *referenceMaskIn, mat44 *transformationMatrixIn, size_t bytesIn): - reference(referenceIn), - floating(floatingIn), + reference(NiftiImage(referenceIn, NiftiImage::Copy::Acquire)), + floating(NiftiImage(floatingIn, NiftiImage::Copy::Acquire)), referenceMask(referenceMaskIn), transformationMatrix(transformationMatrixIn) { if (!referenceIn || !floatingIn) NR_FATAL_ERROR("referenceIn or floatingIn can't be nullptr"); AllocateWarped(); AllocateDeformationField(bytesIn); - activeVoxelNumber = NiftiImage::calcVoxelNumber(reference, 3); + activeVoxelNumber = reference.nVoxelsPerVolume(); if (!referenceMask) { referenceMaskManaged.reset(new int[activeVoxelNumber]()); referenceMask = referenceMaskManaged.get(); } } /* *************************************************************** */ -Content::~Content() { - DeallocateWarped(); - DeallocateDeformationField(); -} -/* *************************************************************** */ void Content::AllocateWarped() { - warped = nifti_copy_nim_info(reference); - warped->dim[0] = warped->ndim = floating->ndim; - warped->dim[4] = warped->nt = floating->nt; - warped->pixdim[4] = warped->dt = 1; - warped->nvox = NiftiImage::calcVoxelNumber(warped, warped->ndim); + warped = NiftiImage(reference, NiftiImage::Copy::ImageInfo); + warped.setDim(NiftiDim::NDim, floating->ndim); + warped.setDim(NiftiDim::T, floating->nt); + warped.setPixDim(NiftiDim::T, 1); warped->datatype = floating->datatype; warped->nbyper = floating->nbyper; - warped->data = calloc(warped->nvox, warped->nbyper); -} -/* *************************************************************** */ -void Content::DeallocateWarped() { - if (warped) { - nifti_image_free(warped); - warped = nullptr; - } + warped.realloc(); } /* *************************************************************** */ void Content::AllocateDeformationField(size_t bytes) { - deformationField = nifti_copy_nim_info(reference); - deformationField->dim[0] = deformationField->ndim = 5; - if (reference->dim[0] == 2) - deformationField->dim[3] = deformationField->nz = 1; - deformationField->dim[4] = deformationField->nt = 1; - deformationField->pixdim[4] = deformationField->dt = 1; - deformationField->dim[5] = deformationField->nu = reference->nz > 1 ? 3 : 2; - deformationField->pixdim[5] = deformationField->du = 1; - deformationField->dim[6] = deformationField->nv = 1; - deformationField->pixdim[6] = deformationField->dv = 1; - deformationField->dim[7] = deformationField->nw = 1; - deformationField->pixdim[7] = deformationField->dw = 1; - deformationField->nvox = NiftiImage::calcVoxelNumber(deformationField, deformationField->ndim); + deformationField = NiftiImage(reference, NiftiImage::Copy::ImageInfo); + deformationField.setDim(NiftiDim::NDim, 5); + if (reference->ndim == 2) + deformationField.setDim(NiftiDim::Z, 1); + deformationField.setDim(NiftiDim::T, 1); + deformationField.setPixDim(NiftiDim::T, 1); + deformationField.setDim(NiftiDim::U, reference->nz > 1 ? 3 : 2); + deformationField.setPixDim(NiftiDim::U, 1); + deformationField.setDim(NiftiDim::V, 1); + deformationField.setPixDim(NiftiDim::V, 1); + deformationField.setDim(NiftiDim::W, 1); + deformationField.setPixDim(NiftiDim::W, 1); deformationField->nbyper = (int)bytes; if (bytes == 4) deformationField->datatype = NIFTI_TYPE_FLOAT32; @@ -67,21 +53,13 @@ void Content::AllocateDeformationField(size_t bytes) { else NR_FATAL_ERROR("Only float or double are expected for the deformation field"); deformationField->intent_code = NIFTI_INTENT_VECTOR; - memset(deformationField->intent_name, 0, sizeof(deformationField->intent_name)); - strcpy(deformationField->intent_name, "NREG_TRANS"); + deformationField.setIntentName("NREG_TRANS"s); // First create a displacement field filled with 0 to obtain an identity disp deformationField->intent_p1 = DISP_FIELD; deformationField->scl_slope = 1; deformationField->scl_inter = 0; - deformationField->data = calloc(deformationField->nvox, deformationField->nbyper); + deformationField.realloc(); // Convert to an identity deformation field reg_getDeformationFromDisplacement(deformationField); } /* *************************************************************** */ -void Content::DeallocateDeformationField() { - if (deformationField) { - nifti_image_free(deformationField); - deformationField = nullptr; - } -} -/* *************************************************************** */ diff --git a/reg-lib/Content.h b/reg-lib/Content.h index 7beb9e4a..7dd16957 100644 --- a/reg-lib/Content.h +++ b/reg-lib/Content.h @@ -5,23 +5,23 @@ class Content { public: Content() = delete; // Can't be initialised without reference and floating images - Content(nifti_image *referenceIn, - nifti_image *floatingIn, + Content(NiftiImage& referenceIn, + NiftiImage& floatingIn, int *referenceMaskIn = nullptr, mat44 *transformationMatrixIn = nullptr, size_t bytesIn = sizeof(float)); - virtual ~Content(); + virtual ~Content() = default; virtual bool IsCurrentComputationDoubleCapable() { return true; } // Getters virtual size_t GetActiveVoxelNumber() { return activeVoxelNumber; } - virtual nifti_image* GetReference() { return reference; } - virtual nifti_image* GetFloating() { return floating; } - virtual nifti_image* GetDeformationField() { return deformationField; } + virtual NiftiImage& GetReference() { return reference; } + virtual NiftiImage& GetFloating() { return floating; } + virtual NiftiImage& GetDeformationField() { return deformationField; } virtual int* GetReferenceMask() { return referenceMask; } virtual mat44* GetTransformationMatrix() { return transformationMatrix; } - virtual nifti_image* GetWarped() { return warped; } + virtual NiftiImage& GetWarped() { return warped; } // Methods for transferring data from nifti to device virtual void UpdateDeformationField() {} @@ -37,19 +37,17 @@ class Content { protected: size_t activeVoxelNumber = 0; - nifti_image *reference = nullptr; - nifti_image *floating = nullptr; - nifti_image *deformationField = nullptr; + NiftiImage reference; + NiftiImage floating; + NiftiImage deformationField; int *referenceMask = nullptr; unique_ptr referenceMaskManaged; mat44 *transformationMatrix = nullptr; - nifti_image *warped = nullptr; + NiftiImage warped; private: void AllocateWarped(); - void DeallocateWarped(); void AllocateDeformationField(size_t bytes); - void DeallocateDeformationField(); #ifdef NR_TESTING public: @@ -57,8 +55,8 @@ class Content { protected: #endif // Functions for testing - virtual void SetDeformationField(nifti_image *deformationFieldIn) { DeallocateDeformationField(); deformationField = deformationFieldIn; } + virtual void SetDeformationField(NiftiImage&& deformationFieldIn) { deformationField = std::move(deformationFieldIn); } virtual void SetReferenceMask(int *referenceMaskIn) { referenceMask = referenceMaskIn; } virtual void SetTransformationMatrix(mat44 *transformationMatrixIn) { transformationMatrix = transformationMatrixIn; } - virtual void SetWarped(nifti_image *warpedIn) { DeallocateWarped(); warped = warpedIn; } + virtual void SetWarped(NiftiImage&& warpedIn) { warped = std::move(warpedIn); } }; diff --git a/reg-lib/ContentCreator.h b/reg-lib/ContentCreator.h index 050bdba8..52b586f8 100644 --- a/reg-lib/ContentCreator.h +++ b/reg-lib/ContentCreator.h @@ -4,8 +4,8 @@ class ContentCreator { public: - virtual Content* Create(nifti_image *reference, - nifti_image *floating, + virtual Content* Create(NiftiImage& reference, + NiftiImage& floating, int *referenceMask = nullptr, mat44 *transformationMatrix = nullptr, size_t bytes = sizeof(float)) { diff --git a/reg-lib/DefContent.cpp b/reg-lib/DefContent.cpp index 6885153e..b8fc2a57 100644 --- a/reg-lib/DefContent.cpp +++ b/reg-lib/DefContent.cpp @@ -2,66 +2,30 @@ #include "_reg_resampling.h" /* *************************************************************** */ -DefContent::DefContent(nifti_image *referenceIn, - nifti_image *floatingIn, - nifti_image *localWeightSimIn, +DefContent::DefContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, + NiftiImage *localWeightSimIn, int *referenceMaskIn, mat44 *transformationMatrixIn, size_t bytesIn): Content(referenceIn, floatingIn, referenceMaskIn, transformationMatrixIn, bytesIn) { - AllocateWarpedGradient(); - AllocateVoxelBasedMeasureGradient(); - AllocateLocalWeightSim(localWeightSimIn); -} -/* *************************************************************** */ -DefContent::~DefContent() { - DeallocateWarpedGradient(); - DeallocateVoxelBasedMeasureGradient(); - DeallocateLocalWeightSim(); -} -/* *************************************************************** */ -void DefContent::AllocateLocalWeightSim(nifti_image *localWeightSimIn) { - if (!localWeightSimIn) return; - localWeightSim = nifti_copy_nim_info(reference); - localWeightSim->dim[0] = localWeightSim->ndim = localWeightSimIn->dim[0]; - localWeightSim->dim[4] = localWeightSim->nt = localWeightSimIn->dim[4]; - localWeightSim->dim[5] = localWeightSim->nu = localWeightSimIn->dim[5]; - localWeightSim->nvox = NiftiImage::calcVoxelNumber(localWeightSim, localWeightSim->ndim); - localWeightSim->data = malloc(localWeightSim->nvox * localWeightSim->nbyper); + warpedGradient = NiftiImage(deformationField, NiftiImage::Copy::ImageInfoAndAllocData); + voxelBasedMeasureGradient = NiftiImage(deformationField, NiftiImage::Copy::ImageInfoAndAllocData); + if (localWeightSimIn && *localWeightSimIn) + AllocateLocalWeightSim(*localWeightSimIn); +} +/* *************************************************************** */ +void DefContent::AllocateLocalWeightSim(NiftiImage& localWeightSimIn) { + localWeightSim = NiftiImage(reference, NiftiImage::Copy::ImageInfo); + localWeightSim.setDim(NiftiDim::NDim, localWeightSimIn->dim[0]); + localWeightSim.setDim(NiftiDim::T, localWeightSimIn->dim[4]); + localWeightSim.setDim(NiftiDim::U, localWeightSimIn->dim[5]); + localWeightSim.realloc(); reg_getDeformationFromDisplacement(voxelBasedMeasureGradient); reg_resampleImage(localWeightSimIn, localWeightSim, voxelBasedMeasureGradient, nullptr, 1, 0); } /* *************************************************************** */ -void DefContent::DeallocateLocalWeightSim() { - if (localWeightSim) { - nifti_image_free(localWeightSim); - localWeightSim = nullptr; - } -} -/* *************************************************************** */ -void DefContent::AllocateWarpedGradient() { - warpedGradient = nifti_dup(*deformationField, false); -} -/* *************************************************************** */ -void DefContent::DeallocateWarpedGradient() { - if (warpedGradient) { - nifti_image_free(warpedGradient); - warpedGradient = nullptr; - } -} -/* *************************************************************** */ -void DefContent::AllocateVoxelBasedMeasureGradient() { - voxelBasedMeasureGradient = nifti_dup(*deformationField, false); -} -/* *************************************************************** */ -void DefContent::DeallocateVoxelBasedMeasureGradient() { - if (voxelBasedMeasureGradient) { - nifti_image_free(voxelBasedMeasureGradient); - voxelBasedMeasureGradient = nullptr; - } -} -/* *************************************************************** */ void DefContent::ZeroVoxelBasedMeasureGradient() { - memset(voxelBasedMeasureGradient->data, 0, voxelBasedMeasureGradient->nvox * voxelBasedMeasureGradient->nbyper); + memset(voxelBasedMeasureGradient->data, 0, voxelBasedMeasureGradient.totalBytes()); } /* *************************************************************** */ diff --git a/reg-lib/DefContent.h b/reg-lib/DefContent.h index a5ccab6f..d6ee0313 100644 --- a/reg-lib/DefContent.h +++ b/reg-lib/DefContent.h @@ -5,18 +5,17 @@ class DefContent: public virtual Content { public: DefContent() = delete; - DefContent(nifti_image *referenceIn, - nifti_image *floatingIn, - nifti_image *localWeightSimIn = nullptr, + DefContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, + NiftiImage *localWeightSimIn = nullptr, int *referenceMaskIn = nullptr, mat44 *transformationMatrixIn = nullptr, size_t bytesIn = sizeof(float)); - virtual ~DefContent(); // Getters - virtual nifti_image* GetLocalWeightSim() { return localWeightSim; } - virtual nifti_image* GetVoxelBasedMeasureGradient() { return voxelBasedMeasureGradient; } - virtual nifti_image* GetWarpedGradient() { return warpedGradient; } + virtual NiftiImage& GetLocalWeightSim() { return localWeightSim; } + virtual NiftiImage& GetVoxelBasedMeasureGradient() { return voxelBasedMeasureGradient; } + virtual NiftiImage& GetWarpedGradient() { return warpedGradient; } // Methods for transferring data from nifti to device virtual void UpdateVoxelBasedMeasureGradient() {} @@ -26,15 +25,10 @@ class DefContent: public virtual Content { virtual void ZeroVoxelBasedMeasureGradient(); protected: - nifti_image *localWeightSim = nullptr; - nifti_image *voxelBasedMeasureGradient = nullptr; - nifti_image *warpedGradient = nullptr; + NiftiImage localWeightSim; + NiftiImage voxelBasedMeasureGradient; + NiftiImage warpedGradient; private: - void AllocateLocalWeightSim(nifti_image*); - void DeallocateLocalWeightSim(); - void AllocateVoxelBasedMeasureGradient(); - void DeallocateVoxelBasedMeasureGradient(); - void AllocateWarpedGradient(); - void DeallocateWarpedGradient(); + void AllocateLocalWeightSim(NiftiImage&); }; diff --git a/reg-lib/DefContentCreator.h b/reg-lib/DefContentCreator.h index dce3ba86..e8302616 100644 --- a/reg-lib/DefContentCreator.h +++ b/reg-lib/DefContentCreator.h @@ -5,9 +5,9 @@ class DefContentCreator: public ContentCreator { public: - virtual DefContent* Create(nifti_image *reference, - nifti_image *floating, - nifti_image *localWeightSim = nullptr, + virtual DefContent* Create(NiftiImage& reference, + NiftiImage& floating, + NiftiImage *localWeightSim = nullptr, int *referenceMask = nullptr, mat44 *transformationMatrix = nullptr, size_t bytes = sizeof(float)) { diff --git a/reg-lib/F3d2ContentCreator.h b/reg-lib/F3d2ContentCreator.h index 106b5ede..6141d6bb 100644 --- a/reg-lib/F3d2ContentCreator.h +++ b/reg-lib/F3d2ContentCreator.h @@ -5,11 +5,11 @@ class F3d2ContentCreator: public ContentCreator { public: - virtual std::pair Create(nifti_image *reference, - nifti_image *floating, - nifti_image *controlPointGrid, - nifti_image *controlPointGridBw, - nifti_image *localWeightSim = nullptr, + virtual std::pair Create(NiftiImage& reference, + NiftiImage& floating, + NiftiImage& controlPointGrid, + NiftiImage& controlPointGridBw, + NiftiImage *localWeightSim = nullptr, int *referenceMask = nullptr, int *floatingMask = nullptr, mat44 *transformationMatrix = nullptr, diff --git a/reg-lib/F3dContent.cpp b/reg-lib/F3dContent.cpp index 6dee6030..e0a36c6d 100644 --- a/reg-lib/F3dContent.cpp +++ b/reg-lib/F3dContent.cpp @@ -1,37 +1,22 @@ #include "F3dContent.h" /* *************************************************************** */ -F3dContent::F3dContent(nifti_image *referenceIn, - nifti_image *floatingIn, - nifti_image *controlPointGridIn, - nifti_image *localWeightSimIn, +F3dContent::F3dContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, + NiftiImage& controlPointGridIn, + NiftiImage *localWeightSimIn, int *referenceMaskIn, mat44 *transformationMatrixIn, size_t bytesIn): DefContent(referenceIn, floatingIn, localWeightSimIn, referenceMaskIn, transformationMatrixIn, bytesIn), Content(referenceIn, floatingIn, referenceMaskIn, transformationMatrixIn, bytesIn), - controlPointGrid(controlPointGridIn) { + controlPointGrid(NiftiImage(controlPointGridIn, NiftiImage::Copy::Acquire)) { if (!controlPointGridIn) NR_FATAL_ERROR("controlPointGridIn can't be nullptr"); - AllocateTransformationGradient(); -} -/* *************************************************************** */ -F3dContent::~F3dContent() { - DeallocateTransformationGradient(); -} -/* *************************************************************** */ -void F3dContent::AllocateTransformationGradient() { - transformationGradient = nifti_dup(*controlPointGrid, false); -} -/* *************************************************************** */ -void F3dContent::DeallocateTransformationGradient() { - if (transformationGradient != nullptr) { - nifti_image_free(transformationGradient); - transformationGradient = nullptr; - } + transformationGradient = NiftiImage(controlPointGrid, NiftiImage::Copy::ImageInfoAndAllocData); } /* *************************************************************** */ void F3dContent::ZeroTransformationGradient() { - memset(transformationGradient->data, 0, transformationGradient->nvox * transformationGradient->nbyper); + memset(transformationGradient->data, 0, transformationGradient.totalBytes()); } /* *************************************************************** */ diff --git a/reg-lib/F3dContent.h b/reg-lib/F3dContent.h index f09157c0..c36f5634 100644 --- a/reg-lib/F3dContent.h +++ b/reg-lib/F3dContent.h @@ -5,18 +5,17 @@ class F3dContent: public virtual DefContent { public: F3dContent() = delete; - F3dContent(nifti_image *referenceIn, - nifti_image *floatingIn, - nifti_image *controlPointGridIn, - nifti_image *localWeightSimIn = nullptr, + F3dContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, + NiftiImage& controlPointGridIn, + NiftiImage *localWeightSimIn = nullptr, int *referenceMaskIn = nullptr, mat44 *transformationMatrixIn = nullptr, size_t bytesIn = sizeof(float)); - virtual ~F3dContent(); // Getters - virtual nifti_image* GetControlPointGrid() { return controlPointGrid; } - virtual nifti_image* GetTransformationGradient() { return transformationGradient; } + virtual NiftiImage& GetControlPointGrid() { return controlPointGrid; } + virtual NiftiImage& GetTransformationGradient() { return transformationGradient; } // Methods for transferring data from nifti to device virtual void UpdateControlPointGrid() {} @@ -26,10 +25,6 @@ class F3dContent: public virtual DefContent { virtual void ZeroTransformationGradient(); protected: - nifti_image *controlPointGrid = nullptr; - nifti_image *transformationGradient = nullptr; - -private: - void AllocateTransformationGradient(); - void DeallocateTransformationGradient(); + NiftiImage controlPointGrid; + NiftiImage transformationGradient; }; \ No newline at end of file diff --git a/reg-lib/F3dContentCreator.h b/reg-lib/F3dContentCreator.h index d57657b0..2ee586dc 100644 --- a/reg-lib/F3dContentCreator.h +++ b/reg-lib/F3dContentCreator.h @@ -5,10 +5,10 @@ class F3dContentCreator: public ContentCreator { public: - virtual F3dContent* Create(nifti_image *reference, - nifti_image *floating, - nifti_image *controlPointGrid, - nifti_image *localWeightSim = nullptr, + virtual F3dContent* Create(NiftiImage& reference, + NiftiImage& floating, + NiftiImage& controlPointGrid, + NiftiImage *localWeightSim = nullptr, int *referenceMask = nullptr, mat44 *transformationMatrix = nullptr, size_t bytes = sizeof(float)) { diff --git a/reg-lib/MeasureCreator.cpp b/reg-lib/MeasureCreator.cpp index 6ff56f71..473f8359 100644 --- a/reg-lib/MeasureCreator.cpp +++ b/reg-lib/MeasureCreator.cpp @@ -38,8 +38,8 @@ void MeasureCreator::Initialise(reg_measure& measure, DefContent& con, DefConten con.GetVoxelBasedMeasureGradient(), con.GetLocalWeightSim(), conBw ? conBw->GetReferenceMask() : nullptr, - conBw ? conBw->GetWarped() : nullptr, - conBw ? conBw->GetWarpedGradient() : nullptr, - conBw ? conBw->GetVoxelBasedMeasureGradient() : nullptr); + conBw ? static_cast(conBw->GetWarped()) : nullptr, + conBw ? static_cast(conBw->GetWarpedGradient()) : nullptr, + conBw ? static_cast(conBw->GetVoxelBasedMeasureGradient()) : nullptr); } /* *************************************************************** */ diff --git a/reg-lib/Platform.cpp b/reg-lib/Platform.cpp index 482089fa..8a88df01 100755 --- a/reg-lib/Platform.cpp +++ b/reg-lib/Platform.cpp @@ -119,27 +119,27 @@ Optimiser* Platform::CreateOptimiser(F3dContent& con, F3dContent *conBw) const { Optimiser *optimiser; nifti_image *controlPointGrid = con.F3dContent::GetControlPointGrid(); - nifti_image *controlPointGridBw = conBw ? conBw->F3dContent::GetControlPointGrid() : nullptr; + nifti_image *controlPointGridBw = conBw ? static_cast(conBw->F3dContent::GetControlPointGrid()) : nullptr; Type *controlPointGridData, *transformationGradientData; Type *controlPointGridDataBw = nullptr, *transformationGradientDataBw = nullptr; if (platformType == PlatformType::Cpu) { optimiser = useConjGradient ? new ConjugateGradient() : new Optimiser(); - controlPointGridData = (Type*)controlPointGrid->data; - transformationGradientData = (Type*)con.GetTransformationGradient()->data; + controlPointGridData = static_cast(controlPointGrid->data); + transformationGradientData = static_cast(con.GetTransformationGradient()->data); if (conBw) { - controlPointGridDataBw = (Type*)controlPointGridBw->data; - transformationGradientDataBw = (Type*)conBw->GetTransformationGradient()->data; + controlPointGridDataBw = static_cast(controlPointGridBw->data); + transformationGradientDataBw = static_cast(conBw->GetTransformationGradient()->data); } } #ifdef USE_CUDA else if (platformType == PlatformType::Cuda) { optimiser = dynamic_cast*>(useConjGradient ? new CudaConjugateGradient() : new CudaOptimiser()); - controlPointGridData = (Type*)dynamic_cast(con).GetControlPointGridCuda(); - transformationGradientData = (Type*)dynamic_cast(con).GetTransformationGradientCuda(); + controlPointGridData = reinterpret_cast(dynamic_cast(con).GetControlPointGridCuda()); + transformationGradientData = reinterpret_cast(dynamic_cast(con).GetTransformationGradientCuda()); if (conBw) { - controlPointGridDataBw = (Type*)dynamic_cast(conBw)->GetControlPointGridCuda(); - transformationGradientDataBw = (Type*)dynamic_cast(conBw)->GetTransformationGradientCuda(); + controlPointGridDataBw = reinterpret_cast(dynamic_cast(conBw)->GetControlPointGridCuda()); + transformationGradientDataBw = reinterpret_cast(dynamic_cast(conBw)->GetTransformationGradientCuda()); } } #endif diff --git a/reg-lib/_reg_aladin.cpp b/reg-lib/_reg_aladin.cpp index 032aeb97..959c7470 100644 --- a/reg-lib/_reg_aladin.cpp +++ b/reg-lib/_reg_aladin.cpp @@ -287,8 +287,8 @@ void reg_aladin::UpdateTransformationMatrix(int type) { } /* *************************************************************** */ template -void reg_aladin::InitAladinContent(nifti_image *ref, - nifti_image *flo, +void reg_aladin::InitAladinContent(NiftiImage& ref, + NiftiImage& flo, int *mask, mat44 *transMat, size_t bytes, @@ -384,7 +384,7 @@ NiftiImage reg_aladin::GetFinalWarpedImage() { reg_aladin::GetWarpedImage(3, this->warpedPaddingValue); // cubic spline interpolation - NiftiImage warpedImage(this->con->GetWarped(), NiftiImage::Copy::Image); + NiftiImage warpedImage(this->con->GetWarped()); warpedImage->cal_min = this->inputFloating->cal_min; warpedImage->cal_max = this->inputFloating->cal_max; warpedImage->scl_slope = this->inputFloating->scl_slope; @@ -397,8 +397,8 @@ NiftiImage reg_aladin::GetFinalWarpedImage() { /* *************************************************************** */ template void reg_aladin::DebugPrintLevelInfoStart() { - const nifti_image *ref = this->con->Content::GetReference(); - const nifti_image *flo = this->con->Content::GetFloating(); + const NiftiImage& ref = this->con->Content::GetReference(); + const NiftiImage& flo = this->con->Content::GetFloating(); NR_VERBOSE("Current level " << this->currentLevel + 1 << " / " << this->numberOfLevels); NR_VERBOSE("Reference image size:\t" << ref->nx << "x" << ref->ny << "x" << ref->nz << " voxels\t" << ref->dx << "x" << ref->dy << "x" << ref->dz << " mm"); diff --git a/reg-lib/_reg_aladin.h b/reg-lib/_reg_aladin.h index 9096688d..17f544b1 100644 --- a/reg-lib/_reg_aladin.h +++ b/reg-lib/_reg_aladin.h @@ -117,8 +117,8 @@ class reg_aladin { void *paramsProgressCallback; //platform factory methods - virtual void InitAladinContent(nifti_image *ref, - nifti_image *flo, + virtual void InitAladinContent(NiftiImage& ref, + NiftiImage& flo, int *mask, mat44 *transMat, size_t bytes, diff --git a/reg-lib/_reg_aladin_sym.cpp b/reg-lib/_reg_aladin_sym.cpp index 62cdd753..2fa2ff18 100644 --- a/reg-lib/_reg_aladin_sym.cpp +++ b/reg-lib/_reg_aladin_sym.cpp @@ -157,8 +157,8 @@ void reg_aladin_sym::UpdateTransformationMatrix(int type) { } /* *************************************************************** */ template -void reg_aladin_sym::InitAladinContent(nifti_image *ref, - nifti_image *flo, +void reg_aladin_sym::InitAladinContent(NiftiImage& ref, + NiftiImage& flo, int *mask, mat44 *transMat, size_t bytes, @@ -203,8 +203,8 @@ void reg_aladin_sym::DeallocateKernels() { /* *************************************************************** */ template void reg_aladin_sym::DebugPrintLevelInfoStart() { - const nifti_image *ref = this->con->Content::GetReference(); - const nifti_image *flo = this->con->Content::GetFloating(); + const NiftiImage& ref = this->con->Content::GetReference(); + const NiftiImage& flo = this->con->Content::GetFloating(); NR_VERBOSE("Current level " << this->currentLevel + 1 << " / " << this->numberOfLevels); NR_VERBOSE("Reference image size:\t" << ref->nx << "x" << ref->ny << "x" << ref->nz << " voxels\t" << ref->dx << "x" << ref->dy << "x" << ref->dz << " mm"); diff --git a/reg-lib/_reg_aladin_sym.h b/reg-lib/_reg_aladin_sym.h index 6da18e76..028c5cb1 100644 --- a/reg-lib/_reg_aladin_sym.h +++ b/reg-lib/_reg_aladin_sym.h @@ -21,8 +21,8 @@ class reg_aladin_sym: public reg_aladin { unique_ptr backCon; unique_ptr bAffineTransformation3DKernel, bConvolutionKernel, bBlockMatchingKernel, bLtsKernel, bResamplingKernel; - virtual void InitAladinContent(nifti_image *ref, - nifti_image *flo, + virtual void InitAladinContent(NiftiImage& ref, + NiftiImage& flo, int *mask, mat44 *transMat, size_t bytes, diff --git a/reg-lib/_reg_f3d.cpp b/reg-lib/_reg_f3d.cpp index 1f005525..8aa5f57b 100644 --- a/reg-lib/_reg_f3d.cpp +++ b/reg-lib/_reg_f3d.cpp @@ -82,29 +82,29 @@ void reg_f3d::SetSpacing(unsigned i, T s) { } /* *************************************************************** */ template -void reg_f3d::InitContent(nifti_image *reference, nifti_image *floating, int *mask) { +void reg_f3d::InitContent(NiftiImage& reference, NiftiImage& floating, int *mask) { unique_ptr contentCreator{ dynamic_cast(this->platform->CreateContentCreator(ContentType::F3d)) }; - this->con.reset(contentCreator->Create(reference, floating, controlPointGrid, this->localWeightSimInput, mask, this->affineTransformation.get(), sizeof(T))); + this->con.reset(contentCreator->Create(reference, floating, controlPointGrid, &this->localWeightSimInput, mask, this->affineTransformation.get(), sizeof(T))); this->compute.reset(this->platform->CreateCompute(*this->con)); } /* *************************************************************** */ template T reg_f3d::InitCurrentLevel(int currentLevel) { // Set the current input images - nifti_image *reference, *floating; + NiftiImage reference, floating; int *mask; if (currentLevel < 0) { // Settings for GetWarpedImage() // Use CPU for warping since CUDA isn't supporting Cubic interpolation // TODO Remove this when CUDA supports Cubic interpolation this->SetPlatformType(PlatformType::Cpu); - reference = this->inputReference; - floating = this->inputFloating; + reference = NiftiImage(this->inputReference, NiftiImage::Copy::Acquire); + floating = NiftiImage(this->inputFloating, NiftiImage::Copy::Acquire); mask = nullptr; } else { const int index = this->usePyramid ? currentLevel : 0; - reference = this->referencePyramid[index]; - floating = this->floatingPyramid[index]; + reference = NiftiImage(this->referencePyramid[index], NiftiImage::Copy::Acquire); + floating = NiftiImage(this->floatingPyramid[index], NiftiImage::Copy::Acquire); mask = this->maskPyramid[index].get(); } @@ -411,8 +411,8 @@ T reg_f3d::NormaliseGradient() { /* *************************************************************** */ template void reg_f3d::DisplayCurrentLevelParameters(int currentLevel) { - const nifti_image *reference = this->con->Content::GetReference(); - const nifti_image *floating = this->con->Content::GetFloating(); + const NiftiImage& reference = this->con->Content::GetReference(); + const NiftiImage& floating = this->con->Content::GetFloating(); NR_VERBOSE("Current level: " << currentLevel + 1 << " / " << this->levelNumber); NR_VERBOSE("Maximum iteration number: " << this->maxIterationNumber); NR_VERBOSE("Current reference image"); @@ -507,11 +507,11 @@ vector reg_f3d::GetWarpedImage() { InitCurrentLevel(-1); this->WarpFloatingImage(3); // cubic spline interpolation - NiftiImage warpedImage = NiftiImage(this->con->GetWarped(), NiftiImage::Copy::Image); + NiftiImage warpedImage = std::move(this->con->GetWarped()); DeinitCurrentLevel(-1); NR_FUNC_CALLED(); - return { warpedImage }; + return { std::move(warpedImage) }; } /* *************************************************************** */ template diff --git a/reg-lib/_reg_f3d.h b/reg-lib/_reg_f3d.h index a7a793ca..e13fbe1f 100644 --- a/reg-lib/_reg_f3d.h +++ b/reg-lib/_reg_f3d.h @@ -33,7 +33,7 @@ class reg_f3d: public reg_base { double bestWBE; double bestWLE; - void InitContent(nifti_image*, nifti_image*, int*); + void InitContent(NiftiImage&, NiftiImage&, int*); virtual T InitCurrentLevel(int) override; virtual void DeinitCurrentLevel(int) override; virtual T NormaliseGradient() override; diff --git a/reg-lib/_reg_f3d2.cpp b/reg-lib/_reg_f3d2.cpp index eaa7a6f0..47dae410 100644 --- a/reg-lib/_reg_f3d2.cpp +++ b/reg-lib/_reg_f3d2.cpp @@ -38,10 +38,10 @@ void reg_f3d2::SetInverseConsistencyWeight(T w) { } /* *************************************************************** */ template -void reg_f3d2::InitContent(nifti_image *reference, nifti_image *floating, int *referenceMask, int *floatingMask) { +void reg_f3d2::InitContent(NiftiImage& reference, NiftiImage& floating, int *referenceMask, int *floatingMask) { unique_ptr contentCreator{ dynamic_cast(this->platform->CreateContentCreator(ContentType::F3d2)) }; auto&& [con, conBw] = contentCreator->Create(reference, floating, this->controlPointGrid, controlPointGridBw, - this->localWeightSimInput, referenceMask, floatingMask, + &this->localWeightSimInput, referenceMask, floatingMask, this->affineTransformation.get(), affineTransformationBw.get(), sizeof(T)); this->con.reset(con); this->conBw.reset(conBw); @@ -52,21 +52,21 @@ void reg_f3d2::InitContent(nifti_image *reference, nifti_image *floating, int template T reg_f3d2::InitCurrentLevel(int currentLevel) { // Set the current input images - nifti_image *reference, *floating; + NiftiImage reference, floating; int *referenceMask, *floatingMask; if (currentLevel < 0) { // Settings for GetWarpedImage() // Use CPU for warping since CUDA isn't supporting Cubic interpolation // TODO Remove this when CUDA supports Cubic interpolation this->SetPlatformType(PlatformType::Cpu); - reference = this->inputReference; - floating = this->inputFloating; + reference = NiftiImage(this->inputReference, NiftiImage::Copy::Acquire); + floating = NiftiImage(this->inputFloating, NiftiImage::Copy::Acquire); referenceMask = nullptr; floatingMask = nullptr; } else { const int index = this->usePyramid ? currentLevel : 0; - reference = this->referencePyramid[index]; - floating = this->floatingPyramid[index]; + reference = NiftiImage(this->referencePyramid[index], NiftiImage::Copy::Acquire); + floating = NiftiImage(this->floatingPyramid[index], NiftiImage::Copy::Acquire); referenceMask = this->maskPyramid[index].get(); floatingMask = floatingMaskPyramid[index].get(); } @@ -666,15 +666,12 @@ vector reg_f3d2::GetWarpedImage() { WarpFloatingImage(3); // cubic spline interpolation F3dContent& con = dynamic_cast(*this->con); - vector warpedImage{ - NiftiImage(con.GetWarped(), NiftiImage::Copy::Image), - NiftiImage(conBw->GetWarped(), NiftiImage::Copy::Image) - }; + vector warpedImages{ std::move(con.GetWarped()), std::move(conBw->GetWarped()) }; DeinitCurrentLevel(-1); NR_FUNC_CALLED(); - return warpedImage; + return warpedImages; } /* *************************************************************** */ template class reg_f3d2; diff --git a/reg-lib/_reg_f3d2.h b/reg-lib/_reg_f3d2.h index c11c857e..12f83917 100644 --- a/reg-lib/_reg_f3d2.h +++ b/reg-lib/_reg_f3d2.h @@ -55,7 +55,7 @@ class reg_f3d2: public reg_f3d { virtual void PrintCurrentObjFunctionValue(T) override; virtual void UpdateBestObjFunctionValue() override; virtual double GetObjectiveFunctionValue() override; - void InitContent(nifti_image*, nifti_image*, int*, int*); + void InitContent(NiftiImage&, NiftiImage&, int*, int*); virtual T InitCurrentLevel(int) override; virtual void DeinitCurrentLevel(int) override; virtual void UpdateParameters(float) override; diff --git a/reg-lib/cl/ClAladinContent.cpp b/reg-lib/cl/ClAladinContent.cpp index 45fac34e..62456523 100644 --- a/reg-lib/cl/ClAladinContent.cpp +++ b/reg-lib/cl/ClAladinContent.cpp @@ -2,8 +2,8 @@ #include "_reg_tools.h" /* *************************************************************** */ -ClAladinContent::ClAladinContent(nifti_image *referenceIn, - nifti_image *floatingIn, +ClAladinContent::ClAladinContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, int *referenceMaskIn, mat44 *transformationMatrixIn, size_t bytesIn, @@ -36,29 +36,28 @@ void ClAladinContent::InitVars() { totalBlockClmem = nullptr; maskClmem = nullptr; - if (reference != nullptr && reference->nbyper != NIFTI_TYPE_FLOAT32) + if (reference && reference->nbyper != NIFTI_TYPE_FLOAT32) reg_tools_changeDatatype(reference); - if (floating != nullptr && floating->nbyper != NIFTI_TYPE_FLOAT32) { + if (floating && floating->nbyper != NIFTI_TYPE_FLOAT32) { reg_tools_changeDatatype(floating); - if (warped != nullptr) + if (warped) reg_tools_changeDatatype(warped); } sContext = &ClContextSingleton::GetInstance(); clContext = sContext->GetContext(); commandQueue = sContext->GetCommandQueue(); - //numBlocks = (blockMatchingParams != nullptr) ? blockMatchingParams->blockNumber[0] * blockMatchingParams->blockNumber[1] * blockMatchingParams->blockNumber[2] : 0; } /* *************************************************************** */ void ClAladinContent::AllocateClPtrs() { - if (warped != nullptr) { + if (warped) { warpedImageClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, warped->nvox * sizeof(float), warped->data, &errNum); sContext->CheckErrNum(errNum, "ClAladinContent::AllocateClPtrs failed to allocate memory (warpedImageClmem): "); } - if (deformationField != nullptr) { + if (deformationField) { deformationFieldClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, sizeof(float) * deformationField->nvox, deformationField->data, &errNum); sContext->CheckErrNum(errNum, "ClAladinContent::AllocateClPtrs failed to allocate memory (deformationFieldClmem): "); } - if (floating != nullptr) { + if (floating) { floatingImageClmem = clCreateBuffer(clContext, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * floating->nvox, floating->data, &errNum); sContext->CheckErrNum(errNum, "ClAladinContent::AllocateClPtrs failed to allocate memory (floating): "); @@ -68,7 +67,7 @@ void ClAladinContent::AllocateClPtrs() { sContext->CheckErrNum(errNum, "ClContent::AllocateClPtrs failed to allocate memory (floMatClmem): "); free(sourceIJKMatrix_h); } - if (reference != nullptr) { + if (reference) { referenceImageClmem = clCreateBuffer(clContext, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * reference->nvox, reference->data, &errNum); @@ -80,22 +79,22 @@ void ClAladinContent::AllocateClPtrs() { sContext->CheckErrNum(errNum, "ClContent::AllocateClPtrs failed to allocate memory (refMatClmem): "); free(targetMat); } - if (blockMatchingParams != nullptr) { - if (blockMatchingParams->referencePosition != nullptr) { + if (blockMatchingParams) { + if (blockMatchingParams->referencePosition) { //targetPositionClmem referencePositionClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim * sizeof(float), blockMatchingParams->referencePosition, &errNum); sContext->CheckErrNum(errNum, "ClContent::AllocateClPtrs failed to allocate memory (referencePositionClmem): "); } - if (blockMatchingParams->warpedPosition != nullptr) { + if (blockMatchingParams->warpedPosition) { //resultPositionClmem warpedPositionClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim * sizeof(float), blockMatchingParams->warpedPosition, &errNum); sContext->CheckErrNum(errNum, "ClContent::AllocateClPtrs failed to allocate memory (warpedPositionClmem): "); } - if (blockMatchingParams->totalBlock != nullptr) { + if (blockMatchingParams->totalBlock) { //totalBlockClmem totalBlockClmem = clCreateBuffer(clContext, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, blockMatchingParams->totalBlockNumber * sizeof(int), @@ -103,19 +102,19 @@ void ClAladinContent::AllocateClPtrs() { sContext->CheckErrNum(errNum, "ClContent::AllocateClPtrs failed to allocate memory (activeBlockClmem): "); } } - if (referenceMask != nullptr && reference != nullptr) { + if (referenceMask && reference) { maskClmem = clCreateBuffer(clContext, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, - NiftiImage::calcVoxelNumber(reference, 3) * sizeof(int), referenceMask, &errNum); + reference.nVoxelsPerVolume() * sizeof(int), referenceMask, &errNum); sContext->CheckErrNum(errNum, "ClContent::AllocateClPtrs failed to allocate memory (clCreateBuffer): "); } } /* *************************************************************** */ -nifti_image* ClAladinContent::GetWarped() { +NiftiImage& ClAladinContent::GetWarped() { DownloadImage(warped, warpedImageClmem, warped->datatype); return warped; } /* *************************************************************** */ -nifti_image* ClAladinContent::GetDeformationField() { +NiftiImage& ClAladinContent::GetDeformationField() { errNum = clEnqueueReadBuffer(commandQueue, deformationFieldClmem, CL_TRUE, 0, deformationField->nvox * sizeof(float), deformationField->data, 0, nullptr, nullptr); //CLCONTEXT sContext->CheckErrNum(errNum, "Get: failed deformationField: "); return deformationField; @@ -133,48 +132,48 @@ void ClAladinContent::SetTransformationMatrix(mat44 *transformationMatrixIn) { AladinContent::SetTransformationMatrix(transformationMatrixIn); } /* *************************************************************** */ -void ClAladinContent::SetDeformationField(nifti_image *deformationFieldIn) { - if (deformationField != nullptr) +void ClAladinContent::SetDeformationField(NiftiImage&& deformationFieldIn) { + if (deformationField) clReleaseMemObject(deformationFieldClmem); - AladinContent::SetDeformationField(deformationFieldIn); + AladinContent::SetDeformationField(std::move(deformationFieldIn)); deformationFieldClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, deformationField->nvox * sizeof(float), deformationField->data, &errNum); sContext->CheckErrNum(errNum, "ClAladinContent::SetDeformationField failed to allocate memory (deformationFieldClmem): "); } /* *************************************************************** */ void ClAladinContent::SetReferenceMask(int *referenceMaskIn) { - if (referenceMask != nullptr) + if (referenceMask) clReleaseMemObject(maskClmem); AladinContent::SetReferenceMask(referenceMaskIn); maskClmem = clCreateBuffer(clContext, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, reference->nvox * sizeof(int), referenceMask, &errNum); sContext->CheckErrNum(errNum, "ClAladinContent::SetReferenceMask failed to allocate memory (maskClmem): "); } /* *************************************************************** */ -void ClAladinContent::SetWarped(nifti_image *warpedIn) { +void ClAladinContent::SetWarped(NiftiImage&& warpedIn) { if (warpedIn->nbyper != NIFTI_TYPE_FLOAT32) - reg_tools_changeDatatype(warpedIn); - if (warped != nullptr) + warpedIn.changeDatatype(NIFTI_TYPE_FLOAT32); + if (warped) clReleaseMemObject(warpedImageClmem); - AladinContent::SetWarped(warpedIn); - warpedImageClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, warpedIn->nvox * sizeof(float), warpedIn->data, &errNum); + AladinContent::SetWarped(std::move(warpedIn)); + warpedImageClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, warped->nvox * sizeof(float), warped->data, &errNum); sContext->CheckErrNum(errNum, "ClAladinContent::SetWarped failed to allocate memory (warpedImageClmem): "); } /* *************************************************************** */ void ClAladinContent::SetBlockMatchingParams(_reg_blockMatchingParam* bmp) { AladinContent::SetBlockMatchingParams(bmp); - if (blockMatchingParams->referencePosition != nullptr) { + if (blockMatchingParams->referencePosition) { clReleaseMemObject(referencePositionClmem); //referencePositionClmem referencePositionClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim * sizeof(float), blockMatchingParams->referencePosition, &errNum); sContext->CheckErrNum(errNum, "ClAladinContent::SetBlockMatchingParams failed to allocate memory (referencePositionClmem): "); } - if (blockMatchingParams->warpedPosition != nullptr) { + if (blockMatchingParams->warpedPosition) { clReleaseMemObject(warpedPositionClmem); //warpedPositionClmem warpedPositionClmem = clCreateBuffer(clContext, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim * sizeof(float), blockMatchingParams->warpedPosition, &errNum); sContext->CheckErrNum(errNum, "ClAladinContent::SetBlockMatchingParams failed to allocate memory (warpedPositionClmem): "); } - if (blockMatchingParams->totalBlock != nullptr) { + if (blockMatchingParams->totalBlock) { clReleaseMemObject(totalBlockClmem); //totalBlockClmem totalBlockClmem = clCreateBuffer(clContext, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, blockMatchingParams->totalBlockNumber * sizeof(int), blockMatchingParams->totalBlock, &errNum); @@ -222,8 +221,7 @@ cl_mem ClAladinContent::GetFloMatClmem() { return floMatClmem; } /* *************************************************************** */ -template -void ClAladinContent::FillImageData(nifti_image *image, cl_mem memoryObject, int datatype) { +void ClAladinContent::DownloadImage(NiftiImage& image, cl_mem memoryObject, int datatype) { const size_t size = image->nvox; unique_ptr buffer(new float[size]); @@ -231,62 +229,33 @@ void ClAladinContent::FillImageData(nifti_image *image, cl_mem memoryObject, int size * sizeof(float), buffer.get(), 0, nullptr, nullptr); sContext->CheckErrNum(errNum, "Error reading warped buffer."); - free(image->data); - image->datatype = datatype; - image->nbyper = sizeof(DataType); - image->data = malloc(size * image->nbyper); - DataType *data = static_cast(image->data); - for (size_t i = 0; i < size; ++i) - data[i] = static_cast(NiftiImage::clampData(image, buffer[i])); -} -/* *************************************************************** */ -void ClAladinContent::DownloadImage(nifti_image *image, cl_mem memoryObject, int datatype) { - switch (datatype) { - case NIFTI_TYPE_FLOAT32: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_FLOAT64: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT8: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT8: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT16: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT16: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT32: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT32: - FillImageData(image, memoryObject, datatype); - break; - default: - NR_FATAL_ERROR("Unsupported type"); - } + std::visit([&](auto&& dataType) { + using DataType = std::decay_t; + image->datatype = datatype; + image->nbyper = sizeof(DataType); + image.realloc(); + DataType *data = static_cast(image->data); + for (size_t i = 0; i < size; ++i) + data[i] = static_cast(image.clampData(buffer[i])); + }, image.getDataType()); } /* *************************************************************** */ void ClAladinContent::FreeClPtrs() { - if (reference != nullptr) { + if (reference) { clReleaseMemObject(referenceImageClmem); clReleaseMemObject(refMatClmem); } - if (floating != nullptr) { + if (floating) { clReleaseMemObject(floatingImageClmem); clReleaseMemObject(floMatClmem); } - if (warped != nullptr) + if (warped) clReleaseMemObject(warpedImageClmem); - if (deformationField != nullptr) + if (deformationField) clReleaseMemObject(deformationFieldClmem); - if (referenceMask != nullptr) + if (referenceMask) clReleaseMemObject(maskClmem); - if (blockMatchingParams != nullptr) { + if (blockMatchingParams) { clReleaseMemObject(totalBlockClmem); clReleaseMemObject(referencePositionClmem); clReleaseMemObject(warpedPositionClmem); diff --git a/reg-lib/cl/ClAladinContent.h b/reg-lib/cl/ClAladinContent.h index 7de5039b..3a76552b 100644 --- a/reg-lib/cl/ClAladinContent.h +++ b/reg-lib/cl/ClAladinContent.h @@ -12,8 +12,8 @@ class ClAladinContent: public AladinContent { public: //constructors - ClAladinContent(nifti_image *referenceIn, - nifti_image *floatingIn, + ClAladinContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, int *referenceMaskIn = nullptr, mat44 *transformationMatrixIn = nullptr, size_t bytesIn = sizeof(float), @@ -38,8 +38,8 @@ class ClAladinContent: public AladinContent { // CPU getters with data downloaded from device virtual _reg_blockMatchingParam* GetBlockMatchingParams() override; - virtual nifti_image* GetDeformationField() override; - virtual nifti_image* GetWarped() override; + virtual NiftiImage& GetDeformationField() override; + virtual NiftiImage& GetWarped() override; private: void InitVars(); @@ -62,8 +62,7 @@ class ClAladinContent: public AladinContent { cl_mem refMatClmem; cl_mem floMatClmem; - template void FillImageData(nifti_image *image, cl_mem memoryObject, int datatype); - void DownloadImage(nifti_image *image, cl_mem memoryObject, int datatype); + void DownloadImage(NiftiImage& image, cl_mem memoryObject, int datatype); #ifdef NR_TESTING public: @@ -72,8 +71,8 @@ class ClAladinContent: public AladinContent { #endif // Functions for testing virtual void SetTransformationMatrix(mat44 *transformationMatrixIn) override; - virtual void SetWarped(nifti_image *warpedIn) override; - virtual void SetDeformationField(nifti_image *deformationFieldIn) override; + virtual void SetWarped(NiftiImage&& warpedIn) override; + virtual void SetDeformationField(NiftiImage&& deformationFieldIn) override; virtual void SetReferenceMask(int *referenceMaskIn) override; virtual void SetBlockMatchingParams(_reg_blockMatchingParam *bmp) override; }; diff --git a/reg-lib/cl/ClAladinContentCreator.h b/reg-lib/cl/ClAladinContentCreator.h index 84442142..ed688de7 100644 --- a/reg-lib/cl/ClAladinContentCreator.h +++ b/reg-lib/cl/ClAladinContentCreator.h @@ -5,8 +5,8 @@ class ClAladinContentCreator: public AladinContentCreator { 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), diff --git a/reg-lib/cl/ClBlockMatchingKernel.cpp b/reg-lib/cl/ClBlockMatchingKernel.cpp index 06002aa9..1553d36e 100644 --- a/reg-lib/cl/ClBlockMatchingKernel.cpp +++ b/reg-lib/cl/ClBlockMatchingKernel.cpp @@ -57,7 +57,6 @@ ClBlockMatchingKernel::ClBlockMatchingKernel(Content *conIn) : BlockMatchingKern //get cpu ptrs reference = con->AladinContent::GetReference(); params = con->AladinContent::GetBlockMatchingParams(); - } /* *************************************************************** */ void ClBlockMatchingKernel::Calculate() { diff --git a/reg-lib/cpu/CpuResampleImageKernel.cpp b/reg-lib/cpu/CpuResampleImageKernel.cpp index 1544e9d5..d723e023 100644 --- a/reg-lib/cpu/CpuResampleImageKernel.cpp +++ b/reg-lib/cpu/CpuResampleImageKernel.cpp @@ -13,7 +13,7 @@ CpuResampleImageKernel::CpuResampleImageKernel(Content *conIn) : ResampleImageKe void CpuResampleImageKernel::Calculate(int interp, float paddingValue, bool *dtiTimePoint, - mat33 * jacMat) { + mat33 *jacMat) { reg_resampleImage(floatingImage, warpedImage, deformationField, diff --git a/reg-lib/cpu/_reg_localTrans.cpp b/reg-lib/cpu/_reg_localTrans.cpp index d070bee1..418b310d 100755 --- a/reg-lib/cpu/_reg_localTrans.cpp +++ b/reg-lib/cpu/_reg_localTrans.cpp @@ -3613,15 +3613,14 @@ void reg_spline_getDefFieldFromVelocityGrid(nifti_image *velocityFieldGrid, } else NR_FATAL_ERROR("The provided input image is not a spline parametrised transformation"); } /* *************************************************************** */ -void reg_spline_getIntermediateDefFieldFromVelGrid(nifti_image *velocityFieldGrid, - nifti_image **deformationField) { +void reg_spline_getIntermediateDefFieldFromVelGrid(NiftiImage& velocityFieldGrid, + NiftiImage deformationFields[]) { // Check if the velocity field is actually a velocity field if (velocityFieldGrid->intent_p1 == SPLINE_VEL_GRID) { // Create an image to store the flow field - nifti_image *flowField = nifti_dup(*deformationField[0], false); + NiftiImage flowField(deformationFields[0], NiftiImage::Copy::ImageInfoAndAllocData); flowField->intent_code = NIFTI_INTENT_VECTOR; - memset(flowField->intent_name, 0, 16); - strcpy(flowField->intent_name, "NREG_TRANS"); + flowField.setIntentName("NREG_TRANS"s); flowField->intent_p1 = DEF_VEL_FIELD; flowField->intent_p2 = velocityFieldGrid->intent_p2; if (velocityFieldGrid->num_ext > 0 && flowField->ext_list == nullptr) @@ -3630,11 +3629,11 @@ void reg_spline_getIntermediateDefFieldFromVelGrid(nifti_image *velocityFieldGri // Generate the velocity field reg_spline_getFlowFieldFromVelocityGrid(velocityFieldGrid, flowField); // Remove the affine component from the flow field - nifti_image *affineOnly = nullptr; + NiftiImage affineOnly; if (flowField->num_ext > 0) { if (flowField->ext_list[0].edata != nullptr) { // Create a field that contains the affine component only - affineOnly = nifti_dup(*deformationField[0], false); + affineOnly = NiftiImage(deformationFields[0], NiftiImage::Copy::ImageInfoAndAllocData); reg_affine_getDeformationField(reinterpret_cast(flowField->ext_list[0].edata), affineOnly, false); reg_tools_subtractImageFromImage(flowField, affineOnly, flowField); } @@ -3647,45 +3646,38 @@ void reg_spline_getIntermediateDefFieldFromVelGrid(nifti_image *velocityFieldGri float scalingValue = pow(2.0f, std::abs((float)squaringNumber)); if (velocityFieldGrid->intent_p2 < 0) // backward deformation field is scaled down - reg_tools_divideValueToImage(flowField, deformationField[0], -scalingValue); + reg_tools_divideValueToImage(flowField, deformationFields[0], -scalingValue); else // forward deformation field is scaled down - reg_tools_divideValueToImage(flowField, deformationField[0], scalingValue); - - // Deallocate the allocated flow field - nifti_image_free(flowField); - flowField = nullptr; + reg_tools_divideValueToImage(flowField, deformationFields[0], scalingValue); // Conversion from displacement to deformation - reg_getDeformationFromDisplacement(deformationField[0]); + reg_getDeformationFromDisplacement(deformationFields[0]); // The deformation field is squared for (unsigned short i = 0; i < squaringNumber; ++i) { // The computed scaled deformation field is copied over - memcpy(deformationField[i + 1]->data, deformationField[i]->data, - deformationField[i]->nvox * deformationField[i]->nbyper); + deformationFields[i + 1].copyData(deformationFields[i]); // The deformation field is applied to itself - reg_defField_compose(deformationField[i], // to apply - deformationField[i + 1], // to update + reg_defField_compose(deformationFields[i], // to apply + deformationFields[i + 1], // to update nullptr); NR_DEBUG("Squaring (composition) step " << i + 1 << "/" << squaringNumber); } // The affine conponent of the transformation is restored - if (affineOnly != nullptr) { + if (affineOnly) { for (unsigned short i = 0; i <= squaringNumber; ++i) { - reg_getDisplacementFromDeformation(deformationField[i]); - reg_tools_addImageToImage(deformationField[i], affineOnly, deformationField[i]); - deformationField[i]->intent_p1 = DEF_FIELD; - deformationField[i]->intent_p2 = 0; + reg_getDisplacementFromDeformation(deformationFields[i]); + reg_tools_addImageToImage(deformationFields[i], affineOnly, deformationFields[i]); + deformationFields[i]->intent_p1 = DEF_FIELD; + deformationFields[i]->intent_p2 = 0; } - nifti_image_free(affineOnly); - affineOnly = nullptr; } // If required an affine component is composed if (velocityFieldGrid->num_ext > 1) { for (unsigned short i = 0; i <= squaringNumber; ++i) { reg_affine_getDeformationField(reinterpret_cast(velocityFieldGrid->ext_list[1].edata), - deformationField[i], + deformationFields[i], true); } } diff --git a/reg-lib/cpu/_reg_localTrans.h b/reg-lib/cpu/_reg_localTrans.h index c2a06195..955a704e 100755 --- a/reg-lib/cpu/_reg_localTrans.h +++ b/reg-lib/cpu/_reg_localTrans.h @@ -166,8 +166,8 @@ void reg_spline_getDefFieldFromVelocityGrid(nifti_image *velocityFieldGrid, nifti_image *deformationFieldImage, const bool updateStepNumber); /* *************************************************************** */ -void reg_spline_getIntermediateDefFieldFromVelGrid(nifti_image *velocityFieldGrid, - nifti_image **deformationFieldImage); +void reg_spline_getIntermediateDefFieldFromVelGrid(NiftiImage& velocityFieldGrid, + NiftiImage deformationFields[]); /* *************************************************************** */ void reg_spline_getFlowFieldFromVelocityGrid(nifti_image *velocityFieldGrid, nifti_image *flowField); diff --git a/reg-lib/cuda/CudaAladinContent.cpp b/reg-lib/cuda/CudaAladinContent.cpp index e3bf130e..9e1f94ed 100644 --- a/reg-lib/cuda/CudaAladinContent.cpp +++ b/reg-lib/cuda/CudaAladinContent.cpp @@ -4,8 +4,8 @@ #include /* *************************************************************** */ -CudaAladinContent::CudaAladinContent(nifti_image *referenceIn, - nifti_image *floatingIn, +CudaAladinContent::CudaAladinContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, int *referenceMaskIn, mat44 *transformationMatrixIn, size_t bytesIn, @@ -42,19 +42,17 @@ void CudaAladinContent::InitVars() { mask_d = nullptr; floIJKMat_d = nullptr; - if (reference != nullptr && reference->nbyper != NIFTI_TYPE_FLOAT32) - reg_tools_changeDatatype(reference); - if (floating != nullptr && floating->nbyper != NIFTI_TYPE_FLOAT32) { - reg_tools_changeDatatype(floating); - if (warped != nullptr) - reg_tools_changeDatatype(warped); + if (reference && reference->nbyper != NIFTI_TYPE_FLOAT32) + reference.changeDatatype(NIFTI_TYPE_FLOAT32); + if (floating && floating->nbyper != NIFTI_TYPE_FLOAT32) { + floating.changeDatatype(NIFTI_TYPE_FLOAT32); + if (warped) + warped.changeDatatype(NIFTI_TYPE_FLOAT32); } - - //numBlocks = (blockMatchingParams->activeBlock != nullptr) ? blockMatchingParams->blockNumber[0] * blockMatchingParams->blockNumber[1] * blockMatchingParams->blockNumber[2] : 0; } /* *************************************************************** */ void CudaAladinContent::AllocateCuPtrs() { - if (transformationMatrix != nullptr) { + if (transformationMatrix) { Cuda::Allocate(&transformationMatrix_d, sizeof(mat44) / sizeof(float)); float *tmpMat_h = (float*)malloc(sizeof(mat44)); @@ -63,11 +61,11 @@ void CudaAladinContent::AllocateCuPtrs() { free(tmpMat_h); } - if (referenceMask != nullptr) { + if (referenceMask) { Cuda::Allocate(&mask_d, reference->nvox); Cuda::TransferNiftiToDevice(mask_d, referenceMask, reference->nvox); } - if (reference != nullptr) { + if (reference) { Cuda::Allocate(&referenceImageArray_d, reference->nvox); Cuda::Allocate(&referenceMat_d, sizeof(mat44) / sizeof(float)); @@ -78,15 +76,15 @@ void CudaAladinContent::AllocateCuPtrs() { Cuda::TransferNiftiToDevice(referenceMat_d, targetMat, sizeof(mat44) / sizeof(float)); free(targetMat); } - if (warped != nullptr) { + if (warped) { Cuda::Allocate(&warpedImageArray_d, warped->nvox); Cuda::TransferNiftiToDevice(warpedImageArray_d, warped); } - if (deformationField != nullptr) { + if (deformationField) { Cuda::Allocate(&deformationFieldArray_d, deformationField->nvox); Cuda::TransferNiftiToDevice(deformationFieldArray_d, deformationField); } - if (floating != nullptr) { + if (floating) { Cuda::Allocate(&floatingImageArray_d, floating->nvox); Cuda::Allocate(&floIJKMat_d, sizeof(mat44) / sizeof(float)); @@ -98,28 +96,28 @@ void CudaAladinContent::AllocateCuPtrs() { free(sourceIJKMatrix_h); } - if (blockMatchingParams != nullptr) { - if (blockMatchingParams->referencePosition != nullptr) { + if (blockMatchingParams) { + if (blockMatchingParams->referencePosition) { Cuda::Allocate(&referencePosition_d, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim); Cuda::TransferFromHostToDevice(referencePosition_d, blockMatchingParams->referencePosition, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim); } - if (blockMatchingParams->warpedPosition != nullptr) { + if (blockMatchingParams->warpedPosition) { Cuda::Allocate(&warpedPosition_d, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim); Cuda::TransferFromHostToDevice(warpedPosition_d, blockMatchingParams->warpedPosition, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim); } - if (blockMatchingParams->totalBlock != nullptr) { + if (blockMatchingParams->totalBlock) { Cuda::Allocate(&totalBlock_d, blockMatchingParams->totalBlockNumber); Cuda::TransferNiftiToDevice(totalBlock_d, blockMatchingParams->totalBlock, blockMatchingParams->totalBlockNumber); } } } /* *************************************************************** */ -nifti_image* CudaAladinContent::GetWarped() { +NiftiImage& CudaAladinContent::GetWarped() { DownloadImage(warped, warpedImageArray_d, warped->datatype); return warped; } /* *************************************************************** */ -nifti_image* CudaAladinContent::GetDeformationField() { +NiftiImage& CudaAladinContent::GetDeformationField() { Cuda::TransferFromDeviceToHost((float*)deformationField->data, deformationFieldArray_d, deformationField->nvox); return deformationField; } @@ -131,7 +129,7 @@ _reg_blockMatchingParam* CudaAladinContent::GetBlockMatchingParams() { } /* *************************************************************** */ void CudaAladinContent::SetTransformationMatrix(mat44 *transformationMatrixIn) { - if (transformationMatrix != nullptr) + if (transformationMatrix) Cuda::Free(transformationMatrix_d); AladinContent::SetTransformationMatrix(transformationMatrixIn); @@ -143,28 +141,33 @@ void CudaAladinContent::SetTransformationMatrix(mat44 *transformationMatrixIn) { free(tmpMat_h); } /* *************************************************************** */ -void CudaAladinContent::SetDeformationField(nifti_image *deformationFieldIn) { - if (deformationField != nullptr) +void CudaAladinContent::SetDeformationField(NiftiImage&& deformationFieldIn) { + if (deformationField) Cuda::Free(deformationFieldArray_d); - AladinContent::SetDeformationField(deformationFieldIn); + AladinContent::SetDeformationField(std::move(deformationFieldIn)); Cuda::Allocate(&deformationFieldArray_d, deformationField->nvox); Cuda::TransferNiftiToDevice(deformationFieldArray_d, deformationField); } /* *************************************************************** */ void CudaAladinContent::SetReferenceMask(int *referenceMaskIn) { - if (referenceMask != nullptr) + if (referenceMask) Cuda::Free(mask_d); AladinContent::SetReferenceMask(referenceMaskIn); Cuda::Allocate(&mask_d, reference->nvox); Cuda::TransferNiftiToDevice(mask_d, referenceMaskIn, reference->nvox); } /* *************************************************************** */ -void CudaAladinContent::SetWarped(nifti_image *warped) { - if (warped != nullptr) +void CudaAladinContent::SetWarped(NiftiImage&& warpedIn) { + AladinContent::SetWarped(std::move(warpedIn)); + if (warpedImageArray_d) { Cuda::Free(warpedImageArray_d); - AladinContent::SetWarped(warped); - reg_tools_changeDatatype(warped); + warpedImageArray_d = nullptr; + } + if (!warped) return; + + if (warped->nbyper != NIFTI_TYPE_FLOAT32) + warped.changeDatatype(NIFTI_TYPE_FLOAT32); Cuda::Allocate(&warpedImageArray_d, warped->nvox); Cuda::TransferNiftiToDevice(warpedImageArray_d, warped); @@ -172,19 +175,19 @@ void CudaAladinContent::SetWarped(nifti_image *warped) { /* *************************************************************** */ void CudaAladinContent::SetBlockMatchingParams(_reg_blockMatchingParam* bmp) { AladinContent::SetBlockMatchingParams(bmp); - if (blockMatchingParams->referencePosition != nullptr) { + if (blockMatchingParams->referencePosition) { Cuda::Free(referencePosition_d); //referencePosition Cuda::Allocate(&referencePosition_d, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim); Cuda::TransferFromHostToDevice(referencePosition_d, blockMatchingParams->referencePosition, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim); } - if (blockMatchingParams->warpedPosition != nullptr) { + if (blockMatchingParams->warpedPosition) { Cuda::Free(warpedPosition_d); //warpedPosition Cuda::Allocate(&warpedPosition_d, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim); Cuda::TransferFromHostToDevice(warpedPosition_d, blockMatchingParams->warpedPosition, blockMatchingParams->activeBlockNumber * blockMatchingParams->dim); } - if (blockMatchingParams->totalBlock != nullptr) { + if (blockMatchingParams->totalBlock) { Cuda::Free(totalBlock_d); //activeBlock Cuda::Allocate(&totalBlock_d, blockMatchingParams->totalBlockNumber); @@ -192,51 +195,21 @@ void CudaAladinContent::SetBlockMatchingParams(_reg_blockMatchingParam* bmp) { } } /* *************************************************************** */ -template -void CudaAladinContent::FillImageData(nifti_image *image, float *memoryObject, int datatype) { +void CudaAladinContent::DownloadImage(NiftiImage& image, float *memoryObject, int datatype) { const size_t size = image->nvox; unique_ptr buffer(new float[size]); Cuda::TransferFromDeviceToHost(buffer.get(), memoryObject, size); - free(image->data); - image->datatype = datatype; - image->nbyper = sizeof(DataType); - image->data = malloc(size * image->nbyper); - DataType *data = static_cast(image->data); - for (size_t i = 0; i < size; ++i) - data[i] = static_cast(NiftiImage::clampData(image, buffer[i])); -} -/* *************************************************************** */ -void CudaAladinContent::DownloadImage(nifti_image *image, float *memoryObject, int datatype) { - switch (datatype) { - case NIFTI_TYPE_FLOAT32: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_FLOAT64: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT8: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT8: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT16: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT16: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT32: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT32: - FillImageData(image, memoryObject, datatype); - break; - default: - NR_FATAL_ERROR("Unsupported type"); - } + std::visit([&](auto&& dataType) { + using DataType = std::decay_t; + image->datatype = datatype; + image->nbyper = sizeof(DataType); + image.realloc(); + DataType *data = static_cast(image->data); + for (size_t i = 0; i < size; ++i) + data[i] = static_cast(image.clampData(buffer[i])); + }, image.getDataType()); } /* *************************************************************** */ float* CudaAladinContent::GetReferenceImageArray_d() { @@ -284,33 +257,33 @@ int* CudaAladinContent::GetMask_d() { } /* *************************************************************** */ void CudaAladinContent::FreeCuPtrs() { - if (transformationMatrix_d != nullptr) + if (transformationMatrix_d) Cuda::Free(transformationMatrix_d); - if (referenceImageArray_d != nullptr) + if (referenceImageArray_d) Cuda::Free(referenceImageArray_d); - if (referenceMat_d != nullptr) + if (referenceMat_d) Cuda::Free(referenceMat_d); - if (floatingImageArray_d != nullptr) + if (floatingImageArray_d) Cuda::Free(floatingImageArray_d); - if (floIJKMat_d != nullptr) + if (floIJKMat_d) Cuda::Free(floIJKMat_d); - if (warpedImageArray_d != nullptr) + if (warpedImageArray_d) Cuda::Free(warpedImageArray_d); - if (deformationFieldArray_d != nullptr) + if (deformationFieldArray_d) Cuda::Free(deformationFieldArray_d); - if (mask_d != nullptr) + if (mask_d) Cuda::Free(mask_d); - if (totalBlock_d != nullptr) + if (totalBlock_d) Cuda::Free(totalBlock_d); - if (referencePosition_d != nullptr) + if (referencePosition_d) Cuda::Free(referencePosition_d); - if (warpedPosition_d != nullptr) + if (warpedPosition_d) Cuda::Free(warpedPosition_d); } /* *************************************************************** */ diff --git a/reg-lib/cuda/CudaAladinContent.h b/reg-lib/cuda/CudaAladinContent.h index 6521829a..ac649f21 100644 --- a/reg-lib/cuda/CudaAladinContent.h +++ b/reg-lib/cuda/CudaAladinContent.h @@ -6,8 +6,8 @@ class CudaAladinContent: public AladinContent { public: - CudaAladinContent(nifti_image *referenceIn, - nifti_image *floatingIn, + CudaAladinContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, int *referenceMaskIn = nullptr, mat44 *transformationMatrixIn = nullptr, size_t bytesIn = sizeof(float), @@ -34,8 +34,8 @@ class CudaAladinContent: public AladinContent { // CPU getters with data downloaded from device virtual _reg_blockMatchingParam* GetBlockMatchingParams() override; - virtual nifti_image* GetDeformationField() override; - virtual nifti_image* GetWarped() override; + virtual NiftiImage& GetDeformationField() override; + virtual NiftiImage& GetWarped() override; private: void InitVars(); @@ -54,8 +54,7 @@ class CudaAladinContent: public AladinContent { float *referenceMat_d; float *floIJKMat_d; - template void FillImageData(nifti_image *image, float *memoryObject, int datatype); - void DownloadImage(nifti_image *image, float *memoryObject, int datatype); + void DownloadImage(NiftiImage& image, float *memoryObject, int datatype); #ifdef NR_TESTING public: @@ -64,8 +63,8 @@ class CudaAladinContent: public AladinContent { #endif // Functions for testing virtual void SetTransformationMatrix(mat44 *transformationMatrixIn) override; - virtual void SetWarped(nifti_image *warpedIn) override; - virtual void SetDeformationField(nifti_image *deformationFieldIn) override; + virtual void SetWarped(NiftiImage&& warpedIn) override; + virtual void SetDeformationField(NiftiImage&& deformationFieldIn) override; virtual void SetReferenceMask(int *referenceMaskIn) override; virtual void SetBlockMatchingParams(_reg_blockMatchingParam *bmp) override; }; diff --git a/reg-lib/cuda/CudaAladinContentCreator.h b/reg-lib/cuda/CudaAladinContentCreator.h index 7da8c0fd..4928f685 100644 --- a/reg-lib/cuda/CudaAladinContentCreator.h +++ b/reg-lib/cuda/CudaAladinContentCreator.h @@ -5,8 +5,8 @@ class CudaAladinContentCreator: public AladinContentCreator { 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), diff --git a/reg-lib/cuda/CudaCompute.cu b/reg-lib/cuda/CudaCompute.cu index 08493d4a..db63657c 100644 --- a/reg-lib/cuda/CudaCompute.cu +++ b/reg-lib/cuda/CudaCompute.cu @@ -10,7 +10,7 @@ /* *************************************************************** */ void CudaCompute::ResampleImage(int interpolation, float paddingValue) { CudaContent& con = dynamic_cast(this->con); - const nifti_image *floating = con.Content::GetFloating(); + const NiftiImage& floating = con.Content::GetFloating(); auto resampleImage = floating->nz > 1 ? Cuda::ResampleImage : Cuda::ResampleImage; resampleImage(floating, con.GetFloatingCuda(), @@ -52,7 +52,7 @@ double CudaCompute::CorrectFolding(bool approx) { /* *************************************************************** */ double CudaCompute::ApproxBendingEnergy() { CudaF3dContent& con = dynamic_cast(this->con); - const nifti_image *controlPointGrid = con.F3dContent::GetControlPointGrid(); + const NiftiImage& controlPointGrid = con.F3dContent::GetControlPointGrid(); auto approxBendingEnergy = controlPointGrid->nz > 1 ? Cuda::ApproxBendingEnergy : Cuda::ApproxBendingEnergy; return approxBendingEnergy(controlPointGrid, con.GetControlPointGridCuda()); @@ -60,7 +60,7 @@ double CudaCompute::ApproxBendingEnergy() { /* *************************************************************** */ void CudaCompute::ApproxBendingEnergyGradient(float weight) { CudaF3dContent& con = dynamic_cast(this->con); - nifti_image *controlPointGrid = con.F3dContent::GetControlPointGrid(); + NiftiImage& controlPointGrid = con.F3dContent::GetControlPointGrid(); auto approxBendingEnergyGradient = controlPointGrid->nz > 1 ? Cuda::ApproxBendingEnergyGradient : Cuda::ApproxBendingEnergyGradient; approxBendingEnergyGradient(controlPointGrid, @@ -71,7 +71,7 @@ void CudaCompute::ApproxBendingEnergyGradient(float weight) { /* *************************************************************** */ double CudaCompute::ApproxLinearEnergy() { CudaF3dContent& con = dynamic_cast(this->con); - const nifti_image *controlPointGrid = con.F3dContent::GetControlPointGrid(); + const NiftiImage& controlPointGrid = con.F3dContent::GetControlPointGrid(); auto approxLinearEnergy = controlPointGrid->nz > 1 ? Cuda::ApproxLinearEnergy : Cuda::ApproxLinearEnergy; return approxLinearEnergy(controlPointGrid, con.GetControlPointGridCuda()); @@ -79,7 +79,7 @@ double CudaCompute::ApproxLinearEnergy() { /* *************************************************************** */ void CudaCompute::ApproxLinearEnergyGradient(float weight) { CudaF3dContent& con = dynamic_cast(this->con); - const nifti_image *controlPointGrid = con.F3dContent::GetControlPointGrid(); + const NiftiImage& controlPointGrid = con.F3dContent::GetControlPointGrid(); auto approxLinearEnergyGradient = controlPointGrid->nz > 1 ? Cuda::ApproxLinearEnergyGradient : Cuda::ApproxLinearEnergyGradient; approxLinearEnergyGradient(controlPointGrid, con.GetControlPointGridCuda(), con.GetTransformationGradientCuda(), weight); @@ -145,9 +145,9 @@ void CudaCompute::UpdateControlPointPosition(float *currentDof, const bool optimiseX, const bool optimiseY, const bool optimiseZ) { - const nifti_image *controlPointGrid = dynamic_cast(con).F3dContent::GetControlPointGrid(); + const NiftiImage& controlPointGrid = dynamic_cast(con).F3dContent::GetControlPointGrid(); const bool optZ = optimiseZ && controlPointGrid->nz > 1; - const size_t nVoxels = NiftiImage::calcVoxelNumber(controlPointGrid, 3); + const size_t nVoxels = controlPointGrid.nVoxelsPerVolume(); auto bestDofTexturePtr = Cuda::CreateTextureObject(reinterpret_cast(bestDof), nVoxels, cudaChannelFormatKindFloat, 4); auto gradientTexturePtr = Cuda::CreateTextureObject(reinterpret_cast(gradient), nVoxels, cudaChannelFormatKindFloat, 4); @@ -172,7 +172,7 @@ void CudaCompute::UpdateControlPointPosition(float *currentDof, /* *************************************************************** */ void CudaCompute::GetImageGradient(int interpolation, float paddingValue, int activeTimePoint) { CudaDefContent& con = dynamic_cast(this->con); - const nifti_image *floating = con.Content::GetFloating(); + const NiftiImage& floating = con.Content::GetFloating(); auto getImageGradient = floating->nz > 1 ? Cuda::GetImageGradient : Cuda::GetImageGradient; getImageGradient(floating, con.GetFloatingCuda(), @@ -187,8 +187,8 @@ void CudaCompute::GetImageGradient(int interpolation, float paddingValue, int ac double CudaCompute::GetMaximalLength(bool optimiseX, bool optimiseY, bool optimiseZ) { if (!optimiseX && !optimiseY && !optimiseZ) return 0; CudaF3dContent& con = dynamic_cast(this->con); - nifti_image *transGrad = con.F3dContent::GetTransformationGradient(); - const size_t voxelsPerVolume = NiftiImage::calcVoxelNumber(transGrad, 3); + NiftiImage& transGrad = con.F3dContent::GetTransformationGradient(); + const size_t voxelsPerVolume = transGrad.nVoxelsPerVolume(); if (transGrad->nz <= 1) optimiseZ = false; return Cuda::GetMaximalLength(con.GetTransformationGradientCuda(), voxelsPerVolume, optimiseX, optimiseY, optimiseZ); } @@ -196,8 +196,8 @@ double CudaCompute::GetMaximalLength(bool optimiseX, bool optimiseY, bool optimi void CudaCompute::NormaliseGradient(double maxGradLength, bool optimiseX, bool optimiseY, bool optimiseZ) { if (maxGradLength == 0 || (!optimiseX && !optimiseY && !optimiseZ)) return; CudaF3dContent& con = dynamic_cast(this->con); - nifti_image *transGrad = con.F3dContent::GetTransformationGradient(); - const size_t voxelsPerVolume = NiftiImage::calcVoxelNumber(transGrad, 3); + NiftiImage& transGrad = con.F3dContent::GetTransformationGradient(); + const size_t voxelsPerVolume = transGrad.nVoxelsPerVolume(); if (transGrad->nz <= 1) optimiseZ = false; Cuda::NormaliseGradient(con.GetTransformationGradientCuda(), voxelsPerVolume, maxGradLength, optimiseX, optimiseY, optimiseZ); } @@ -224,8 +224,8 @@ void CudaCompute::GetDefFieldFromVelocityGrid(const bool updateStepNumber) { updateStepNumber); } /* *************************************************************** */ -void CudaCompute::ConvolveImage(const nifti_image *image, float4 *imageCuda) { - const nifti_image *controlPointGrid = dynamic_cast(con).F3dContent::GetControlPointGrid(); +void CudaCompute::ConvolveImage(const NiftiImage& image, float4 *imageCuda) { + const NiftiImage& controlPointGrid = dynamic_cast(con).F3dContent::GetControlPointGrid(); constexpr ConvKernelType kernelType = ConvKernelType::Cubic; float currentNodeSpacing[3]; currentNodeSpacing[0] = currentNodeSpacing[1] = currentNodeSpacing[2] = controlPointGrid->dx; @@ -260,7 +260,7 @@ void CudaCompute::ConvolveImage(const nifti_image *image, float4 *imageCuda) { void CudaCompute::VoxelCentricToNodeCentric(float weight) { CudaF3dContent& con = dynamic_cast(this->con); const mat44 *reorientation = Content::GetIJKMatrix(*con.Content::GetFloating()); - const nifti_image *transGrad = con.F3dContent::GetTransformationGradient(); + const NiftiImage& transGrad = con.F3dContent::GetTransformationGradient(); auto voxelCentricToNodeCentric = transGrad->nz > 1 ? Cuda::VoxelCentricToNodeCentric : Cuda::VoxelCentricToNodeCentric; voxelCentricToNodeCentric(transGrad, @@ -281,21 +281,21 @@ void CudaCompute::ConvolveVoxelBasedMeasureGradient(float weight) { void CudaCompute::ExponentiateGradient(Content& conBwIn) { CudaF3dContent& con = dynamic_cast(this->con); CudaF3dContent& conBw = dynamic_cast(conBwIn); - nifti_image *deformationField = con.Content::GetDeformationField(); - nifti_image *voxelBasedMeasureGradient = con.DefContent::GetVoxelBasedMeasureGradient(); + NiftiImage& deformationField = con.Content::GetDeformationField(); + NiftiImage& voxelBasedMeasureGradient = con.DefContent::GetVoxelBasedMeasureGradient(); float4 *voxelBasedMeasureGradientCuda = con.GetVoxelBasedMeasureGradientCuda(); - nifti_image *controlPointGridBw = conBw.F3dContent::GetControlPointGrid(); + NiftiImage& controlPointGridBw = conBw.F3dContent::GetControlPointGrid(); float4 *controlPointGridBwCuda = conBw.GetControlPointGridCuda(); mat44 *affineTransformationBw = conBw.Content::GetTransformationMatrix(); const int compNum = std::abs(static_cast(controlPointGridBw->intent_p2)); // The number of composition /* Allocate a temporary gradient image to store the backward gradient */ - const size_t voxelGradNumber = NiftiImage::calcVoxelNumber(voxelBasedMeasureGradient, 3); + const size_t voxelGradNumber = voxelBasedMeasureGradient.nVoxelsPerVolume(); NiftiImage warped(voxelBasedMeasureGradient, NiftiImage::Copy::ImageInfo); thrust::device_vector warpedCudaVec(voxelGradNumber); // Create all deformation field images needed for resampling - const size_t defFieldNumber = NiftiImage::calcVoxelNumber(deformationField, 3); + const size_t defFieldNumber = deformationField.nVoxelsPerVolume(); vector defFields(compNum + 1, NiftiImage(deformationField, NiftiImage::Copy::ImageInfo)); vector> defFieldCudaVecs(compNum + 1, thrust::device_vector(defFieldNumber)); @@ -341,8 +341,8 @@ void CudaCompute::UpdateVelocityField(float scale, bool optimiseX, bool optimise if (!optimiseX && !optimiseY && !optimiseZ) return; CudaF3dContent& con = dynamic_cast(this->con); - const nifti_image *controlPointGrid = con.F3dContent::GetControlPointGrid(); - const size_t voxelNumber = NiftiImage::calcVoxelNumber(controlPointGrid, 3); + const NiftiImage& controlPointGrid = con.F3dContent::GetControlPointGrid(); + const size_t voxelNumber = controlPointGrid.nVoxelsPerVolume(); auto scaledGradientCudaPtr = ScaleGradient(con.GetTransformationGradientCuda(), voxelNumber, scale); // Reset the gradient along the axes if appropriate @@ -365,11 +365,11 @@ void CudaCompute::SymmetriseVelocityFields(Content& conBwIn) { CudaF3dContent& con = dynamic_cast(this->con); CudaF3dContent& conBw = dynamic_cast(conBwIn); - nifti_image *controlPointGrid = con.F3dContent::GetControlPointGrid(); - nifti_image *controlPointGridBw = conBw.F3dContent::GetControlPointGrid(); + NiftiImage& controlPointGrid = con.F3dContent::GetControlPointGrid(); + NiftiImage& controlPointGridBw = conBw.F3dContent::GetControlPointGrid(); float4 *controlPointGridCuda = con.GetControlPointGridCuda(); float4 *controlPointGridBwCuda = conBw.GetControlPointGridCuda(); - const size_t voxelNumber = NiftiImage::calcVoxelNumber(controlPointGrid, 3); + const size_t voxelNumber = controlPointGrid.nVoxelsPerVolume(); // In order to ensure symmetry, the forward and backward velocity fields // are averaged in both image spaces: reference and floating @@ -395,9 +395,9 @@ void CudaCompute::SymmetriseVelocityFields(Content& conBwIn) { Cuda::GetDeformationFromDisplacement(controlPointGridBw, controlPointGridBwCuda); } /* *************************************************************** */ -void CudaCompute::DefFieldCompose(const nifti_image *defField) { +void CudaCompute::DefFieldCompose(const NiftiImage& defField) { CudaContent& con = dynamic_cast(this->con); - const size_t voxelNumber = NiftiImage::calcVoxelNumber(defField, 3); + const size_t voxelNumber = defField.nVoxelsPerVolume(); thrust::device_vector defFieldCuda(voxelNumber); Cuda::TransferNiftiToDevice(defFieldCuda.data().get(), defField); auto defFieldCompose = defField->nz > 1 ? Cuda::DefFieldCompose : Cuda::DefFieldCompose; @@ -406,7 +406,7 @@ void CudaCompute::DefFieldCompose(const nifti_image *defField) { /* *************************************************************** */ NiftiImage CudaCompute::ResampleGradient(int interpolation, float padding) { CudaDefContent& con = dynamic_cast(this->con); - const nifti_image *voxelBasedMeasureGradient = con.DefContent::GetVoxelBasedMeasureGradient(); + const NiftiImage& voxelBasedMeasureGradient = con.DefContent::GetVoxelBasedMeasureGradient(); auto resampleGradient = voxelBasedMeasureGradient->nz > 1 ? Cuda::ResampleGradient : Cuda::ResampleGradient; resampleGradient(voxelBasedMeasureGradient, con.GetVoxelBasedMeasureGradientCuda(), @@ -418,7 +418,7 @@ NiftiImage CudaCompute::ResampleGradient(int interpolation, float padding) { con.GetActiveVoxelNumber(), interpolation, padding); - return NiftiImage(con.GetWarpedGradient(), NiftiImage::Copy::Image); + return con.GetWarpedGradient(); } /* *************************************************************** */ void CudaCompute::GetAffineDeformationField(bool compose) { diff --git a/reg-lib/cuda/CudaCompute.h b/reg-lib/cuda/CudaCompute.h index 0982623d..ff342a0c 100644 --- a/reg-lib/cuda/CudaCompute.h +++ b/reg-lib/cuda/CudaCompute.h @@ -35,11 +35,11 @@ class CudaCompute: public Compute { #ifndef NR_TESTING protected: #endif - virtual void DefFieldCompose(const nifti_image *defField) override; + virtual void DefFieldCompose(const NiftiImage& defField) override; virtual NiftiImage ResampleGradient(int interpolation, float padding) override; virtual void VoxelCentricToNodeCentric(float weight) override; private: - void ConvolveImage(const nifti_image*, float4*); + void ConvolveImage(const NiftiImage&, float4*); Cuda::UniquePtr ScaleGradient(const float4*, const size_t, const float); }; diff --git a/reg-lib/cuda/CudaContent.cpp b/reg-lib/cuda/CudaContent.cpp index a1f02b0e..51428168 100644 --- a/reg-lib/cuda/CudaContent.cpp +++ b/reg-lib/cuda/CudaContent.cpp @@ -1,8 +1,8 @@ #include "CudaContent.h" /* *************************************************************** */ -CudaContent::CudaContent(nifti_image *referenceIn, - nifti_image *floatingIn, +CudaContent::CudaContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, int *referenceMaskIn, mat44 *transformationMatrixIn, size_t bytesIn): @@ -24,7 +24,7 @@ CudaContent::~CudaContent() { /* *************************************************************** */ void CudaContent::AllocateReference() { if (reference->nbyper != NIFTI_TYPE_FLOAT32) - reg_tools_changeDatatype(reference); + reference.changeDatatype(NIFTI_TYPE_FLOAT32); Cuda::Allocate(&referenceCuda, reference->nvox); referenceCudaManaged.reset(referenceCuda); Cuda::TransferNiftiToDevice(referenceCuda, reference); @@ -32,7 +32,7 @@ void CudaContent::AllocateReference() { /* *************************************************************** */ void CudaContent::AllocateFloating() { if (floating->nbyper != NIFTI_TYPE_FLOAT32) - reg_tools_changeDatatype(floating); + floating.changeDatatype(NIFTI_TYPE_FLOAT32); Cuda::Allocate(&floatingCuda, floating->nvox); floatingCudaManaged.reset(floatingCuda); Cuda::TransferNiftiToDevice(floatingCuda, floating); @@ -65,13 +65,13 @@ bool CudaContent::IsCurrentComputationDoubleCapable() { return CudaContext::GetInstance().IsCardDoubleCapable(); } /* *************************************************************** */ -nifti_image* CudaContent::GetDeformationField() { +NiftiImage& CudaContent::GetDeformationField() { Cuda::TransferFromDeviceToNifti(deformationField, deformationFieldCuda); return deformationField; } /* *************************************************************** */ -void CudaContent::SetDeformationField(nifti_image *deformationFieldIn) { - Content::SetDeformationField(deformationFieldIn); +void CudaContent::SetDeformationField(NiftiImage&& deformationFieldIn) { + Content::SetDeformationField(std::move(deformationFieldIn)); DeallocateDeformationField(); if (!deformationField) return; @@ -94,7 +94,7 @@ void CudaContent::SetReferenceMask(int *referenceMaskIn) { activeVoxelNumber = 0; if (!referenceMask) return; - const size_t voxelNumber = NiftiImage::calcVoxelNumber(reference, 3); + const size_t voxelNumber = reference.nVoxelsPerVolume(); thrust::host_vector mask(voxelNumber); int *maskPtr = mask.data(); for (size_t i = 0; i < voxelNumber; i++) { @@ -125,17 +125,19 @@ void CudaContent::SetTransformationMatrix(mat44 *transformationMatrixIn) { free(transformationMatrixCptr); } /* *************************************************************** */ -nifti_image* CudaContent::GetWarped() { +NiftiImage& CudaContent::GetWarped() { DownloadImage(warped, warpedCuda, warped->datatype); return warped; } /* *************************************************************** */ -void CudaContent::SetWarped(nifti_image *warpedIn) { - Content::SetWarped(warpedIn); +void CudaContent::SetWarped(NiftiImage&& warpedIn) { + Content::SetWarped(std::move(warpedIn)); DeallocateWarped(); if (!warped) return; - reg_tools_changeDatatype(warped); + if (warped->nbyper != NIFTI_TYPE_FLOAT32) + warped.changeDatatype(NIFTI_TYPE_FLOAT32); + AllocateWarped(); Cuda::TransferNiftiToDevice(warpedCuda, warped); } @@ -144,50 +146,20 @@ void CudaContent::UpdateWarped() { Cuda::TransferNiftiToDevice(warpedCuda, warped); } /* *************************************************************** */ -template -void CudaContent::FillImageData(nifti_image *image, float *memoryObject, int datatype) { +void CudaContent::DownloadImage(NiftiImage& image, float *memoryObject, int datatype) { const size_t size = image->nvox; unique_ptr buffer(new float[size]); Cuda::TransferFromDeviceToHost(buffer.get(), memoryObject, size); - free(image->data); - image->datatype = datatype; - image->nbyper = sizeof(DataType); - image->data = malloc(size * image->nbyper); - DataType *data = static_cast(image->data); - for (size_t i = 0; i < size; ++i) - data[i] = static_cast(NiftiImage::clampData(image, buffer[i])); -} -/* *************************************************************** */ -void CudaContent::DownloadImage(nifti_image *image, float *memoryObject, int datatype) { - switch (datatype) { - case NIFTI_TYPE_FLOAT32: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_FLOAT64: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT8: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT8: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT16: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT16: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_UINT32: - FillImageData(image, memoryObject, datatype); - break; - case NIFTI_TYPE_INT32: - FillImageData(image, memoryObject, datatype); - break; - default: - NR_FATAL_ERROR("Unsupported type"); - } + std::visit([&](auto&& dataType) { + using DataType = std::decay_t; + image->datatype = datatype; + image->nbyper = sizeof(DataType); + image.realloc(); + DataType *data = static_cast(image->data); + for (size_t i = 0; i < size; ++i) + data[i] = static_cast(image.clampData(buffer[i])); + }, image.getDataType()); } /* *************************************************************** */ diff --git a/reg-lib/cuda/CudaContent.h b/reg-lib/cuda/CudaContent.h index d5225ba6..f3deee15 100644 --- a/reg-lib/cuda/CudaContent.h +++ b/reg-lib/cuda/CudaContent.h @@ -6,8 +6,8 @@ class CudaContent: public virtual Content { public: CudaContent() = delete; - CudaContent(nifti_image *referenceIn, - nifti_image *floatingIn, + CudaContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, int *referenceMaskIn = nullptr, mat44 *transformationMatrixIn = nullptr, size_t bytesIn = sizeof(float)); @@ -16,8 +16,8 @@ class CudaContent: public virtual Content { virtual bool IsCurrentComputationDoubleCapable() override; // Getters - virtual nifti_image* GetDeformationField() override; - virtual nifti_image* GetWarped() override; + virtual NiftiImage& GetDeformationField() override; + virtual NiftiImage& GetWarped() override; virtual float* GetReferenceCuda() { return referenceCuda; } virtual float* GetFloatingCuda() { return floatingCuda; } virtual float4* GetDeformationFieldCuda() { return deformationFieldCuda; } @@ -46,8 +46,7 @@ class CudaContent: public virtual Content { void DeallocateDeformationField(); void AllocateWarped(); void DeallocateWarped(); - template void FillImageData(nifti_image *image, float *memoryObject, int datatype); - void DownloadImage(nifti_image *image, float *memoryObject, int datatype); + void DownloadImage(NiftiImage& image, float *memoryObject, int datatype); void SetReferenceCuda(float *referenceCudaIn) { referenceCudaManaged = nullptr; referenceCuda = referenceCudaIn; } void SetFloatingCuda(float *floatingCudaIn) { floatingCudaManaged = nullptr; floatingCuda = floatingCudaIn; } @@ -60,8 +59,8 @@ class CudaContent: public virtual Content { protected: #endif // Functions for testing - virtual void SetDeformationField(nifti_image *deformationFieldIn) override; + virtual void SetDeformationField(NiftiImage&& deformationFieldIn) override; virtual void SetReferenceMask(int *referenceMaskIn) override; virtual void SetTransformationMatrix(mat44 *transformationMatrixIn) override; - virtual void SetWarped(nifti_image *warpedIn) override; + virtual void SetWarped(NiftiImage&& warpedIn) override; }; diff --git a/reg-lib/cuda/CudaContentCreator.h b/reg-lib/cuda/CudaContentCreator.h index 2bd82113..a889c67c 100644 --- a/reg-lib/cuda/CudaContentCreator.h +++ b/reg-lib/cuda/CudaContentCreator.h @@ -5,8 +5,8 @@ class CudaContentCreator: public ContentCreator { public: - virtual Content* Create(nifti_image *reference, - nifti_image *floating, + virtual Content* Create(NiftiImage& reference, + NiftiImage& floating, int *referenceMask = nullptr, mat44 *transformationMatrix = nullptr, size_t bytes = sizeof(float)) override { diff --git a/reg-lib/cuda/CudaDefContent.cpp b/reg-lib/cuda/CudaDefContent.cpp index 72f1c88c..bae8967f 100644 --- a/reg-lib/cuda/CudaDefContent.cpp +++ b/reg-lib/cuda/CudaDefContent.cpp @@ -1,9 +1,9 @@ #include "CudaDefContent.h" /* *************************************************************** */ -CudaDefContent::CudaDefContent(nifti_image *referenceIn, - nifti_image *floatingIn, - nifti_image *localWeightSimIn, +CudaDefContent::CudaDefContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, + NiftiImage *localWeightSimIn, int *referenceMaskIn, mat44 *transformationMatrixIn, size_t bytesIn): @@ -56,12 +56,12 @@ void CudaDefContent::DeallocateVoxelBasedMeasureGradient() { } } /* *************************************************************** */ -nifti_image* CudaDefContent::GetLocalWeightSim() { +NiftiImage& CudaDefContent::GetLocalWeightSim() { Cuda::TransferFromDeviceToNifti(localWeightSim, localWeightSimCuda); return localWeightSim; } /* *************************************************************** */ -nifti_image* CudaDefContent::GetVoxelBasedMeasureGradient() { +NiftiImage& CudaDefContent::GetVoxelBasedMeasureGradient() { Cuda::TransferFromDeviceToNifti(voxelBasedMeasureGradient, voxelBasedMeasureGradientCuda); return voxelBasedMeasureGradient; } @@ -70,7 +70,7 @@ void CudaDefContent::UpdateVoxelBasedMeasureGradient() { Cuda::TransferNiftiToDevice(voxelBasedMeasureGradientCuda, voxelBasedMeasureGradient); } /* *************************************************************** */ -nifti_image* CudaDefContent::GetWarpedGradient() { +NiftiImage& CudaDefContent::GetWarpedGradient() { Cuda::TransferFromDeviceToNifti(warpedGradient, warpedGradientCuda); return warpedGradient; } @@ -80,6 +80,6 @@ void CudaDefContent::UpdateWarpedGradient() { } /* *************************************************************** */ void CudaDefContent::ZeroVoxelBasedMeasureGradient() { - cudaMemset(voxelBasedMeasureGradientCuda, 0, NiftiImage::calcVoxelNumber(voxelBasedMeasureGradient, 3) * sizeof(float4)); + cudaMemset(voxelBasedMeasureGradientCuda, 0, voxelBasedMeasureGradient.nVoxelsPerVolume() * sizeof(float4)); } /* *************************************************************** */ diff --git a/reg-lib/cuda/CudaDefContent.h b/reg-lib/cuda/CudaDefContent.h index 76e09b21..42030543 100644 --- a/reg-lib/cuda/CudaDefContent.h +++ b/reg-lib/cuda/CudaDefContent.h @@ -6,18 +6,18 @@ class CudaDefContent: public virtual DefContent, public virtual CudaContent { public: CudaDefContent() = delete; - CudaDefContent(nifti_image *referenceIn, - nifti_image *floatingIn, - nifti_image *localWeightSimIn = nullptr, + CudaDefContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, + NiftiImage *localWeightSimIn = nullptr, int *referenceMaskIn = nullptr, mat44 *transformationMatrixIn = nullptr, size_t bytesIn = sizeof(float)); virtual ~CudaDefContent(); // Getters - virtual nifti_image* GetLocalWeightSim() override; - virtual nifti_image* GetVoxelBasedMeasureGradient() override; - virtual nifti_image* GetWarpedGradient() override; + virtual NiftiImage& GetLocalWeightSim() override; + virtual NiftiImage& GetVoxelBasedMeasureGradient() override; + virtual NiftiImage& GetWarpedGradient() override; virtual float* GetLocalWeightSimCuda() { return localWeightSimCuda; } virtual float4* GetVoxelBasedMeasureGradientCuda() { return voxelBasedMeasureGradientCuda; } virtual float4* GetWarpedGradientCuda() { return warpedGradientCuda; } diff --git a/reg-lib/cuda/CudaDefContentCreator.h b/reg-lib/cuda/CudaDefContentCreator.h index af3fb561..499d2717 100644 --- a/reg-lib/cuda/CudaDefContentCreator.h +++ b/reg-lib/cuda/CudaDefContentCreator.h @@ -5,9 +5,9 @@ class CudaDefContentCreator: public DefContentCreator { public: - virtual DefContent* Create(nifti_image *reference, - nifti_image *floating, - nifti_image *localWeightSim = nullptr, + virtual DefContent* Create(NiftiImage& reference, + NiftiImage& floating, + NiftiImage *localWeightSim = nullptr, int *referenceMask = nullptr, mat44 *transformationMatrix = nullptr, size_t bytes = sizeof(float)) override { diff --git a/reg-lib/cuda/CudaF3d2ContentCreator.h b/reg-lib/cuda/CudaF3d2ContentCreator.h index 347e07cc..fa6da14a 100644 --- a/reg-lib/cuda/CudaF3d2ContentCreator.h +++ b/reg-lib/cuda/CudaF3d2ContentCreator.h @@ -5,11 +5,11 @@ class CudaF3d2ContentCreator: public F3d2ContentCreator { public: - virtual std::pair Create(nifti_image *reference, - nifti_image *floating, - nifti_image *controlPointGrid, - nifti_image *controlPointGridBw, - nifti_image *localWeightSim = nullptr, + virtual std::pair Create(NiftiImage& reference, + NiftiImage& floating, + NiftiImage& controlPointGrid, + NiftiImage& controlPointGridBw, + NiftiImage *localWeightSim = nullptr, int *referenceMask = nullptr, int *floatingMask = nullptr, mat44 *transformationMatrix = nullptr, diff --git a/reg-lib/cuda/CudaF3dContent.cpp b/reg-lib/cuda/CudaF3dContent.cpp index c6722b9e..c673ca8f 100644 --- a/reg-lib/cuda/CudaF3dContent.cpp +++ b/reg-lib/cuda/CudaF3dContent.cpp @@ -1,10 +1,10 @@ #include "CudaF3dContent.h" /* *************************************************************** */ -CudaF3dContent::CudaF3dContent(nifti_image *referenceIn, - nifti_image *floatingIn, - nifti_image *controlPointGridIn, - nifti_image *localWeightSimIn, +CudaF3dContent::CudaF3dContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, + NiftiImage& controlPointGridIn, + NiftiImage *localWeightSimIn, int *referenceMaskIn, mat44 *transformationMatrixIn, size_t bytesIn): @@ -46,7 +46,7 @@ void CudaF3dContent::DeallocateTransformationGradient() { } } /* *************************************************************** */ -nifti_image* CudaF3dContent::GetControlPointGrid() { +NiftiImage& CudaF3dContent::GetControlPointGrid() { Cuda::TransferFromDeviceToNifti(controlPointGrid, controlPointGridCuda); return controlPointGrid; } @@ -55,7 +55,7 @@ void CudaF3dContent::UpdateControlPointGrid() { Cuda::TransferNiftiToDevice(controlPointGridCuda, controlPointGrid); } /* *************************************************************** */ -nifti_image* CudaF3dContent::GetTransformationGradient() { +NiftiImage& CudaF3dContent::GetTransformationGradient() { Cuda::TransferFromDeviceToNifti(transformationGradient, transformationGradientCuda); return transformationGradient; } @@ -65,6 +65,6 @@ void CudaF3dContent::UpdateTransformationGradient() { } /* *************************************************************** */ void CudaF3dContent::ZeroTransformationGradient() { - cudaMemset(transformationGradientCuda, 0, NiftiImage::calcVoxelNumber(transformationGradient, 3) * sizeof(float4)); + cudaMemset(transformationGradientCuda, 0, transformationGradient.nVoxelsPerVolume() * sizeof(float4)); } /* *************************************************************** */ diff --git a/reg-lib/cuda/CudaF3dContent.h b/reg-lib/cuda/CudaF3dContent.h index ca085945..01d63923 100644 --- a/reg-lib/cuda/CudaF3dContent.h +++ b/reg-lib/cuda/CudaF3dContent.h @@ -6,18 +6,18 @@ class CudaF3dContent: public F3dContent, public CudaDefContent { public: CudaF3dContent() = delete; - CudaF3dContent(nifti_image *referenceIn, - nifti_image *floatingIn, - nifti_image *controlPointGridIn, - nifti_image *localWeightSimIn = nullptr, + CudaF3dContent(NiftiImage& referenceIn, + NiftiImage& floatingIn, + NiftiImage& controlPointGridIn, + NiftiImage *localWeightSimIn = nullptr, int *referenceMaskIn = nullptr, mat44 *transformationMatrixIn = nullptr, size_t bytesIn = sizeof(float)); virtual ~CudaF3dContent(); // Getters - virtual nifti_image* GetControlPointGrid() override; - virtual nifti_image* GetTransformationGradient() override; + virtual NiftiImage& GetControlPointGrid() override; + virtual NiftiImage& GetTransformationGradient() override; virtual float4* GetControlPointGridCuda() { return controlPointGridCuda; } virtual float4* GetTransformationGradientCuda() { return transformationGradientCuda; } diff --git a/reg-lib/cuda/CudaF3dContentCreator.h b/reg-lib/cuda/CudaF3dContentCreator.h index 3e741eb6..af972868 100644 --- a/reg-lib/cuda/CudaF3dContentCreator.h +++ b/reg-lib/cuda/CudaF3dContentCreator.h @@ -5,10 +5,10 @@ class CudaF3dContentCreator: public F3dContentCreator { public: - virtual F3dContent* Create(nifti_image *reference, - nifti_image *floating, - nifti_image *controlPointGrid, - nifti_image *localWeightSim = nullptr, + virtual F3dContent* Create(NiftiImage& reference, + NiftiImage& floating, + NiftiImage& controlPointGrid, + NiftiImage *localWeightSim = nullptr, int *referenceMask = nullptr, mat44 *transformationMatrix = nullptr, size_t bytes = sizeof(float)) override { diff --git a/reg-lib/cuda/CudaMeasureCreator.cpp b/reg-lib/cuda/CudaMeasureCreator.cpp index 3795297d..711b66fb 100644 --- a/reg-lib/cuda/CudaMeasureCreator.cpp +++ b/reg-lib/cuda/CudaMeasureCreator.cpp @@ -47,11 +47,11 @@ void CudaMeasureCreator::Initialise(reg_measure& measure, DefContent& con, DefCo cudaCon.GetLocalWeightSimCuda(), cudaConBw ? cudaConBw->Content::GetReferenceMask() : nullptr, cudaConBw ? cudaConBw->GetReferenceMaskCuda() : nullptr, - cudaConBw ? cudaConBw->Content::GetWarped() : nullptr, + cudaConBw ? static_cast(cudaConBw->Content::GetWarped()) : nullptr, cudaConBw ? cudaConBw->GetWarpedCuda() : nullptr, - cudaConBw ? cudaConBw->DefContent::GetWarpedGradient() : nullptr, + cudaConBw ? static_cast(cudaConBw->DefContent::GetWarpedGradient()) : nullptr, cudaConBw ? cudaConBw->GetWarpedGradientCuda() : nullptr, - cudaConBw ? cudaConBw->DefContent::GetVoxelBasedMeasureGradient() : nullptr, + cudaConBw ? static_cast(cudaConBw->DefContent::GetVoxelBasedMeasureGradient()) : nullptr, cudaConBw ? cudaConBw->GetVoxelBasedMeasureGradientCuda() : nullptr); } /* *************************************************************** */ diff --git a/reg-test/reg_test_affineDeformationField.cpp b/reg-test/reg_test_affineDeformationField.cpp index dc551cf8..055ae1e1 100644 --- a/reg-test/reg_test_affineDeformationField.cpp +++ b/reg-test/reg_test_affineDeformationField.cpp @@ -169,17 +169,14 @@ class AffineDeformationFieldTest { // Set the deformation field if composition is required if (defField) - aladinContent->SetDeformationField(NiftiImage(defField).disown()); + aladinContent->SetDeformationField(NiftiImage(defField)); // Do the calculation for Aladin unique_ptr affineDeformKernel{ platform->CreateKernel(AffineDeformationFieldKernel::GetName(), aladinContent.get()) }; affineDeformKernel->castTo()->Calculate(defField); - // Get the result - NiftiImage resDefField(aladinContent->GetDeformationField(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName + " - Aladin", std::move(resDefField), expRes }); + // Save the results for testing + testCases.push_back({ testName + " - Aladin", std::move(aladinContent->GetDeformationField()), expRes }); // Do the calculation also for Compute using Content // Skip OpenCL as it is not supported @@ -192,17 +189,14 @@ class AffineDeformationFieldTest { // Set the deformation field if composition is required if (defField) - content->SetDeformationField(NiftiImage(defField).disown()); + content->SetDeformationField(NiftiImage(defField)); // Do the calculation unique_ptr compute{ platform->CreateCompute(*content) }; compute->GetAffineDeformationField(defField); - // Get the result - resDefField = NiftiImage(content->GetDeformationField(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName, std::move(resDefField), std::move(expRes) }); + // Save the results for testing + testCases.push_back({ testName, std::move(content->GetDeformationField()), std::move(expRes) }); } } } diff --git a/reg-test/reg_test_blockMatching.cpp b/reg-test/reg_test_blockMatching.cpp index 2243ec2d..9283047c 100644 --- a/reg-test/reg_test_blockMatching.cpp +++ b/reg-test/reg_test_blockMatching.cpp @@ -101,14 +101,14 @@ class BMTest { testData.emplace_back(TestData( "BlockMatching 2D", reference2d, - NiftiImage(contentResampling2d->GetWarped()), + contentResampling2d->GetWarped(), mask2d.get() )); contentResampling2d.release(); testData.emplace_back(TestData( "BlockMatching 3D", reference3d, - NiftiImage(contentResampling3d->GetWarped()), + contentResampling3d->GetWarped(), mask3d.get() )); contentResampling3d.release(); @@ -120,7 +120,6 @@ class BMTest { for (auto&& platformType : PlatformTypes) { // Create images NiftiImage referenceTest(reference); - NiftiImage warpedTest(warped); // Create the contents shared_ptr platform{ new Platform(platformType) }; @@ -137,7 +136,7 @@ class BMTest { 100, 1 ) }; - content->SetWarped(warpedTest.disown()); + content->SetWarped(NiftiImage(warped)); // Initialise the block matching unique_ptr bmKernel{ platform->CreateKernel(BlockMatchingKernel::GetName(), content.get()) }; diff --git a/reg-test/reg_test_composeField.cpp b/reg-test/reg_test_composeField.cpp index affaa42d..ba9395bf 100644 --- a/reg-test/reg_test_composeField.cpp +++ b/reg-test/reg_test_composeField.cpp @@ -137,12 +137,10 @@ class ComposeDeformationFieldTest { unique_ptr content{ contentCreator->Create(reference, reference) }; unique_ptr compute{ platform->CreateCompute(*content) }; // Run the compose - content->SetDeformationField(NiftiImage(outDefField).disown()); + content->SetDeformationField(NiftiImage(outDefField)); compute->DefFieldCompose(defField); - // Get the result - NiftiImage resDefField(content->GetDeformationField(), NiftiImage::Copy::Image); - // Save for testing - testCases.push_back({ testName + " "s + platform->GetName(), std::move(resDefField), expDefField }); + // Save the results for testing + testCases.push_back({ testName + " "s + platform->GetName(), std::move(content->GetDeformationField()), expDefField }); } } } diff --git a/reg-test/reg_test_conjugateGradient.cpp b/reg-test/reg_test_conjugateGradient.cpp index 6f39ef3c..b4c2d212 100644 --- a/reg-test/reg_test_conjugateGradient.cpp +++ b/reg-test/reg_test_conjugateGradient.cpp @@ -242,18 +242,14 @@ TEST_CASE_METHOD(ConjugateGradientTest, "Conjugate Gradient", "[unit]") { // Increase the precision for the output NR_COUT << std::fixed << std::setprecision(10); - // Set the control point grid - NiftiImage img = content->GetControlPointGrid(); - // Use bestControlPointGrid to store bestDof during initialisation of the optimiser - img.copyData(bestControlPointGrid); + // Set the control point grid by using bestControlPointGrid to store bestDof during initialisation of the optimiser + content->F3dContent::GetControlPointGrid().copyData(bestControlPointGrid); content->UpdateControlPointGrid(); // Set the transformation gradients - img = content->GetTransformationGradient(); - img.copyData(transGrad); + content->F3dContent::GetTransformationGradient().copyData(transGrad); content->UpdateTransformationGradient(); - img = contentBw->GetTransformationGradient(); - img.copyData(transGradBw); + contentBw->F3dContent::GetTransformationGradient().copyData(transGradBw); contentBw->UpdateTransformationGradient(); // Create a copy of the control point grid for expected results @@ -266,8 +262,7 @@ TEST_CASE_METHOD(ConjugateGradientTest, "Conjugate Gradient", "[unit]") { UpdateControlPointPosition(controlPointGridExpected, bestControlPointGrid, transGrad, scale, optimiseX, optimiseY, optimiseZ); // Check the results - img = content->GetControlPointGrid(); - const auto cppPtr = img.data(); + const auto cppPtr = content->GetControlPointGrid().data(); const auto cppExpPtr = controlPointGridExpected.data(); for (size_t i = 0; i < controlPointGridExpected.nVoxels(); ++i) { const float cppVal = cppPtr[i]; @@ -306,12 +301,10 @@ TEST_CASE_METHOD(ConjugateGradientTest, "Conjugate Gradient", "[unit]") { gradientBwPtr[i] = distr(gen); } // Update the transformation gradients - img = content->GetTransformationGradient(); - img.copyData(transGrad); + content->F3dContent::GetTransformationGradient().copyData(transGrad); content->UpdateTransformationGradient(); if (isSymmetric) { - img = contentBw->GetTransformationGradient(); - img.copyData(transGradBw); + contentBw->F3dContent::GetTransformationGradient().copyData(transGradBw); contentBw->UpdateTransformationGradient(); } @@ -320,13 +313,11 @@ TEST_CASE_METHOD(ConjugateGradientTest, "Conjugate Gradient", "[unit]") { UpdateGradientValues(transGrad, false, isSymmetric, &transGradBw); // Check the results - img = content->GetTransformationGradient(); - const auto gradPtr = img.data(); + const auto gradPtr = content->GetTransformationGradient().data(); const auto gradExpPtr = transGrad.data(); NiftiImageData gradBwPtr, gradExpBwPtr; if (isSymmetric) { - img = contentBw->GetTransformationGradient(); - gradBwPtr = img.data(); + gradBwPtr = contentBw->GetTransformationGradient().data(); gradExpBwPtr = transGradBw.data(); } for (size_t i = 0; i < transGrad.nVoxels(); ++i) { diff --git a/reg-test/reg_test_getDeformationField.cpp b/reg-test/reg_test_getDeformationField.cpp index 8c6e0c67..5856d89e 100644 --- a/reg-test/reg_test_getDeformationField.cpp +++ b/reg-test/reg_test_getDeformationField.cpp @@ -114,10 +114,8 @@ class GetDeformationFieldTest { unique_ptr compute{ platform->CreateCompute(*content) }; // Compute the deformation field compute->GetDeformationField(false, true); // no composition - use bspline - // Retrieve the deformation field - NiftiImage defField(content->GetDeformationField(), NiftiImage::Copy::Image); - // Save for testing - testCases.push_back({ testName + " "s + platform->GetName(), std::move(defField), std::move(expDefField) }); + // Save the results for testing + testCases.push_back({ testName + " "s + platform->GetName(), std::move(content->GetDeformationField()), std::move(expDefField) }); } } @@ -186,12 +184,10 @@ class GetDeformationFieldTest { unique_ptr content{ contentCreator->Create(reference, reference, controlPointGrid) }; unique_ptr compute{ platform->CreateCompute(*content) }; // Compute the deformation field - content->SetDeformationField(defField.disown()); + content->SetDeformationField(std::move(defField)); compute->GetDeformationField(true, true); // with composition - use bspline - // Retrieve the deformation field - defField = NiftiImage(content->GetDeformationField(), NiftiImage::Copy::Image); - // Save for testing - testCases.push_back({ testName + " "s + platform->GetName(), std::move(defField), std::move(expDefField) }); + // Save the results for testing + testCases.push_back({ testName + " "s + platform->GetName(), std::move(content->GetDeformationField()), std::move(expDefField) }); } } } diff --git a/reg-test/reg_test_imageGradient.cpp b/reg-test/reg_test_imageGradient.cpp index 22e53ad1..310451c4 100644 --- a/reg-test/reg_test_imageGradient.cpp +++ b/reg-test/reg_test_imageGradient.cpp @@ -179,7 +179,7 @@ TEST_CASE("Image Gradient", "[unit]") { NR_COUT << std::fixed << std::setprecision(10); // Set the warped gradient image to host the computation - NiftiImage warpedGradient(content->GetWarpedGradient()); + NiftiImage& warpedGradient = content->GetWarpedGradient(); warpedGradient.setDim(NiftiDim::NDim, defField->ndim); warpedGradient.setDim(NiftiDim::X, 1); warpedGradient.setDim(NiftiDim::Y, 1); @@ -188,14 +188,14 @@ TEST_CASE("Image Gradient", "[unit]") { warpedGradient.recalcVoxelNumber(); // Set the deformation field - content->SetDeformationField(defField.disown()); + content->SetDeformationField(std::move(defField)); // Do the computation unique_ptr compute{ platform->CreateCompute(*content) }; compute->GetImageGradient(interp, 0, 0); // Check all values - warpedGradient = content->GetWarpedGradient(); + content->GetWarpedGradient(); const auto warpedGradPtr = warpedGradient.data(); const size_t nVoxels = warpedGradient.nVoxels(); for (size_t i = 0; i < nVoxels; ++i) { diff --git a/reg-test/reg_test_interpolation.cpp b/reg-test/reg_test_interpolation.cpp index c46e817e..de8df81f 100644 --- a/reg-test/reg_test_interpolation.cpp +++ b/reg-test/reg_test_interpolation.cpp @@ -208,10 +208,10 @@ TEST_CASE("Interpolation", "[unit]") { warped.setDim(NiftiDim::Z, 1); warped.setDim(NiftiDim::U, 1); warped.realloc(); - content->SetWarped(warped.disown()); + content->SetWarped(std::move(warped)); // Set the deformation field - content->SetDeformationField(defField.disown()); + content->SetDeformationField(std::move(defField)); // Do the computation if (isAladinContent) { @@ -223,7 +223,7 @@ TEST_CASE("Interpolation", "[unit]") { } // Check all values - warped = content->GetWarped(); + warped = std::move(content->GetWarped()); const auto warpedPtr = warped.data(); const size_t nVoxels = warped.nVoxels(); for (size_t i = 0; i < nVoxels; ++i) { diff --git a/reg-test/reg_test_lncc.cpp b/reg-test/reg_test_lncc.cpp index e1bcd0ad..f37d0b1c 100644 --- a/reg-test/reg_test_lncc.cpp +++ b/reg-test/reg_test_lncc.cpp @@ -150,7 +150,7 @@ class LnccTest { // Initialise the warped image using the nearest-neighbour interpolation unique_ptr compute{ platform->CreateCompute(*content) }; compute->ResampleImage(0, 0); - content->SetWarped(floating.disown()); + content->SetWarped(NiftiImage(floating)); // Create the measure creator unique_ptr measureCreator{ platform->CreateMeasureCreator() }; // Use LNCC as a measure @@ -159,7 +159,7 @@ class LnccTest { measure_lncc->SetTimePointWeight(0, 1.0); // weight initially set to default value of 1.0 measureCreator->Initialise(*measure_lncc, *content); const double lncc = measure_lncc->GetSimilarityMeasureValue(); - // Save for testing + // Save the results for testing testCases.push_back({ testName, lncc, expLncc }); } } diff --git a/reg-test/reg_test_nmi.cpp b/reg-test/reg_test_nmi.cpp index d3a2770e..c0c7f9f2 100644 --- a/reg-test/reg_test_nmi.cpp +++ b/reg-test/reg_test_nmi.cpp @@ -81,7 +81,7 @@ class NmiTest { // Create the content unique_ptr content{ contentCreator->Create(reference, floating) }; // Initialise the warped image using floating image - content->SetWarped(floating.disown()); + content->SetWarped(NiftiImage(floating)); // Create the measure creator unique_ptr measureCreator{ platform->CreateMeasureCreator() }; // Use NMI as a measure @@ -89,7 +89,7 @@ class NmiTest { measure_nmi->SetTimePointWeight(0, 1.0); // weight initially set to default value of 1.0 measureCreator->Initialise(*measure_nmi, *content); const double nmi = measure_nmi->GetSimilarityMeasureValue(); - + // Save the results for testing testCases.push_back({ testName + " " + platform->GetName(), nmi, expected }); } } diff --git a/reg-test/reg_test_nmi_gradient.cpp b/reg-test/reg_test_nmi_gradient.cpp index 5342e1b1..2ab098b9 100644 --- a/reg-test/reg_test_nmi_gradient.cpp +++ b/reg-test/reg_test_nmi_gradient.cpp @@ -81,7 +81,7 @@ class NmiGradientTest { // Create the content unique_ptr content{ contentCreator->Create(reference, floating) }; // Add some displacements to the deformation field to avoid grid effect - nifti_image *defField = content->Content::GetDeformationField(); + NiftiImage& defField = content->Content::GetDeformationField(); float *defPtr = static_cast(defField->data); for (size_t index = 0; index < defField->nvox; ++index) defPtr[index] += 0.1f; @@ -101,9 +101,9 @@ class NmiGradientTest { // Compute the NMI gradient measure_nmi->GetVoxelBasedSimilarityMeasureGradient(0); // Create an image to store the gradient values - NiftiImage gradientImage(content->GetVoxelBasedMeasureGradient(), NiftiImage::Copy::Image); + NiftiImage gradientImage(content->GetVoxelBasedMeasureGradient()); // Create an image to store the expected gradient values - NiftiImage expectedGradientImage(content->GetDeformationField(), NiftiImage::Copy::Image); + NiftiImage expectedGradientImage(content->GetDeformationField()); // Apply perturbations to each value in the deformation field float *gradPtr = static_cast(expectedGradientImage->data); constexpr float delta = 0.00001f; diff --git a/reg-test/reg_test_normaliseGradient.cpp b/reg-test/reg_test_normaliseGradient.cpp index 08e9b0d3..f77cd68c 100644 --- a/reg-test/reg_test_normaliseGradient.cpp +++ b/reg-test/reg_test_normaliseGradient.cpp @@ -98,8 +98,7 @@ class NormaliseGradientTest { unique_ptr content{ contentCreator->Create(reference, reference, controlPointGrid) }; // Set the transformation gradient image to host the computation - NiftiImage transGrad = content->GetTransformationGradient(); - transGrad.copyData(expTransGrad); + content->F3dContent::GetTransformationGradient().copyData(expTransGrad); content->UpdateTransformationGradient(); // Calculate the maximal length @@ -111,11 +110,8 @@ class NormaliseGradientTest { compute->NormaliseGradient(expMaxLength, optimiseX, optimiseY, optimiseZ); NormaliseGradient(expTransGrad, expMaxLength, optimiseX, optimiseY, optimiseZ); - // Get the results - transGrad = NiftiImage(content->GetTransformationGradient(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName, maxLength, expMaxLength, std::move(transGrad), std::move(expTransGrad) }); + // Save the results for testing + testCases.push_back({ testName, maxLength, expMaxLength, std::move(content->GetTransformationGradient()), std::move(expTransGrad) }); } } } diff --git a/reg-test/reg_test_regr_approxBendingEnergyGradient.cpp b/reg-test/reg_test_regr_approxBendingEnergyGradient.cpp index a2a01bdf..ac30c7c0 100644 --- a/reg-test/reg_test_regr_approxBendingEnergyGradient.cpp +++ b/reg-test/reg_test_regr_approxBendingEnergyGradient.cpp @@ -112,12 +112,9 @@ class ApproxBendingEnergyGradientTest { computeCpu->ApproxBendingEnergyGradient(weight); computeCuda->ApproxBendingEnergyGradient(weight); - // Get the transformation gradients - NiftiImage transGradCpu(contentCpu->GetTransformationGradient(), NiftiImage::Copy::Image); - NiftiImage transGradCuda(contentCuda->GetTransformationGradient(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName, approxBendingEnergyCpu, approxBendingEnergyCuda, std::move(transGradCpu), std::move(transGradCuda) }); + // Save the results for testing + testCases.push_back({ testName, approxBendingEnergyCpu, approxBendingEnergyCuda, + std::move(contentCpu->GetTransformationGradient()), std::move(contentCuda->GetTransformationGradient()) }); } } }; diff --git a/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp b/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp index 530d404b..ab66f84f 100644 --- a/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp +++ b/reg-test/reg_test_regr_approxLinearEnergyGradient.cpp @@ -114,12 +114,9 @@ class ApproxLinearEnergyGradientTest { computeCpu->ApproxLinearEnergyGradient(weight); computeCuda->ApproxLinearEnergyGradient(weight); - // Get the transformation gradients - NiftiImage transGradCpu(contentCpu->GetTransformationGradient(), NiftiImage::Copy::Image); - NiftiImage transGradCuda(contentCuda->GetTransformationGradient(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName, approxLinearEnergyCpu, approxLinearEnergyCuda, std::move(transGradCpu), std::move(transGradCuda) }); + // Save the results for testing + testCases.push_back({ testName, approxLinearEnergyCpu, approxLinearEnergyCuda, + std::move(contentCpu->GetTransformationGradient()), std::move(contentCuda->GetTransformationGradient()) }); } } }; diff --git a/reg-test/reg_test_regr_blockMatching.cpp b/reg-test/reg_test_regr_blockMatching.cpp index 8676f005..076a9945 100644 --- a/reg-test/reg_test_regr_blockMatching.cpp +++ b/reg-test/reg_test_regr_blockMatching.cpp @@ -98,8 +98,8 @@ class BMTest { ) }; // Initialise the warped images - contentCpu->SetWarped(warpedCpu.disown()); - contentCuda->SetWarped(warpedCuda.disown()); + contentCpu->SetWarped(std::move(warpedCpu)); + contentCuda->SetWarped(std::move(warpedCuda)); // Initialise the block matching unique_ptr kernelCpu{ new CpuBlockMatchingKernel(contentCpu.get()) }; diff --git a/reg-test/reg_test_regr_exponentiateGradient.cpp b/reg-test/reg_test_regr_exponentiateGradient.cpp index db24ff79..69e0ee40 100644 --- a/reg-test/reg_test_regr_exponentiateGradient.cpp +++ b/reg-test/reg_test_regr_exponentiateGradient.cpp @@ -115,25 +115,25 @@ class ExponentiateGradientTest { unique_ptr contentBwCuda{ new CudaF3dContent(referenceBwCuda, referenceBwCuda, cppBwCuda) }; // Set the deformation fields - contentCpu->SetDeformationField(defFieldCpu.disown()); - contentCuda->SetDeformationField(defFieldCuda.disown()); + contentCpu->SetDeformationField(std::move(defFieldCpu)); + contentCuda->SetDeformationField(std::move(defFieldCuda)); // Set the voxel-based measure gradient images - NiftiImage voxelGrad = contentCpu->GetVoxelBasedMeasureGradient(); - voxelGrad->sform_code = voxelBasedGrad->sform_code; - voxelGrad->qto_ijk = voxelBasedGrad->qto_ijk; - voxelGrad->qto_xyz = voxelBasedGrad->qto_xyz; - voxelGrad->sto_ijk = voxelBasedGrad->sto_ijk; - voxelGrad->sto_xyz = voxelBasedGrad->sto_xyz; - voxelGrad.copyData(voxelBasedGrad); + NiftiImage& voxelGradCpu = contentCpu->GetVoxelBasedMeasureGradient(); + voxelGradCpu->sform_code = voxelBasedGrad->sform_code; + voxelGradCpu->qto_ijk = voxelBasedGrad->qto_ijk; + voxelGradCpu->qto_xyz = voxelBasedGrad->qto_xyz; + voxelGradCpu->sto_ijk = voxelBasedGrad->sto_ijk; + voxelGradCpu->sto_xyz = voxelBasedGrad->sto_xyz; + voxelGradCpu.copyData(voxelBasedGrad); contentCpu->UpdateVoxelBasedMeasureGradient(); - voxelGrad = contentCuda->DefContent::GetVoxelBasedMeasureGradient(); - voxelGrad->sform_code = voxelBasedGrad->sform_code; - voxelGrad->qto_ijk = voxelBasedGrad->qto_ijk; - voxelGrad->qto_xyz = voxelBasedGrad->qto_xyz; - voxelGrad->sto_ijk = voxelBasedGrad->sto_ijk; - voxelGrad->sto_xyz = voxelBasedGrad->sto_xyz; - voxelGrad.copyData(voxelBasedGrad); + NiftiImage& voxelGradCuda = contentCuda->DefContent::GetVoxelBasedMeasureGradient(); + voxelGradCuda->sform_code = voxelBasedGrad->sform_code; + voxelGradCuda->qto_ijk = voxelBasedGrad->qto_ijk; + voxelGradCuda->qto_xyz = voxelBasedGrad->qto_xyz; + voxelGradCuda->sto_ijk = voxelBasedGrad->sto_ijk; + voxelGradCuda->sto_xyz = voxelBasedGrad->sto_xyz; + voxelGradCuda.copyData(voxelBasedGrad); contentCuda->UpdateVoxelBasedMeasureGradient(); // Create the computes @@ -144,12 +144,8 @@ class ExponentiateGradientTest { computeCpu->ExponentiateGradient(*contentBwCpu); computeCuda->ExponentiateGradient(*contentBwCuda); - // Get the results - NiftiImage voxelGradCpu(contentCpu->GetVoxelBasedMeasureGradient(), NiftiImage::Copy::Image); - NiftiImage voxelGradCuda(contentCuda->GetVoxelBasedMeasureGradient(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName, std::move(voxelGradCpu), std::move(voxelGradCuda) }); + // Save the results for testing + testCases.push_back({ testName, std::move(contentCpu->GetVoxelBasedMeasureGradient()), std::move(contentCuda->GetVoxelBasedMeasureGradient()) }); } } }; diff --git a/reg-test/reg_test_regr_getDeformationField.cpp b/reg-test/reg_test_regr_getDeformationField.cpp index a4e8cc11..14daeb93 100644 --- a/reg-test/reg_test_regr_getDeformationField.cpp +++ b/reg-test/reg_test_regr_getDeformationField.cpp @@ -72,14 +72,13 @@ class GetDeformationFieldTest { testName += " "s + platform->GetName() + " Composition="s + std::to_string(composition) + " Bspline="s + std::to_string(bspline); unique_ptr content{ contentCreator->Create(reference, reference, controlPointGrid) }; unique_ptr compute{ platform->CreateCompute(*content) }; - NiftiImage expDefField(content->Content::GetDeformationField(), NiftiImage::Copy::Image); + NiftiImage expDefField(content->Content::GetDeformationField()); // Compute the deformation field compute->GetDeformationField(composition, bspline); - NiftiImage defField(content->GetDeformationField(), NiftiImage::Copy::Image); // Compute the expected deformation field GetDeformationField(controlPointGrid, expDefField, content->GetReferenceMask(), composition, bspline); - // Save for testing - testCases.push_back({ std::move(testName), std::move(defField), std::move(expDefField) }); + // Save the results for testing + testCases.push_back({ std::move(testName), std::move(content->GetDeformationField()), std::move(expDefField) }); } } } diff --git a/reg-test/reg_test_regr_kernelConvolution.cpp b/reg-test/reg_test_regr_kernelConvolution.cpp index a65e4879..6d5f6a63 100644 --- a/reg-test/reg_test_regr_kernelConvolution.cpp +++ b/reg-test/reg_test_regr_kernelConvolution.cpp @@ -121,8 +121,8 @@ class KernelConvolutionTest { ) }; // Use deformation fields to store images - contentCpu->SetDeformationField(imageCpu.disown()); - contentCuda->SetDeformationField(imageCuda.disown()); + contentCpu->SetDeformationField(std::move(imageCpu)); + contentCuda->SetDeformationField(std::move(imageCuda)); // Create the kernel convolution function for CUDA auto cudaKernelConvolution = Cuda::KernelConvolution; @@ -136,12 +136,8 @@ class KernelConvolutionTest { reg_tools_kernelConvolution(contentCpu->GetDeformationField(), sigmaValues.data(), ConvKernelType(kernelType), nullptr, activeTimePoints, activeAxes); cudaKernelConvolution(contentCuda->Content::GetDeformationField(), contentCuda->GetDeformationFieldCuda(), sigmaValues.data(), activeTimePoints, activeAxes); - // Get the images - imageCpu = NiftiImage(contentCpu->GetDeformationField(), NiftiImage::Copy::Image); - imageCuda = NiftiImage(contentCuda->GetDeformationField(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName, std::move(imageCpu), std::move(imageCuda) }); + // Save the results for testing + testCases.push_back({ testName, std::move(contentCpu->GetDeformationField()), std::move(contentCuda->GetDeformationField()) }); } } }; diff --git a/reg-test/reg_test_regr_lts.cpp b/reg-test/reg_test_regr_lts.cpp index 681a8ffc..d215df69 100644 --- a/reg-test/reg_test_regr_lts.cpp +++ b/reg-test/reg_test_regr_lts.cpp @@ -110,8 +110,8 @@ class LtsTest { ) }; // Initialise the warped images - contentCpu->SetWarped(warpedCpu.disown()); - contentCuda->SetWarped(warpedCuda.disown()); + contentCpu->SetWarped(std::move(warpedCpu)); + contentCuda->SetWarped(std::move(warpedCuda)); // Initialise the block matching and run it on the CPU unique_ptr bmKernelCpu{ new CpuBlockMatchingKernel(contentCpu.get()) }; diff --git a/reg-test/reg_test_regr_measure.cpp b/reg-test/reg_test_regr_measure.cpp index 08b25515..3deba91b 100644 --- a/reg-test/reg_test_regr_measure.cpp +++ b/reg-test/reg_test_regr_measure.cpp @@ -113,8 +113,8 @@ class MeasureTest { NiftiImage localWeightSimCpu(localWeightSim), localWeightSimCuda(localWeightSim); // Create the contents - auto contentsCpu = contentCreatorCpu->Create(referenceCpu, floatingCpu, controlPointGridCpu, controlPointGridCpuBw, localWeightSimCpu, nullptr, nullptr, nullptr, nullptr, sizeof(float)); - auto contentsCuda = contentCreatorCuda->Create(referenceCuda, floatingCuda, controlPointGridCuda, controlPointGridCudaBw, localWeightSimCuda, nullptr, nullptr, nullptr, nullptr, sizeof(float)); + auto contentsCpu = contentCreatorCpu->Create(referenceCpu, floatingCpu, controlPointGridCpu, controlPointGridCpuBw, &localWeightSimCpu, nullptr, nullptr, nullptr, nullptr, sizeof(float)); + auto contentsCuda = contentCreatorCuda->Create(referenceCuda, floatingCuda, controlPointGridCuda, controlPointGridCudaBw, &localWeightSimCuda, nullptr, nullptr, nullptr, nullptr, sizeof(float)); if (!isSymmetric) { delete contentsCpu.second; delete contentsCuda.second; @@ -184,12 +184,9 @@ class MeasureTest { measureCuda->GetVoxelBasedSimilarityMeasureGradient(t); } - // Get the voxel-based similarity measure gradients - NiftiImage voxelBasedGradCpu(contentCpu->GetVoxelBasedMeasureGradient(), NiftiImage::Copy::Image); - NiftiImage voxelBasedGradCuda(contentCuda->GetVoxelBasedMeasureGradient(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName, simMeasureCpu, simMeasureCuda, std::move(voxelBasedGradCpu), std::move(voxelBasedGradCuda) }); + // Save the results for testing + testCases.push_back({ testName, simMeasureCpu, simMeasureCuda, std::move(contentCpu->GetVoxelBasedMeasureGradient()), + std::move(contentCuda->GetVoxelBasedMeasureGradient()) }); } } }; diff --git a/reg-test/reg_test_regr_resampleGradient.cpp b/reg-test/reg_test_regr_resampleGradient.cpp index 0eadbce3..5a24573b 100644 --- a/reg-test/reg_test_regr_resampleGradient.cpp +++ b/reg-test/reg_test_regr_resampleGradient.cpp @@ -92,25 +92,25 @@ class ResampleGradientTest { unique_ptr contentCuda{ new CudaDefContent(referenceCuda, referenceCuda) }; // Set the deformation fields - contentCpu->SetDeformationField(defFieldCpu.disown()); - contentCuda->SetDeformationField(defFieldCuda.disown()); + contentCpu->SetDeformationField(std::move(defFieldCpu)); + contentCuda->SetDeformationField(std::move(defFieldCuda)); // Set the voxel-based measure gradient images - NiftiImage voxelGrad = contentCpu->GetVoxelBasedMeasureGradient(); - voxelGrad->sform_code = voxelBasedGrad->sform_code; - voxelGrad->qto_ijk = voxelBasedGrad->qto_ijk; - voxelGrad->qto_xyz = voxelBasedGrad->qto_xyz; - voxelGrad->sto_ijk = voxelBasedGrad->sto_ijk; - voxelGrad->sto_xyz = voxelBasedGrad->sto_xyz; - voxelGrad.copyData(voxelBasedGrad); + NiftiImage& voxelGradCpu = contentCpu->GetVoxelBasedMeasureGradient(); + voxelGradCpu->sform_code = voxelBasedGrad->sform_code; + voxelGradCpu->qto_ijk = voxelBasedGrad->qto_ijk; + voxelGradCpu->qto_xyz = voxelBasedGrad->qto_xyz; + voxelGradCpu->sto_ijk = voxelBasedGrad->sto_ijk; + voxelGradCpu->sto_xyz = voxelBasedGrad->sto_xyz; + voxelGradCpu.copyData(voxelBasedGrad); contentCpu->UpdateVoxelBasedMeasureGradient(); - voxelGrad = contentCuda->DefContent::GetVoxelBasedMeasureGradient(); - voxelGrad->sform_code = voxelBasedGrad->sform_code; - voxelGrad->qto_ijk = voxelBasedGrad->qto_ijk; - voxelGrad->qto_xyz = voxelBasedGrad->qto_xyz; - voxelGrad->sto_ijk = voxelBasedGrad->sto_ijk; - voxelGrad->sto_xyz = voxelBasedGrad->sto_xyz; - voxelGrad.copyData(voxelBasedGrad); + NiftiImage& voxelGradCuda = contentCuda->DefContent::GetVoxelBasedMeasureGradient(); + voxelGradCuda->sform_code = voxelBasedGrad->sform_code; + voxelGradCuda->qto_ijk = voxelBasedGrad->qto_ijk; + voxelGradCuda->qto_xyz = voxelBasedGrad->qto_xyz; + voxelGradCuda->sto_ijk = voxelBasedGrad->sto_ijk; + voxelGradCuda->sto_xyz = voxelBasedGrad->sto_xyz; + voxelGradCuda.copyData(voxelBasedGrad); contentCuda->UpdateVoxelBasedMeasureGradient(); // Create the computes @@ -121,7 +121,7 @@ class ResampleGradientTest { NiftiImage warpedCpu = computeCpu->ResampleGradient(1, -2.f); NiftiImage warpedCuda = computeCuda->ResampleGradient(1, -2.f); - // Save for testing + // Save the results for testing testCases.push_back({ testName, std::move(warpedCpu), std::move(warpedCuda) }); } } diff --git a/reg-test/reg_test_regr_symmetriseVelocityFields.cpp b/reg-test/reg_test_regr_symmetriseVelocityFields.cpp index d7149814..0abbd3b4 100644 --- a/reg-test/reg_test_regr_symmetriseVelocityFields.cpp +++ b/reg-test/reg_test_regr_symmetriseVelocityFields.cpp @@ -103,11 +103,11 @@ class SymmetriseVelocityFieldsTest { computeCpu->SymmetriseVelocityFields(*contentBwCpu); computeCuda->SymmetriseVelocityFields(*contentBwCuda); - // Get the results of CUDA since CPU results are already inplace - contentCuda->GetControlPointGrid(); - contentBwCuda->GetControlPointGrid(); + // Since CPU results are already inplace, get CUDA results by destructing their contents + contentCuda = nullptr; + contentBwCuda = nullptr; - // Save for testing + // Save the results for testing testCases.push_back({ testName, std::move(cppCpu), std::move(cppBwCpu), std::move(cppCuda), std::move(cppBwCuda) }); } } diff --git a/reg-test/reg_test_regr_updateVelocityField.cpp b/reg-test/reg_test_regr_updateVelocityField.cpp index 5e85062b..4d40e18a 100644 --- a/reg-test/reg_test_regr_updateVelocityField.cpp +++ b/reg-test/reg_test_regr_updateVelocityField.cpp @@ -78,11 +78,9 @@ class UpdateVelocityFieldTest { unique_ptr contentCuda{ new CudaF3dContent(referenceCuda, referenceCuda, cppCuda) }; // Set the transformation gradient image to host the computation - NiftiImage transGradCpu = contentCpu->GetTransformationGradient(); - transGradCpu.copyData(transGrad); + contentCpu->GetTransformationGradient().copyData(transGrad); contentCpu->UpdateTransformationGradient(); - NiftiImage transGradCuda = contentCuda->GetTransformationGradient(); - transGradCuda.copyData(transGrad); + contentCuda->F3dContent::GetTransformationGradient().copyData(transGrad); contentCuda->UpdateTransformationGradient(); // Create the computes @@ -93,12 +91,9 @@ class UpdateVelocityFieldTest { computeCpu->UpdateVelocityField(scale, optimiseX, optimiseY, optimiseZ); computeCuda->UpdateVelocityField(scale, optimiseX, optimiseY, optimiseZ); - // Get the results - transGradCpu = NiftiImage(contentCpu->GetTransformationGradient(), NiftiImage::Copy::Image); - transGradCuda = NiftiImage(contentCuda->GetTransformationGradient(), NiftiImage::Copy::Image); - - // Save for testing - testCases.push_back({ testName, std::move(transGradCpu), std::move(transGradCuda) }); + // Save the results for testing + testCases.push_back({ testName, std::move(contentCpu->GetTransformationGradient()), + std::move(contentCuda->GetTransformationGradient()) }); } } } diff --git a/reg-test/reg_test_voxelCentricToNodeCentric.cpp b/reg-test/reg_test_voxelCentricToNodeCentric.cpp index 7d807217..b588989a 100644 --- a/reg-test/reg_test_voxelCentricToNodeCentric.cpp +++ b/reg-test/reg_test_voxelCentricToNodeCentric.cpp @@ -90,11 +90,11 @@ class VoxelCentricToNodeCentricTest { unique_ptr content{ contentCreator->Create(reference, reference, controlPointGrid) }; // Set the matrices required for computation - nifti_image *floating = content->Content::GetFloating(); + NiftiImage& floating = content->Content::GetFloating(); if (floating->sform_code > 0) floating->sto_ijk = matrices[0]; else floating->qto_ijk = matrices[0]; - NiftiImage transGrad = content->F3dContent::GetTransformationGradient(); + NiftiImage& transGrad = content->F3dContent::GetTransformationGradient(); static int sfc = 0; transGrad->sform_code = sfc++ % 2; if (transGrad->sform_code > 0) @@ -104,7 +104,7 @@ class VoxelCentricToNodeCentricTest { nifti_add_extension(transGrad, reinterpret_cast(&invMatrix), sizeof(invMatrix), NIFTI_ECODE_IGNORE); // Set the voxel-based measure gradient to host the computation - NiftiImage voxelGrad = content->F3dContent::GetVoxelBasedMeasureGradient(); + NiftiImage& voxelGrad = content->F3dContent::GetVoxelBasedMeasureGradient(); if (voxelGrad->sform_code > 0) voxelGrad->sto_ijk = matrices[3]; else voxelGrad->qto_ijk = matrices[3]; @@ -119,25 +119,26 @@ class VoxelCentricToNodeCentricTest { // Extract the node-based NMI gradient from the voxel-based NMI gradient unique_ptr compute{ platform->CreateCompute(*content) }; compute->VoxelCentricToNodeCentric(weight); - transGrad = NiftiImage(content->GetTransformationGradient(), NiftiImage::Copy::Image); - testCases.push_back({ testName + " "s + platform->GetName() + " Weight="s + std::to_string(weight), std::move(transGrad), std::move(expTransGrad) }); + // Save the results for testing + testCases.push_back({ testName + " "s + platform->GetName() + " Weight="s + std::to_string(weight), + std::move(content->GetTransformationGradient()), std::move(expTransGrad) }); } } } template - void VoxelCentricToNodeCentric(const nifti_image *floating, NiftiImage& nodeGrad, const NiftiImage& voxelGrad, float weight) { + void VoxelCentricToNodeCentric(const NiftiImage& floating, NiftiImage& nodeGrad, const NiftiImage& voxelGrad, float weight) { const mat44 *voxelToMillimetre = floating->sform_code > 0 ? &floating->sto_ijk : &floating->qto_ijk; const bool is3d = nodeGrad->nz > 1; - const size_t nodeNumber = NiftiImage::calcVoxelNumber(nodeGrad, 3); + const size_t nodeNumber = nodeGrad.nVoxelsPerVolume(); auto nodePtr = nodeGrad.data(); auto nodePtrX = nodePtr.begin(); auto nodePtrY = nodePtrX + nodeNumber; auto nodePtrZ = nodePtrY + nodeNumber; - const size_t voxelNumber = NiftiImage::calcVoxelNumber(voxelGrad, 3); + const size_t voxelNumber = voxelGrad.nVoxelsPerVolume(); auto voxelPtr = voxelGrad.data(); auto voxelPtrX = voxelPtr.begin(); auto voxelPtrY = voxelPtrX + voxelNumber;