diff --git a/ants/lib/CMakeLists.txt b/ants/lib/CMakeLists.txt index 2cecc696..60da6aed 100644 --- a/ants/lib/CMakeLists.txt +++ b/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}) diff --git a/ants/lib/LOCAL_ANTSIntegrateVelocityField.cxx b/ants/lib/LOCAL_ANTSIntegrateVelocityField.cxx deleted file mode 100644 index bce8a79b..00000000 --- a/ants/lib/LOCAL_ANTSIntegrateVelocityField.cxx +++ /dev/null @@ -1,101 +0,0 @@ -#include -#include - -#include -#include -#include -#include - -#include "itkTimeVaryingVelocityFieldIntegrationImageFilter.h" -#include "itkTimeVaryingVelocityFieldTransform.h" -#include "itkImageFileWriter.h" -#include "itkImage.h" - -// NEED THIS INCLUDE FOR WRAPPING/UNWRAPPING -#include "LOCAL_antsImage.h" - -// NEED THIS FOR PYBIND11 -namespace py = pybind11; - - -template< class ImageType > -void integrateVelocityField( - py::capsule & r_refimg, - std::string r_velocity, - std::string r_deformation, - float r_t0, - float r_t1, - float r_dt ) -{ - using PixelType = typename ImageType::PixelType; - enum { Dimension = ImageType::ImageDimension }; - PixelType starttime = static_cast< PixelType >( r_t0 ); - PixelType finishtime = static_cast< PixelType >( r_t1 ); - PixelType dT = static_cast< PixelType >( r_dt ); - using VectorType = itk::Vector; - - typedef typename ImageType::Pointer ImagePointerType; - typename ImageType::Pointer input = as< ImageType >( r_refimg ); - - using DisplacementFieldType = itk::Image; - using TimeVaryingVelocityFieldType = itk::Image; - using tvt = TimeVaryingVelocityFieldType; - typename tvt::Pointer timeVaryingVelocity; - - using ReaderType = itk::ImageFileReader; - typename ReaderType::Pointer reader = ReaderType::New(); - reader->SetFileName( r_velocity.c_str() ); - reader->Update(); - timeVaryingVelocity = reader->GetOutput(); - - VectorType zero; - zero.Fill(0); - typename DisplacementFieldType::Pointer deformation = DisplacementFieldType::New(); - deformation->SetRegions( input->GetLargestPossibleRegion() ); - deformation->SetSpacing( input->GetSpacing() ); - deformation->SetOrigin( input->GetOrigin() ); - deformation->SetDirection( input->GetDirection() ); - deformation->Allocate(); - deformation->FillBuffer( zero ); - - if( starttime < 0 ) - { - starttime = 0; - } - if( starttime > 1 ) - { - starttime = 1; - } - if( finishtime < 0 ) - { - finishtime = 0; - } - if( finishtime > 1 ) - { - finishtime = 1; - } - - using IntegratorType = itk::TimeVaryingVelocityFieldIntegrationImageFilter; - typename IntegratorType::Pointer integrator = IntegratorType::New(); - integrator->SetInput( timeVaryingVelocity ); - integrator->SetLowerTimeBound( starttime ); - integrator->SetUpperTimeBound( finishtime ); - integrator->SetNumberOfIntegrationSteps( static_cast( std::round( 1.0 / dT ) ) ); - integrator->Update(); - - using WriterType = itk::ImageFileWriter; - typename WriterType::Pointer writer = WriterType::New(); - writer->SetFileName( r_deformation.c_str() ); - writer->SetInput( integrator->GetOutput() ); - writer->Update(); -// WriteImage( integrator->GetOutput(), r_deformation.c_str() ); - return; -} - - -// DECLARE OUR FUNCTION FOR APPROPRIATE TYPES -PYBIND11_MODULE(integrateVelocityField, m) -{ - m.def("integrateVelocityField2D", &integrateVelocityField>); - m.def("integrateVelocityField3D", &integrateVelocityField>); -} diff --git a/ants/lib/LOCAL_integrateVelocityField.cxx b/ants/lib/LOCAL_integrateVelocityField.cxx new file mode 100644 index 00000000..ab4c1159 --- /dev/null +++ b/ants/lib/LOCAL_integrateVelocityField.cxx @@ -0,0 +1,123 @@ + +#include +#include + +#include +#include +#include + +#include "itkImage.h" +#include "itkVector.h" +#include "itkTimeVaryingVelocityFieldIntegrationImageFilter.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkImageRegionConstIteratorWithIndex.h" + +#include "LOCAL_antsImage.h" + +namespace py = pybind11; + +template +py::capsule integrateVelocityField( py::capsule & antsVelocityField, + float lowerBound, + float upperBound, + unsigned int numberOfIntegrationSteps ) +{ + using RealType = float; + + using ANTsVelocityFieldType = itk::VectorImage; + using ANTsVelocityFieldPointerType = typename ANTsVelocityFieldType::Pointer; + + using ANTsFieldType = itk::VectorImage; + using ANTsFieldPointerType = typename ANTsFieldType::Pointer; + + using VectorType = itk::Vector; + + using ITKVelocityFieldType = itk::Image; + using ITKVelocityFieldPointerType = typename ITKVelocityFieldType::Pointer; + using ITKFieldType = itk::Image; + + using IteratorType = itk::ImageRegionIteratorWithIndex; + using ConstIteratorType = itk::ImageRegionConstIteratorWithIndex; + + ANTsVelocityFieldPointerType inputVelocityField = as( 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; + 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>); +} + diff --git a/ants/utils/__init__.py b/ants/utils/__init__.py index f532caf6..d5403db5 100644 --- a/ants/utils/__init__.py +++ b/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 * diff --git a/ants/utils/ants_integrate_velocity_field.py b/ants/utils/ants_integrate_velocity_field.py deleted file mode 100644 index d731c49c..00000000 --- a/ants/utils/ants_integrate_velocity_field.py +++ /dev/null @@ -1,79 +0,0 @@ -__all__ = ["integrate_velocity_field"] - -from .process_args import _int_antsProcessArguments -from .. import utils - - -def integrate_velocity_field( - reference_image, - velocity_field_filename, - deformation_field_filename, - time_0=0, - time_1=1, - delta_time=0.01, -): - """ - Integrate a velocityfield - - ANTsR function: `integrateVelocityField` - - Arguments - --------- - reference_image : ANTsImage - Reference image domain, same as velocity field space - - velocity_field_filename : string - Filename to velocity field, output from ants.registration - - deformation_field_filename : string - Filename to output deformation field - - time_0 : scalar - Typically one or zero but can take intermediate values - - time_1 : scalar - Typically one or zero but can take intermediate values - - delta_time : scalar - Time step value in zero to one; typically 0.01 - - Returns - ------- - None - - Example - ------- - >>> import ants - >>> fi = ants.image_read( ants.get_data( "r16" ) ) - >>> mi = ants.image_read( ants.get_data( "r27" ) ) - >>> mytx2 = ants.registration( fi, mi, "TV[2]" ) - >>> ants.integrate_velocity_field( fi, mytx2['velocityfield'][0], "/tmp/def.nii.gz", 0, 1, 0.01 ) - >>> mydef = ants.apply_transforms( fi, mi, ["/tmp/def.nii.gz", mytx2['fwdtransforms'][1]] ) - >>> ants.image_mutual_information(fi,mi) - >>> ants.image_mutual_information(fi,mytx2['warpedmovout']) - >>> ants.image_mutual_information(fi,mydef) - >>> ants.integrate_velocity_field( fi, mytx2['velocityfield'][0], "/tmp/defi.nii.gz", 1, 0, 0.5 ) - >>> mydefi = ants.apply_transforms( mi, fi, [ mytx2['fwdtransforms'][1], "/tmp/defi.nii.gz" ] ) - >>> ants.image_mutual_information(mi,mydefi) - >>> ants.image_mutual_information(mi,mytx2['warpedfixout']) - """ - - libfn = utils.get_lib_fn("integrateVelocityField") - if reference_image.dimension == 2: - libfn.integrateVelocityField2D( - reference_image.pointer, - velocity_field_filename, - deformation_field_filename, - time_0, - time_1, - delta_time, - ) - if reference_image.dimension == 3: - libfn.integrateVelocityField3D( - reference_image.pointer, - velocity_field_filename, - deformation_field_filename, - time_0, - time_1, - delta_time, - ) diff --git a/ants/utils/integrate_velocity_field.py b/ants/utils/integrate_velocity_field.py new file mode 100644 index 00000000..f357c107 --- /dev/null +++ b/ants/utils/integrate_velocity_field.py @@ -0,0 +1,57 @@ + +__all__ = ['integrate_velocity_field'] + +from ..core import ants_image as iio +from .. import utils + + +def integrate_velocity_field(velocity_field, + lower_integration_bound=0.0, + upper_integration_bound=1.0, + number_of_integration_steps=10): + """ + Integrate velocity field. + + Arguments + --------- + velocity_field : ANTsImage velocity field + time-varying displacement field + + lower_integration_bound: float + Lower time bound for integration in [0, 1] + + upper_integration_bound: float + Upper time bound for integration in [0, 1] + + number_of_integation_steps: integer + Number of integration steps used in the Runge-Kutta solution + + Example + ------- + >>> import ants + >>> fi = ants.image_read( ants.get_data( "r16" ) ) + >>> mi = ants.image_read( ants.get_data( "r27" ) ) + >>> reg = ants.registration(fi, mi, "TV[2]") + >>> velocity_field = ants.image_read(reg['velocityfield'][0]) + >>> field = ants.integrate_velocity_field(velocity_field, 0.0, 1.0, 10) + """ + + if lower_integration_bound < 0.0: + lower_integration_bound = 0.0 + elif lower_integration_bound > 1.0: + lower_integration_bound = 1.0 + + if upper_integration_bound < 0.0: + upper_integration_bound = 0.0 + elif upper_integration_bound > 1.0: + upper_integration_bound = 1.0 + + libfn = utils.get_lib_fn('integrateVelocityFieldD%i' % (velocity_field.dimension-1)) + integrated_field = libfn(velocity_field.pointer, lower_integration_bound, + upper_integration_bound, number_of_integration_steps) + + new_image = iio.ANTsImage(pixeltype='float', dimension=(velocity_field.dimension-1), + components=(velocity_field.dimension-1), pointer=integrated_field).clone('float') + return new_image + +