B spline2

src/Registration/NiftyMoMo/BSplineTransformation.cpp

@@ -1337,7 +1337,7 @@ double BSplineTransformation::GetConstraintValue()
 //------------------------------------------------------------------
 // BSplineTransformation::GetDVFGradientWRTTransformationParameters
 //------------------------------------------------------------------
-BSplineTransformation::PrecisionType* BSplineTransformation::GetDVFGradientWRTTransformationParameters( nifti_image* denseDVFIn, nifti_image* sourceImage )
+BSplineTransformation::PrecisionType* BSplineTransformation::GetDVFGradientWRTTransformationParameters( nifti_image* denseDVFIn )
 {
 #ifndef NDEBUG
   std::cout << "Called BSplineTransformation::GetDVFGradientWRTTransformationParameters()" << std::endl;
@@ -1369,12 +1369,6 @@ BSplineTransformation::PrecisionType* BSplineTransformation::GetDVFGradientWRTTr
     nmm_exit( 1, __FILE__, __LINE__ );
   }
 
-  // Need to reorientate from index gradient to real-world gradient
-  // Note: Performing the reorientation here is way more efficient, since
-  //       only the transformation parameters need to be touched (and not)
-  //       the complete DVF
-  this->ReorientateVectorImage( outDVFGradWRTTrafoParams, sourceImage->sto_ijk );
-
   // Copy over the data pointer from the the image and detach it. Then delete the image.
   PrecisionType* outDVFGradWRTTrafoParamData = (PrecisionType*) outDVFGradWRTTrafoParams->data;
   outDVFGradWRTTrafoParams->data = nullptr;

src/Registration/NiftyMoMo/Transformation.cpp

@@ -134,6 +134,12 @@ void Transformation::GetImageGradientWRTDVF( nifti_image* sourceImage, nifti_ima
                         this->interpolation, 
                         this->warpedPaddingValue, 0 );
 
+  // Need to reorientate from index gradient to real-world gradient
+  // Note: Performing the reorientation here is way more efficient, since
+  //       only the transformation parameters need to be touched (and not)
+  //       the complete DVF
+  this->ReorientateVectorImage( outWarpedGradientImage, sourceImage->sto_ijk );
+
   return;
 }
 

src/Registration/NiftyMoMo/include/sirf/NiftyMoMo/BSplineTransformation.h

@@ -87,7 +87,7 @@ class BSplineTransformation :public Transformation
    *  \param denseDVFIn Deformation vector field 
    *  \param sourceImage The image for which the gradient is calculated
    */
-  virtual BSplineTransformation::PrecisionType* GetDVFGradientWRTTransformationParameters( nifti_image* denseDVFIn, nifti_image* sourceImage );
+  virtual BSplineTransformation::PrecisionType* GetDVFGradientWRTTransformationParameters( nifti_image* denseDVFIn );
 
   /** Calculate the gradient of the constraint term (regularisation) for the transformation with the current parameters.
    */

src/Registration/NiftyMoMo/include/sirf/NiftyMoMo/Transformation.h

@@ -87,7 +87,7 @@ class Transformation
    *  \return pointer to the internally allocated output data of length numberOfParameters. Needs to be freed outside 
    *          of this class since no internal references will be kept.
    */
