+++ /dev/null
-#include <cstdlib>
-#include <ctime>
-#include <iostream>
-#include <limits>
-#include <string>
-
-#include <itkImage.h>
-
-#include <itkImageFileReader.h>
-#include <itkMinimumMaximumImageCalculator.h>
-#include <itkShiftScaleImageFilter.h>
-
-#include <itkSignedMaurerDistanceMapImageFilter.h>
-
-#include <itkImageFileWriter.h>
-
-#include <fpa/Image/DijkstraWithEndPointDetection.h>
-#include <fpa/Base/Functors/InvertCostFunction.h>
-
-#include <fpa/IO/MinimumSpanningTreeWriter.h>
-#include <fpa/IO/UniqueValuesContainerWriter.h>
-#include <fpa/IO/MatrixValuesContainerWriter.h>
-
-// -------------------------------------------------------------------------
-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$
-
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / blobdiff