B spline2

src/Registration/NiftyMoMo/CMakeLists.txt

@@ -32,6 +32,3 @@ TARGET_INCLUDE_DIRECTORIES(NiftyMoMo BEFORE PUBLIC "${NIFTYREG_INCLUDE_DIRS}")
 
 # Install
 INSTALL(TARGETS NiftyMoMo DESTINATION ${CMAKE_INSTALL_PREFIX}/lib)
-INSTALL (
-    DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/include/
-    DESTINATION ${CMAKE_INSTALL_PREFIX}/include)

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

src/Registration/cReg/cReg.cpp

@@ -26,6 +26,8 @@ limitations under the License.
 #include "sirf/Reg/NiftiImageData3DTensor.h"
 #include "sirf/Reg/NiftiImageData3DDisplacement.h"
 #include "sirf/Reg/NiftiImageData3DDeformation.h"
+#include "sirf/Reg/NiftiImageData3DBSpline.h"
+#include "sirf/Reg/ControlPointGridToDeformationConverter.h"
 #include "sirf/Reg/NiftyAladinSym.h"
 #include "sirf/Reg/NiftyF3dSym.h"
 #include "sirf/Reg/NiftyResample.h"
@@ -68,6 +70,10 @@ void* cReg_newObject(const char* name)
             return newObjectHandle(std::shared_ptr<NiftiImageData3DDisplacement<float> >(new NiftiImageData3DDisplacement<float>));
         if (strcmp(name, "NiftiImageData3DDeformation") == 0)
             return newObjectHandle(std::shared_ptr<NiftiImageData3DDeformation<float> >(new NiftiImageData3DDeformation<float>));
+        if (strcmp(name, "NiftiImageData3DBSpline") == 0)
+            return newObjectHandle(std::shared_ptr<NiftiImageData3DBSpline<float> >(new NiftiImageData3DBSpline<float>));
+        if (strcmp(name, "ControlPointGridToDeformationConverter") == 0)
+            return newObjectHandle(std::shared_ptr<ControlPointGridToDeformationConverter<float> >(new ControlPointGridToDeformationConverter<float>));
         if (strcmp(name, "NiftyAladinSym") == 0)
             return newObjectHandle(std::shared_ptr<NiftyAladinSym<float> >(new NiftyAladinSym<float>));
         if (strcmp(name, "NiftyF3dSym") == 0)
@@ -159,6 +165,11 @@ void* cReg_objectFromFile(const char* name, const char* filename)
                 sptr(new NiftiImageData3DDeformation<float>(filename));
             return newObjectHandle(sptr);
         }
+        if (strcmp(name, "NiftiImageData3DBSpline") == 0) {
+            std::shared_ptr<NiftiImageData3DBSpline<float> >
+                sptr(new NiftiImageData3DBSpline<float>(filename));
+            return newObjectHandle(sptr);
+        }
         if (strcmp(name, "AffineTransformation") == 0) {
             std::shared_ptr<AffineTransformation<float> >
                 sptr(new AffineTransformation<float>(filename));
@@ -533,6 +544,12 @@ void* cReg_NiftiImageData3DTensor_construct_from_3_components(const char* obj, c
             sptr.reset(new NiftiImageData3DDisplacement<float>(x,y,z));
         else if (strcmp(obj,"NiftiImageData3DDeformation") == 0)
             sptr.reset(new NiftiImageData3DDeformation<float>(x,y,z));
+        else if (strcmp(obj,"NiftiImageData3DBSpline") == 0)
+            sptr.reset(new NiftiImageData3DBSpline<float>(x,y,z));
+        else
+            throw std::runtime_error(
+                    "cReg_NiftiImageData3DTensor_construct_from_3_components, unknown type:" +
+                    std::string(obj));
         return newObjectHandle(sptr);
     }
     CATCH;
@@ -547,6 +564,17 @@ void* cReg_NiftiImageData3DTensor_flip_component(const void *ptr, const int dim)
     }
     CATCH;
 }