-  virtual PrecisionType* GetDVFGradientWRTTransformationParameters( nifti_image* denseDVFIn, nifti_image* sourceImage ) = 0;
+  virtual PrecisionType* GetDVFGradientWRTTransformationParameters( nifti_image* denseDVFIn ) = 0;
 
   /** Get the gradient of the regularisation/constraint with respect to the transformation parameters
    *  \f$ \frac{\partial \mathcal{R}_t}{\partial\mathbf{M}_t} \f$

src/Registration/cReg/CMakeLists.txt

@@ -32,7 +32,10 @@ SET(SOURCES
   "NiftiImageData3D.cpp"
   "NiftiImageData3DTensor.cpp"
   "NiftiImageData3DDeformation.cpp"
-  "NiftiImageData3DDisplacement.cpp")
+  "NiftiImageData3DDisplacement.cpp"
+  "NiftiImageData3DBSpline.cpp"
+  "ControlPointGridToDeformationConverter.cpp"
+  )
 
 # If we're also wrapping to python or matlab, include the c-files
 IF(BUILD_PYTHON OR BUILD_MATLAB)

src/Registration/cReg/ControlPointGridToDeformationConverter.cpp

@@ -0,0 +1,111 @@
+/*
+SyneRBI Synergistic Image Reconstruction Framework (SIRF)
+Copyright 2020 University College London
+
+This is software developed for the Collaborative Computational
+Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
+(http://www.ccpsynerbi.ac.uk/).
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+http://www.apache.org/licenses/LICENSE-2.0
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+*/
+
+/*!
+\file
+\ingroup Registration
+\brief Class for converting control point grids to deformation field transformations.
+
+\author Richard Brown
+\author SyneRBI
+*/
+
+#include "sirf/Reg/ControlPointGridToDeformationConverter.h"
+#include "sirf/Reg/NiftiImageData3DDeformation.h"
+#include "sirf/Reg/NiftiImageData3DBSpline.h"
+#include "sirf/NiftyMoMo/BSplineTransformation.h"
+
+using namespace sirf;
+
+template<class dataType>
+ControlPointGridToDeformationConverter<dataType>::
+ControlPointGridToDeformationConverter()
+{
+    for (unsigned i=0; i<3; ++i)
+        _spacing[i] = std::numeric_limits<float>::quiet_NaN();
+}
+
+template<class dataType>
+void
+ControlPointGridToDeformationConverter<dataType>::
+set_cpg_spacing(const float spacing[3])
+{
+    for (unsigned i=0; i<3; ++i)
+        _spacing[i] = spacing[i];
+}
+
+template<class dataType>
+void
+ControlPointGridToDeformationConverter<dataType>::
+set_reference_image(const NiftiImageData<dataType> &ref)
+{
+    _template_ref_sptr = ref.clone();
+}
+
+template<class dataType>
+NiftiImageData3DDeformation<dataType>
+ControlPointGridToDeformationConverter<dataType>::
+forward(const NiftiImageData3DBSpline<dataType> &cpg) const
+{
+    check_is_set_up();
+//    NiftiImageData3DDeformation<float> dvf;
+//    dvf.create_from_cpp(cpg, *_template_ref_sptr);
+//    return dvf;
+    return cpg.get_as_deformation_field(*_template_ref_sptr);
+}
+
+template<class dataType>
+NiftiImageData3DBSpline<dataType>
+ControlPointGridToDeformationConverter<dataType>::
+backward(const NiftiImageData3DDeformation<dataType> &dvf) const
+{
+    check_is_set_up();
+    // not marked const, so copy
+    float spacing_nonconst[3] = {_spacing[0], _spacing[1], _spacing[2]};
+    // Get raw nifti_image from reference image
+    nifti_image *ref_ptr = _template_ref_sptr->get_raw_nifti_sptr().get();
+    // Create the NiftyMoMo bspline transformation class
+    NiftyMoMo::BSplineTransformation bspline(ref_ptr, 1, spacing_nonconst);
+    // Get cpg_ptr
+    nifti_image *cpg_ptr = bspline.GetTransformationAsImage();
+    // Convert DVF to CPG
+    std::shared_ptr<NiftiImageData3DDeformation<dataType> > dvf_sptr = dvf.clone();
+    cpg_ptr->data = bspline.GetDVFGradientWRTTransformationParameters(dvf_sptr->get_raw_nifti_sptr().get());
+    cpg_ptr->intent_p1 = SPLINE_VEL_GRID;
+    return NiftiImageData3DBSpline<dataType>(*cpg_ptr);
+}
+
+template<class dataType>
+void ControlPointGridToDeformationConverter<dataType>::
+check_is_set_up() const
+{
+    // Has spacing been set?
+    for (unsigned i=0; i<3; ++i)
+        if (std::isnan(_spacing[i]))
+            throw std::runtime_error("ControlPointGridToDeformationConverter: Set CPG spacing.");
+
+    // Has template deformation been set?
+    if (!_template_ref_sptr)
+        throw std::runtime_error("ControlPointGridToDeformationConverter: Set template DVF.");
+}
+
+namespace sirf {
+template class ControlPointGridToDeformationConverter<float>;
+}

