X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;ds=inline;f=appli%2Fexamples%2Fexample_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx;fp=appli%2Fexamples%2Fexample_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx;h=1a8b4769b6eb1684ec685b50f200146c9d39b6cd;hb=b4ed6ddfa7e90e892f07cad9a760961bd4e84e6b;hp=0000000000000000000000000000000000000000;hpb=5e326afe442245572b6c3ec98ebeec8b45f9012f;p=FrontAlgorithms.git diff --git a/appli/examples/example_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx b/appli/examples/example_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx new file mode 100644 index 0000000..1a8b476 --- /dev/null +++ b/appli/examples/example_Image_Dijkstra_EndPointDetection_WithoutVTK.cxx @@ -0,0 +1,249 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 3; +typedef unsigned char TPixel; +typedef float TScalar; + +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::Image< TScalar, Dim > TScalarImage; + +// ------------------------------------------------------------------------- +template< class I > +void ReadImage( typename I::Pointer& image, const std::string& filename ); + +template< class I > +void SaveImage( const I* image, const std::string& filename ); + +template< class I, class O > +void DistanceMap( + const typename I::Pointer& input, typename O::Pointer& output + ); + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 12 ) + { + std::cerr + << "Usage: " << argv[ 0 ] << std::endl + << " input_image" << std::endl + << " seed_x seed_y seed_z" << std::endl + << " output_distancemap" << std::endl + << " output_costmap" << std::endl + << " output_labels" << std::endl + << " output_minimum_spanning_tree" << std::endl + << " output_endpoints" << std::endl + << " output_bifurcations" << std::endl + << " output_branches" << std::endl + << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + double seed_x = std::atof( argv[ 2 ] ); + double seed_y = std::atof( argv[ 3 ] ); + double seed_z = std::atof( argv[ 4 ] ); + std::string distancemap_fn = argv[ 5 ]; + std::string output_costmap_fn = argv[ 6 ]; + std::string output_labels_fn = argv[ 7 ]; + std::string mst_output_fn = argv[ 8 ]; + std::string endpoints_output_fn = argv[ 9 ]; + std::string bifurcations_output_fn = argv[ 10 ]; + std::string branches_output_fn = argv[ 11 ]; + + // Read image + TImage::Pointer input_image; + try + { + ReadImage< TImage >( input_image, input_image_fn ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr + << "Error caught while reading \"" + << input_image_fn << "\": " << err + << std::endl; + return( 1 ); + + } // yrt + + // Get seed + double p[ 3 ]; + TImage::PointType pnt; + TImage::IndexType seed; + pnt[ 0 ] = seed_x; + pnt[ 1 ] = seed_y; + pnt[ 2 ] = seed_z; + input_image->TransformPhysicalPointToIndex( pnt, seed ); + + // Compute squared distance map + TScalarImage::Pointer dmap; + DistanceMap< TImage, TScalarImage >( input_image, dmap ); + + // Prepare cost conversion function + typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction; + TFunction::Pointer function = TFunction::New( ); + + // Prepare Dijkstra filter + typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter; + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( dmap ); + filter->SetConversionFunction( function ); + filter->SetNeighborhoodOrder( 2 ); + filter->StopAtOneFrontOff( ); + filter->AddSeed( seed, TScalar( 0 ) ); + + // Go! + std::time_t start, end; + std::time( &start ); + filter->Update( ); + std::time( &end ); + std::cout + << "Extraction time = " + << std::difftime( end, start ) + << " s." << std::endl; + + // Outputs + const TFilter::TOutputImage* accumulated_costs = filter->GetOutput( ); + const TFilter::TLabelImage* labeled_image = filter->GetLabelImage( ); + const TFilter::TMinimumSpanningTree* mst = filter->GetMinimumSpanningTree( ); + const TFilter::TUniqueVertices* endpoints = filter->GetEndPoints( ); + const TFilter::TUniqueVertices* bifurcations = filter->GetBifurcations( ); + const TFilter::TBranches* branches = filter->GetBranches( ); + unsigned long nBranches = filter->GetNumberOfBranches( ); + + // Save outputs + SaveImage( dmap.GetPointer( ), distancemap_fn ); + SaveImage( accumulated_costs, output_costmap_fn ); + SaveImage( labeled_image, output_labels_fn ); + + fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::Pointer + mst_writer = + fpa::IO::MinimumSpanningTreeWriter< TFilter::TMinimumSpanningTree >::New( ); + mst_writer->SetInput( mst ); + mst_writer->SetFileName( mst_output_fn ); + mst_writer->Update( ); + + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer + endpoints_writer = + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( ); + endpoints_writer->SetInput( endpoints ); + endpoints_writer->SetFileName( endpoints_output_fn ); + endpoints_writer->Update( ); + + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::Pointer + bifurcations_writer = + fpa::IO::UniqueValuesContainerWriter< TFilter::TUniqueVertices >::New( ); + bifurcations_writer->SetInput( bifurcations ); + bifurcations_writer->SetFileName( bifurcations_output_fn ); + bifurcations_writer->Update( ); + + fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::Pointer + branches_writer = + fpa::IO::MatrixValuesContainerWriter< TFilter::TBranches >::New( ); + branches_writer->SetInput( branches ); + branches_writer->SetFileName( branches_output_fn ); + branches_writer->Update( ); + + return( 0 ); +} + +// ------------------------------------------------------------------------- +template< class I > +void ReadImage( typename I::Pointer& image, const std::string& filename ) +{ + typename itk::ImageFileReader< I >::Pointer reader = + itk::ImageFileReader< I >::New( ); + reader->SetFileName( filename ); + reader->Update( ); + + typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax = + itk::MinimumMaximumImageCalculator< I >::New( ); + minmax->SetImage( reader->GetOutput( ) ); + minmax->Compute( ); + double vmin = double( minmax->GetMinimum( ) ); + double vmax = double( minmax->GetMaximum( ) ); + + typename itk::ShiftScaleImageFilter< I, I >::Pointer shift = + itk::ShiftScaleImageFilter< I, I >::New( ); + shift->SetInput( reader->GetOutput( ) ); + shift->SetScale( vmax - vmin ); + shift->SetShift( vmin / ( vmax - vmin ) ); + shift->Update( ); + + image = shift->GetOutput( ); + image->DisconnectPipeline( ); +} + +// ------------------------------------------------------------------------- +template< class I > +void SaveImage( const I* image, const std::string& filename ) +{ + typename itk::ImageFileWriter< I >::Pointer writer = + itk::ImageFileWriter< I >::New( ); + writer->SetInput( image ); + writer->SetFileName( filename ); + try + { + writer->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr + << "Error saving \"" << filename << "\": " << err + << std::endl; + + } // yrt +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void DistanceMap( + const typename I::Pointer& input, typename O::Pointer& output + ) +{ + typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap = + itk::SignedMaurerDistanceMapImageFilter< I, O >::New( ); + dmap->SetInput( input ); + dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) ); + dmap->InsideIsPositiveOn( ); + dmap->SquaredDistanceOn( ); + dmap->UseImageSpacingOn( ); + + std::time_t start, end; + std::time( &start ); + dmap->Update( ); + std::time( &end ); + std::cout + << "Distance map time = " + << std::difftime( end, start ) + << " s." << std::endl; + + output = dmap->GetOutput( ); + output->DisconnectPipeline( ); +} + +// eof - $RCSfile$ +