Merge pull request #414 from ANTsX/integrateVelocityField

ENH:  Integrate velocity field.

ants/lib/CMakeLists.txt

@@ -56,6 +56,7 @@ pybind11_add_module(fsl2antstransform LOCAL_fsl2antstransform.cxx)
 pybind11_add_module(getNeighborhoodMatrix LOCAL_getNeighborhoodMatrix.cxx)
 pybind11_add_module(hausdorffDistance LOCAL_hausdorffDistance.cxx)
 pybind11_add_module(histogramMatchImage LOCAL_histogramMatchImages.cxx)
+pybind11_add_module(integrateVelocityField LOCAL_integrateVelocityField.cxx)
 pybind11_add_module(invertDisplacementField LOCAL_invertDisplacementField.cxx)
 pybind11_add_module(labelOverlapMeasures LOCAL_labelOverlapMeasures.cxx)
 pybind11_add_module(labelStats LOCAL_labelStats.cxx)
@@ -72,7 +73,6 @@ pybind11_add_module(simulateDisplacementField LOCAL_simulateDisplacementField.cx
 pybind11_add_module(sliceImage LOCAL_sliceImage.cxx)
 pybind11_add_module(SmoothImage LOCAL_SmoothImage.cxx)
 pybind11_add_module(weingartenImageCurvature LOCAL_weingartenImageCurvature.cxx)
-pybind11_add_module(integrateVelocityField LOCAL_ANTSIntegrateVelocityField.cxx)
 
 ## WRAP ##
 pybind11_add_module(antsAffineInitializer antscore/antsAffineInitializer.cxx WRAP_antsAffineInitializer.cxx)
@@ -130,6 +130,7 @@ target_link_libraries(fsl2antstransform PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(getNeighborhoodMatrix PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(hausdorffDistance PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(histogramMatchImage PRIVATE ${ITK_LIBRARIES})
+target_link_libraries(integrateVelocityField PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(invertDisplacementField PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(labelOverlapMeasures PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(labelStats PRIVATE ${ITK_LIBRARIES})
@@ -143,7 +144,6 @@ target_link_libraries(reorientImage2 PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(rgbToVector PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(sccaner PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(simulateDisplacementField PRIVATE ${ITK_LIBRARIES})
-target_link_libraries(integrateVelocityField PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(sliceImage PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(SmoothImage PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(weingartenImageCurvature PRIVATE ${ITK_LIBRARIES})

ants/lib/LOCAL_integrateVelocityField.cxx

@@ -0,0 +1,123 @@
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <exception>
+#include <vector>
+#include <string>
+
+#include "itkImage.h"
+#include "itkVector.h"
+#include "itkTimeVaryingVelocityFieldIntegrationImageFilter.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+
+#include "LOCAL_antsImage.h"
+
+namespace py = pybind11;
+
+template <unsigned int Dimension>
+py::capsule integrateVelocityField( py::capsule & antsVelocityField,
+                                    float lowerBound,
+                                    float upperBound,
+                                    unsigned int numberOfIntegrationSteps )
+{
+  using RealType = float;
+
+  using ANTsVelocityFieldType = itk::VectorImage<RealType, Dimension+1>;
+  using ANTsVelocityFieldPointerType = typename ANTsVelocityFieldType::Pointer;
+
+  using ANTsFieldType = itk::VectorImage<RealType, Dimension>;
+  using ANTsFieldPointerType = typename ANTsFieldType::Pointer;
+
+  using VectorType = itk::Vector<RealType, Dimension>;
+
+  using ITKVelocityFieldType = itk::Image<VectorType, Dimension+1>;
+  using ITKVelocityFieldPointerType = typename ITKVelocityFieldType::Pointer;
+  using ITKFieldType = itk::Image<VectorType, Dimension>;
+
+  using IteratorType = itk::ImageRegionIteratorWithIndex<ITKVelocityFieldType>;
+  using ConstIteratorType = itk::ImageRegionConstIteratorWithIndex<ITKFieldType>;
+
+  ANTsVelocityFieldPointerType inputVelocityField = as<ANTsVelocityFieldType>( antsVelocityField );
+
+  typename ITKVelocityFieldType::PointType fieldOrigin;
+  typename ITKVelocityFieldType::SpacingType fieldSpacing;
+  typename ITKVelocityFieldType::SizeType fieldSize;
+  typename ITKVelocityFieldType::DirectionType fieldDirection;
+
+  for( unsigned int d = 0; d < Dimension + 1; d++ )
+    {
+    fieldOrigin[d] = inputVelocityField->GetOrigin()[d];
+    fieldSpacing[d] = inputVelocityField->GetSpacing()[d];
+    fieldSize[d] = inputVelocityField->GetRequestedRegion().GetSize()[d];
+    for( unsigned int e = 0; e < Dimension + 1; e++ )
+      {
+      fieldDirection(d, e) = inputVelocityField->GetDirection()(d, e);
+      }
+    }
+
+  ITKVelocityFieldPointerType inputITKVelocityField = ITKVelocityFieldType::New();
+  inputITKVelocityField->SetOrigin( fieldOrigin );
+  inputITKVelocityField->SetRegions( fieldSize );
+  inputITKVelocityField->SetSpacing( fieldSpacing );
+  inputITKVelocityField->SetDirection( fieldDirection );
+  inputITKVelocityField->Allocate();
+
+  IteratorType It( inputITKVelocityField, 
+                   inputITKVelocityField->GetRequestedRegion() );
+  for( It.GoToBegin(); !It.IsAtEnd(); ++It )
+    {
+    VectorType vector;
+
+    typename ANTsFieldType::PixelType antsVector = inputVelocityField->GetPixel( It.GetIndex() );
+    for( unsigned int d = 0; d < Dimension; d++ )
+      {
+      vector[d] = antsVector[d];
+      }
+    It.Set( vector );
+    }
+
+  using IntegratorType = itk::TimeVaryingVelocityFieldIntegrationImageFilter<ITKVelocityFieldType, ITKFieldType>;
+  typename IntegratorType::Pointer integrator = IntegratorType::New();
+
+  integrator->SetInput( inputITKVelocityField );
+  integrator->SetLowerTimeBound( lowerBound );
+  integrator->SetUpperTimeBound( upperBound );
+  integrator->SetNumberOfIntegrationSteps( numberOfIntegrationSteps );
+  integrator->Update();
+
+  //////////////////////////
+  //
+  //  Now convert back to vector image type.
+  //
+
+  ANTsFieldPointerType antsField = ANTsFieldType::New();
+  antsField->CopyInformation( integrator->GetOutput() );
+  antsField->SetRegions( integrator->GetOutput()->GetRequestedRegion() );
+  antsField->SetVectorLength( Dimension );
+  antsField->Allocate();
+
+  ConstIteratorType ItI( integrator->GetOutput(),
+    integrator->GetOutput()->GetRequestedRegion() );
+  for( ItI.GoToBegin(); !ItI.IsAtEnd(); ++ItI )
+    {
+    VectorType data = ItI.Value();
+
+    typename ANTsFieldType::PixelType antsVector( Dimension );
+    for( unsigned int d = 0; d < Dimension; d++ )
+      {
+      antsVector[d] = data[d];
+      }
+    antsField->SetPixel( ItI.GetIndex(), antsVector );
+    }
+
+  return wrap< ANTsFieldType >( antsField );
+}
+
+PYBIND11_MODULE(integrateVelocityField, m)
+{
+  m.def("integrateVelocityFieldD2", &integrateVelocityField<2>);
+  m.def("integrateVelocityFieldD3", &integrateVelocityField<3>);
+}
+

ants/utils/__init__.py

@@ -17,6 +17,7 @@
 from .image_to_cluster_images import *
 from .iMath import *
 from .impute import *
+from .integrate_velocity_field import *
 from .invariant_image_similarity import *
 from .invert_displacement_field import *
 from .label_clusters import *
@@ -37,5 +38,4 @@
 from .slice_image import *
 from .smooth_image import *
 from .threshold_image import *
-from .ants_integrate_velocity_field import *
 from .weingarten_image_curvature import *