src/Registration/cReg/NiftiImageData.cpp

@@ -36,6 +36,7 @@ limitations under the License.
 #include "sirf/Reg/NiftiImageData3DTensor.h"
 #include "sirf/Reg/NiftiImageData3DDeformation.h"
 #include "sirf/Reg/NiftiImageData3DDisplacement.h"
+#include "sirf/Reg/NiftiImageData3DBSpline.h"
 #include "sirf/Reg/AffineTransformation.h"
 #include "sirf/Reg/NiftyResample.h"
 #include <iomanip>
@@ -549,7 +550,10 @@ void NiftiImageData<dataType>::check_dimensions(const NiftiImageDataType image_t
     if        (image_type == _3D)       { ndim= 3; nt= 1; nu= 1; intent_code = NIFTI_INTENT_NONE;   intent_p1=-1;         }
     else   if (image_type == _3DTensor) { ndim= 5; nt= 1; nu= 3; intent_code = NIFTI_INTENT_VECTOR; intent_p1=-1;         }
     else   if (image_type == _3DDisp)   { ndim= 5; nt= 1; nu= 3; intent_code = NIFTI_INTENT_VECTOR; intent_p1=DISP_FIELD; }
-    else /*if (image_type == _3DDef)*/  { ndim= 5; nt= 1; nu= 3; intent_code = NIFTI_INTENT_VECTOR; intent_p1=DEF_FIELD;  }
+    else   if (image_type == _3DDef)    { ndim= 5; nt= 1; nu= 3; intent_code = NIFTI_INTENT_VECTOR; intent_p1=DEF_FIELD;  }
+    else   if (image_type == _3DBSpl)   { ndim= 5; nt= 1; nu= 3; intent_code = NIFTI_INTENT_VECTOR; intent_p1=SPLINE_VEL_GRID;  }
+    else
+        throw std::runtime_error("NiftiImageData::check_dimensions: Unknown image type");
 
     // Check everthing is as it should be. -1 means we don't care about it
     // (e.g., NiftiImageData3D doesn't care about intent_p1, which is used by NiftyReg for Disp/Def fields)
@@ -571,16 +575,19 @@ void NiftiImageData<dataType>::check_dimensions(const NiftiImageDataType image_t
     else if (typeid(*this) == typeid(NiftiImageData3DTensor<dataType>))       ss << "NiftiImageData3DTensor";
     else if (typeid(*this) == typeid(NiftiImageData3DDisplacement<dataType>)) ss << "NiftiImageData3DDisplacement";
     else if (typeid(*this) == typeid(NiftiImageData3DDeformation<dataType>))  ss << "NiftiImageData3DDeformation";
+    else if (typeid(*this) == typeid(NiftiImageData3DBSpline<dataType>))      ss << "NiftiImageData3DBSpline";
     ss << ".\n\t\tExpected params: ndim = " << ndim << ", nu = " << nu << ", nt = " << nt;
     if      (intent_code == NIFTI_INTENT_NONE)   ss << ", intent_code = None";
     else if (intent_code == NIFTI_INTENT_VECTOR) ss << ", intent_code = Vector";
     if      (intent_p1 == 0) ss << ", intent_p1 = Deformation";
     else if (intent_p1 == 1) ss << ", intent_p1 = Displacement";
+    else if (intent_p1 == SPLINE_VEL_GRID) ss << ", intent_p1 = Control point grid";
     ss << "\n\t\tActual params:   ndim = " << _nifti_image->ndim << ", nu = " << _nifti_image->nu << ", nt = " << _nifti_image->nt;
     if      (_nifti_image->intent_code == NIFTI_INTENT_NONE)   ss << ", intent_code = None";
     else if (_nifti_image->intent_code == NIFTI_INTENT_VECTOR) ss << ", intent_code = Vector";
     if      (intent_p1 != -1 && _nifti_image->intent_p1 == 0)  ss << ", intent_p1 = Deformation";
     else if (intent_p1 != -1 && _nifti_image->intent_p1 == 1)  ss << ", intent_p1 = Displacement";
+    else if (intent_p1 != -1 && _nifti_image->intent_p1 == SPLINE_VEL_GRID)  ss << ", intent_p1 = Control point grid";
     //std::cout << ss.str() << "\n";
     throw std::runtime_error(ss.str());
 }

src/Registration/cReg/NiftiImageData3DBSpline.cpp

@@ -0,0 +1,78 @@
+/*
+SyneRBI Synergistic Image Reconstruction Framework (SIRF)
+Copyright 2017 - 2020 University College London
+
+This is software developed for the Collaborative Computational
+Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
+(http://www.ccpsynerbi.ac.uk/).
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+http://www.apache.org/licenses/LICENSE-2.0
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+*/
+
+/*!
+\file
+\ingroup Registration
+\brief Class for deformation field transformations.
+
+\author Richard Brown
+\author SyneRBI
+*/
+
+#include "sirf/Reg/NiftiImageData3DBSpline.h"
+#include "sirf/Reg/NiftiImageData3DDeformation.h"
+#include "sirf/NiftyMoMo/BSplineTransformation.h"
+
+using namespace sirf;
+
+template<class dataType>
+void NiftiImageData3DBSpline<dataType>::create_from_3D_image(const NiftiImageData<dataType> &image)
+{
+    NiftiImageData3DTensor<dataType>::create_from_3D_image(image);
+    this->_nifti_image->intent_p1 = SPLINE_VEL_GRID;
+}
+
+template<class dataType>
+NiftiImageData3DDeformation<dataType> NiftiImageData3DBSpline<dataType>::get_as_deformation_field(const NiftiImageData<dataType> &ref) const
+{
+    // Get spacing of reference image
+    float spacing[3];
+    for (unsigned i=0; i<3; ++i)
+        spacing[i] = this->_nifti_image->pixdim[i+1];
+    // Create the NiftyMoMo bspline transformation class
+    NiftyMoMo::BSplineTransformation bspline(ref.clone()->get_raw_nifti_sptr().get(), 1, spacing);
+    // Set the CPG
+    bspline.SetParameters(static_cast<float*>(this->_nifti_image->data), false);
+    // Get the DVF
+    nifti_image *output_def_ptr = bspline.GetDeformationVectorField(ref.get_raw_nifti_sptr().get());
+    return NiftiImageData3DDeformation<dataType>(*output_def_ptr);
+}
+
+template<class dataType>
+NiftiImageData3DBSpline<dataType>*
+NiftiImageData3DBSpline<dataType>::get_inverse_impl_nr(const std::shared_ptr<const NiftiImageData<dataType> >) const
+{
+    throw std::runtime_error("NiftiImageData3DBSpline::get_inverse_impl_nr not yet implemented.");
+}
+
+template<class dataType>
+NiftiImageData3DBSpline<dataType>*
+NiftiImageData3DBSpline<dataType>::get_inverse_impl_vtk(const std::shared_ptr<const NiftiImageData<dataType> >) const
+{
+    throw std::runtime_error("NiftiImageData3DBSpline::get_inverse_impl_vtk not yet implemented.");
+#ifndef SIRF_VTK
+    throw std::runtime_error("Build SIRF with VTK support for this functionality");
+#endif
+}
+
+namespace sirf {
+template class NiftiImageData3DBSpline<float>;
+}

src/Registration/cReg/NiftiImageData3DDeformation.cpp

@@ -54,11 +54,13 @@ void NiftiImageData3DDeformation<dataType>::create_from_3D_image(const NiftiImag
 }
 
 template<class dataType>
-void NiftiImageData3DDeformation<dataType>::create_from_cpp(NiftiImageData3DTensor<dataType> &cpp, const NiftiImageData<dataType> &ref)
+void NiftiImageData3DDeformation<dataType>::create_from_cpp(const NiftiImageData3DTensor<dataType> &cpp, const NiftiImageData<dataType> &ref)
 {
     this->create_from_3D_image(ref);
 
-    reg_spline_getDeformationField(cpp.get_raw_nifti_sptr().get(),
+    auto cpp_clone = cpp.clone();
+
+    reg_spline_getDeformationField(cpp_clone->get_raw_nifti_sptr().get(),
                                    this->_nifti_image.get(),
                                    NULL,
                                    false, //composition