+extern "C"
+void* cReg_NiftiImageData3DTensor_get_tensor_component(const void *ptr, const int dim)
+{
+    try {
+        NiftiImageData3DTensor<float>& im = objectFromHandle<NiftiImageData3DTensor<float> >(ptr);
+        std::shared_ptr<NiftiImageData<float> > im_sptr = im.get_tensor_component(dim);
+        std::shared_ptr<NiftiImageData3D<float> > im3D_sptr = std::make_shared<NiftiImageData3D<float> >(*im_sptr);
+        return newObjectHandle(im3D_sptr);
+    }
+    CATCH;
+}
 // -------------------------------------------------------------------------------- //
 //      NiftiImageData3DDeformation
 // -------------------------------------------------------------------------------- //
@@ -610,6 +638,61 @@ void* cReg_NiftiImageData3DDisplacement_create_from_def(const void* def_ptr)
     CATCH;
 }
 
+// -------------------------------------------------------------------------------- //
+//      ControlPointGridToDeformationConverter
+// -------------------------------------------------------------------------------- //
+extern "C"
+void* cReg_CPG2DVF_set_cpg_spacing(const void* converter_ptr, const float spacing_x, const float spacing_y, const float spacing_z)
+{
+    try {
+        ControlPointGridToDeformationConverter<float>& cpg_2_dvf_converter =
+                objectFromHandle<ControlPointGridToDeformationConverter<float> >(converter_ptr);
+        const float spacing[3] = {spacing_x, spacing_y, spacing_z};
+        cpg_2_dvf_converter.set_cpg_spacing(spacing);
+        return new DataHandle;
+    }
+    CATCH;
+}
+extern "C"
+void* cReg_CPG2DVF_set_ref_im(const void* converter_ptr, const void* ref_im_ptr)
+{
+    try {
+        ControlPointGridToDeformationConverter<float>& cpg_2_dvf_converter =
+                objectFromHandle<ControlPointGridToDeformationConverter<float> >(converter_ptr);
+        NiftiImageData<float>& ref_im =
+                objectFromHandle<NiftiImageData<float> >(ref_im_ptr);
+        cpg_2_dvf_converter.set_reference_image(ref_im);
+        return new DataHandle;
+    }
+    CATCH;
+}
+extern "C"
+void* cReg_CPG2DVF_forward(const void* converter_ptr, const void* cpg_ptr)
+{
+    try {
+        ControlPointGridToDeformationConverter<float>& cpg_2_dvf_converter =
+                objectFromHandle<ControlPointGridToDeformationConverter<float> >(converter_ptr);
+        NiftiImageData3DBSpline<float>& cpg =
+                objectFromHandle<NiftiImageData3DBSpline<float> >(cpg_ptr);
+        NiftiImageData3DDeformation<float> def = cpg_2_dvf_converter.forward(cpg);
+        return newObjectHandle(std::make_shared<NiftiImageData3DDeformation<float> >(def));
+    }
+    CATCH;
+}
+extern "C"
+void* cReg_CPG2DVF_backward(const void* converter_ptr, const void* dvf_ptr)
+{
+    try {
+        ControlPointGridToDeformationConverter<float>& cpg_2_dvf_converter =
+                objectFromHandle<ControlPointGridToDeformationConverter<float> >(converter_ptr);
+        NiftiImageData3DDeformation<float>& dvf =
+                objectFromHandle<NiftiImageData3DDeformation<float> >(dvf_ptr);
+        NiftiImageData3DBSpline<float> cpg = cpg_2_dvf_converter.backward(dvf);
+        return newObjectHandle(std::make_shared<NiftiImageData3DBSpline<float> >(cpg));
+    }
+    CATCH;
+}
+
 // -------------------------------------------------------------------------------- //
 //      Registration
 // -------------------------------------------------------------------------------- //
@@ -903,6 +986,8 @@ void* cReg_Transformation_get_as_deformation_field(const void* ptr, const char*
             trans = &objectFromHandle<NiftiImageData3DDisplacement<float> >(ptr);
         else if (strcmp(name,"NiftiImageData3DDeformation") == 0)
             trans = &objectFromHandle<NiftiImageData3DDeformation<float> >(ptr);
+        else if (strcmp(name,"NiftiImageData3DBSpline") == 0)
+            trans = &objectFromHandle<NiftiImageData3DBSpline<float> >(ptr);
         else
             throw std::runtime_error("cReg_Transformation_get_as_deformation_field: type should be affine, disp or def.");