Merge pull request #411 from ANTsX/InvertField

Invert field

ants/lib/CMakeLists.txt

@@ -55,6 +55,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(invertDisplacementField LOCAL_invertDisplacementField.cxx)
 pybind11_add_module(labelOverlapMeasures LOCAL_labelOverlapMeasures.cxx)
 pybind11_add_module(labelStats LOCAL_labelStats.cxx)
 pybind11_add_module(mergeChannels LOCAL_mergeChannels.cxx)
@@ -127,6 +128,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(invertDisplacementField PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(labelOverlapMeasures PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(labelStats PRIVATE ${ITK_LIBRARIES})
 target_link_libraries(mergeChannels PRIVATE ${ITK_LIBRARIES})

ants/lib/LOCAL_invertDisplacementField.cxx

@@ -0,0 +1,83 @@
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include <exception>
+#include <vector>
+#include <string>
+
+#include "itkImage.h"
+#include "itkVector.h"
+#include "itkInvertDisplacementFieldImageFilter.h"
+
+#include "LOCAL_antsImage.h"
+
+namespace py = pybind11;
+
+template <unsigned int Dimension>
+py::capsule invertDisplacementField( py::capsule & antsDisplacementField,
+ py::capsule & antsInverseFieldInitialEstimate,
+ unsigned int maximumNumberOfIterations,
+ float meanErrorToleranceThreshold,
+ float maxErrorToleranceThreshold,
+ bool enforceBoundaryCondition )
+{
+ using RealType = float;
+
+ using ANTsFieldType = itk::VectorImage<RealType, Dimension>;
+ using ANTsFieldPointerType = typename ANTsFieldType::Pointer;
+
+ using VectorType = itk::Vector<RealType, Dimension>;
+
+ using ITKFieldType = itk::Image<VectorType, Dimension>;
+ using ITKFieldPointerType = typename ITKFieldType::Pointer;
+
+ using IteratorType = itk::ImageRegionConstIteratorWithIndex<ITKFieldType>;
+
+ ITKFieldPointerType itkDisplacementField = as< ITKFieldType >( antsDisplacementField );
+ ITKFieldPointerType itkInverseFieldInitialEstimate = as< ITKFieldType >( antsInverseFieldInitialEstimate );
+
+ typedef itk::InvertDisplacementFieldImageFilter<ITKFieldType> InverterType;
+ typename InverterType::Pointer inverter = InverterType::New();
+
+ inverter->SetInput( itkDisplacementField );
+ inverter->SetInverseFieldInitialEstimate( itkInverseFieldInitialEstimate );
+ inverter->SetMaximumNumberOfIterations( maximumNumberOfIterations );
+ inverter->SetMeanErrorToleranceThreshold( meanErrorToleranceThreshold );
+ inverter->SetMaxErrorToleranceThreshold( maxErrorToleranceThreshold );
+ inverter->SetEnforceBoundaryCondition( enforceBoundaryCondition );
+
+ //////////////////////////
+ //
+ // Now convert back to vector image type.
+ //
+
+ ANTsFieldPointerType antsField = ANTsFieldType::New();
+ antsField->CopyInformation( inverter->GetOutput() );
+ antsField->SetRegions( inverter->GetOutput()->GetRequestedRegion() );
+ antsField->SetVectorLength( Dimension );
+ antsField->Allocate();
+
+ IteratorType ItI( inverter->GetOutput(),
+ inverter->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(invertDisplacementField, m)
+{
+ m.def("invertDisplacementFieldF2", &invertDisplacementField<2>);
+ m.def("invertDisplacementFieldF3", &invertDisplacementField<3>);
+}
+

ants/lib/__init__.py

@@ -16,6 +16,7 @@
 from .getNeighborhoodMatrix import *
 from .hausdorffDistance import *
 from .histogramMatchImage import *
+from .invertDisplacementField import *
 # from .invariantImageSimilarity import *
 from .labelOverlapMeasures import *
 from .labelStats import *

ants/registration/interface.py

@@ -1161,8 +1161,17 @@ def registration(
  args.append("[NA,NA]")
  # ------------------------------------------------------------
  elif "antsRegistrationSyN" in type_of_transform:
+
+ do_quick = False
+ if "Quick" in type_of_transform:
+ do_quick = True
+
  subtype_of_transform = "s"
  spline_distance = 26
+ metric_parameter = 4
+ if do_quick:
+ metric_parameter = 32
+
  if "[" in type_of_transform and "]" in type_of_transform:
  subtype_of_transform = type_of_transform.split("[")[1].split(
  "]"
@@ -1172,13 +1181,14 @@ def registration(
  subtype_of_transform = subtype_of_transform_args[0]
  if not ( subtype_of_transform == "b"
  or subtype_of_transform == "br"
- or subtype_of_transform == "bo" ):
- raise ValueError("Extra parameters are only valid for B-spline SyN transform.")
- spline_distance = subtype_of_transform_args[1]
-
- do_quick = False
- if "Quick" in type_of_transform:
- do_quick = True
+ or subtype_of_transform == "bo" 
+ or subtype_of_transform == "s"
+ or subtype_of_transform == "sr"
+ or subtype_of_transform == "so" ):
+ raise ValueError("Extra parameters are only valid for 's' or 'b' SyN transforms.")
+ metric_parameter = subtype_of_transform_args[1]
+ if len(subtype_of_transform_args) > 2:
+ spline_distance = subtype_of_transform_args[2]
 
  do_repro = False
  if "Repro" in type_of_transform:
@@ -1204,16 +1214,16 @@ def registration(
 
  if do_quick == True:
  syn_convergence = "[100x70x50x0,1e-6,10]"
- syn_metric = "MI[%s,%s,1,32]" % (f, m)
+ syn_metric = "MI[%s,%s,1,%s]" % (f, m, metric_parameter)
  else:
  syn_convergence = "[100x70x50x20,1e-6,10]"
- syn_metric = "CC[%s,%s,1,4]" % (f, m)
+ syn_metric = "CC[%s,%s,1,%s]" % (f, m, metric_parameter)
  syn_shrink_factors = "8x4x2x1"
  syn_smoothing_sigmas = "3x2x1x0vox"
 
  if do_quick == True and do_repro == True:
  syn_convergence = "[100x70x50x0,1e-6,10]"
- syn_metric = "CC[%s,%s,1,2]" % (f, m)
+ syn_metric = "CC[%s,%s,1,%s]" % (f, m, metric_parameter)
 
  if random_seed is None and do_repro == True:
  random_seed = str( 1 )

ants/utils/__init__.py

@@ -17,6 +17,7 @@
 from .iMath import *
 from .impute import *
 from .invariant_image_similarity import *
+from .invert_displacement_field import *
 from .label_clusters import *
 from .label_image_centroids import *
 from .label_overlap_measures import *

ants/utils/invert_displacement_field.py

@@ -0,0 +1,52 @@
+
+__all__ = ['invert_displacement_field']
+
+from ..core import ants_image as iio
+from .. import utils
+
+
+def invert_displacement_field(displacement_field, 
+ inverse_field_initial_estimate, 
+ maximum_number_of_iterations=20, 
+ mean_error_tolerance_threshold=0.001, 
+ max_error_tolerance_threshold=0.1,
+ enforce_boundary_condition=True):
+ """
+ Invert displacement field.
+
+ Arguments
+ ---------
+ displacement_field : ANTsImage displacement field
+ displacement field
+
+ inverse_field_initial_estimate : ANTsImage displacement field
+ initial guess
+
+ maximum_number_of_iterations : integer
+ number of iterations
+
+ mean_error_tolerance_threshold : float
+ mean error tolerance threshold
+
+ max_error_tolerance_threshold : float
+ max error tolerance threshold
+
+ enforce_boundary_condition : bool
+ enforce stationary boundary condition
+ 
+
+ Example
+ -------
+ >>> import ants
+ """
+
+ libfn = utils.get_lib_fn('invertDisplacementFieldF%i' % displacement_field.dimension)
+ inverse_field = libfn(displacement_field.pointer, inverse_field_initial_estimate.pointer, 
+ maximum_number_of_iterations, mean_error_tolerance_threshold, 
+ max_error_tolerance_threshold, enforce_boundary_condition)
+
+ new_image = iio.ANTsImage(pixeltype='float', dimension=displacement_field.dimension, 
+ components=displacement_field.dimension, pointer=inverse_field).clone('float')
+ return new_image
+
+