X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2Fexamples%2Fexample_ImageAlgorithm_Skeletonization.cxx;fp=appli%2Fexamples%2Fexample_ImageAlgorithm_Skeletonization.cxx;h=0000000000000000000000000000000000000000;hb=9c4035dc18f4e44a92cda3e3c1c6e583c94ef06b;hp=22b85e773ec9b5ead755e23d76f51a3c8b44fe04;hpb=e832f438861464a8d521464b61be7e31c896adb8;p=FrontAlgorithms.git diff --git a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx deleted file mode 100644 index 22b85e7..0000000 --- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx +++ /dev/null @@ -1,257 +0,0 @@ -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include - -#include -#include -#include -#include - -// ------------------------------------------------------------------------- -const unsigned int Dim = 3; -typedef short TPixel; -typedef float TScalar; -typedef itk::Image< TPixel, Dim > TImage; -typedef itk::Image< TScalar, Dim > TDistanceMap; -typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; -typedef itk::ImageToVTKImageFilter< TDistanceMap > TVTKDistanceMap; - -typedef itk::ImageFileReader< TImage > TImageReader; -typedef itk::ImageFileWriter< TImage > TImageWriter; -typedef itk::ImageFileWriter< TDistanceMap > TDistanceMapWriter; -typedef fpa::Image::RegionGrowWithMultipleThresholds< TImage > TSegmentor; -typedef -itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap > -TDistanceFilter; -typedef -fpa::Image::DijkstraWithSphereBacktracking< TDistanceMap, TScalar > -TDijkstra; - -typedef fpa::VTK::ImageMPR TMPR; -typedef fpa::VTK::Image3DObserver< TSegmentor, vtkRenderWindow > TSegmentorObs; -typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs; - -// ------------------------------------------------------------------------- -int main( int argc, char* argv[] ) -{ - if( argc < 8 ) - { - std::cerr - << "Usage: " << argv[ 0 ] - << " input_image thr_0 thr_1 step" - << " output_segmentation output_distancemap output_dijkstra" - << " visual_debug" - << std::endl; - return( 1 ); - - } // fi - std::string input_image_fn = argv[ 1 ]; - TPixel thr_0 = TPixel( std::atof( argv[ 2 ] ) ); - TPixel thr_1 = TPixel( std::atof( argv[ 3 ] ) ); - unsigned int step = std::atoi( argv[ 4 ] ); - std::string output_segmentation_fn = argv[ 5 ]; - std::string output_distancemap_fn = argv[ 6 ]; - std::string output_dijkstra_fn = argv[ 7 ]; - bool visual_debug = false; - if( argc > 8 ) - visual_debug = ( std::atoi( argv[ 8 ] ) == 1 ); - - // Read image - TImageReader::Pointer input_image_reader = TImageReader::New( ); - input_image_reader->SetFileName( input_image_fn ); - try - { - input_image_reader->Update( ); - } - catch( itk::ExceptionObject& err ) - { - std::cerr << "Error caught: " << err << std::endl; - return( 1 ); - - } // yrt - TImage::ConstPointer input_image = input_image_reader->GetOutput( ); - - // Show image - TVTKImage::Pointer vtk_image = TVTKImage::New( ); - vtk_image->SetInput( input_image ); - vtk_image->Update( ); - - TMPR view; - view.SetBackground( 0.3, 0.2, 0.8 ); - view.SetSize( 800, 800 ); - view.SetImage( vtk_image->GetOutput( ) ); - - // Wait for a seed to be given - while( view.GetNumberOfSeeds( ) == 0 ) - view.Start( ); - - // Get seed from user - double seed[ 3 ]; - view.GetSeed( 0, seed ); - TImage::PointType seed_pnt; - seed_pnt[ 0 ] = seed[ 0 ]; - seed_pnt[ 1 ] = seed[ 1 ]; - seed_pnt[ 2 ] = seed[ 2 ]; - TImage::IndexType seed_idx; - input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx ); - - // Segment input image - TSegmentor::Pointer segmentor = TSegmentor::New( ); - segmentor->AddThresholds( thr_0, thr_1, step ); - segmentor->AddSeed( seed_idx, 0 ); - segmentor->SetInput( input_image ); - segmentor->SetNeighborhoodOrder( 1 ); - segmentor->SetDifferenceThreshold( double( 3 ) ); - segmentor->SetInsideValue( TPixel( 0 ) ); // WARNING: This is to optimize - segmentor->SetOutsideValue( TPixel( 1 ) ); // distance map calculation - if( visual_debug ) - { - // Configure observer - TSegmentorObs::Pointer obs = TSegmentorObs::New( ); - obs->SetRenderWindow( view.GetWindow( ) ); - segmentor->AddObserver( itk::AnyEvent( ), obs ); - segmentor->ThrowEventsOn( ); - } - else - segmentor->ThrowEventsOff( ); - std::clock_t start = std::clock( ); - segmentor->Update( ); - std::clock_t end = std::clock( ); - double seconds = double( end - start ) / double( CLOCKS_PER_SEC ); - std::cout << "Segmentation time = " << seconds << std::endl; - - // Extract distance map - TDistanceFilter::Pointer distanceMap = TDistanceFilter::New( ); - distanceMap->SetInput( segmentor->GetOutput( ) ); - distanceMap->InputIsBinaryOn( ); - distanceMap->SquaredDistanceOn( ); - distanceMap->UseImageSpacingOn( ); - start = std::clock( ); - distanceMap->Update( ); - end = std::clock( ); - seconds = double( end - start ) / double( CLOCKS_PER_SEC ); - std::cout << "Distance map time = " << seconds << std::endl; - - TVTKDistanceMap::Pointer vtk_distanceMap = TVTKDistanceMap::New( ); - vtk_distanceMap->SetInput( distanceMap->GetOutput( ) ); - vtk_distanceMap->Update( ); - - vtkSmartPointer< vtkImageMarchingCubes > mc = - vtkSmartPointer< vtkImageMarchingCubes >::New( ); - mc->SetInputData( vtk_distanceMap->GetOutput( ) ); - mc->SetValue( 0, 1e-2 ); - mc->Update( ); - view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 ); - view.Render( ); - - // Extract paths - TDijkstra::Pointer paths = TDijkstra::New( ); - paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) ); - paths->SetInput( distanceMap->GetOutput( ) ); - paths->SetNeighborhoodOrder( 2 ); - - if( visual_debug ) - { - // Configure observer - TDijkstraObs::Pointer obs = TDijkstraObs::New( ); - obs->SetRenderWindow( view.GetWindow( ) ); - paths->AddObserver( itk::AnyEvent( ), obs ); - paths->ThrowEventsOn( ); - } - else - paths->ThrowEventsOff( ); - start = std::clock( ); - paths->Update( ); - end = std::clock( ); - seconds = double( end - start ) / double( CLOCKS_PER_SEC ); - std::cout << "Paths extraction time = " << seconds << std::endl; - - std::cout << "Final seed: " << paths->GetSeed( 0 ) << std::endl; - - // Create polydata - vtkSmartPointer< vtkPoints > points = - vtkSmartPointer< vtkPoints >::New( ); - vtkSmartPointer< vtkCellArray > cells = - vtkSmartPointer< vtkCellArray >::New( ); - vtkSmartPointer< vtkFloatArray > scalars = - vtkSmartPointer< vtkFloatArray >::New( ); - - const TDijkstra::TVertices& endpoints = paths->GetEndPoints( ); - const TDijkstra::TTree& tree = paths->GetFullTree( ); - TDijkstra::TVertices::const_iterator epIt = endpoints.begin( ); - for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId ) - { - double pd = double( epId ) / double( endpoints.size( ) - 1 ); - - TDijkstra::TVertex idx = *epIt; - do - { - TImage::PointType pnt; - input_image->TransformIndexToPhysicalPoint( idx, pnt ); - - points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - scalars->InsertNextTuple1( pd ); - if( idx != *epIt ) - { - cells->InsertNextCell( 2 ); - cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); - cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); - - } // fi - idx = tree.find( idx )->second.first; - - } while( idx != tree.find( idx )->second.first ); - - } // rof - - vtkSmartPointer< vtkPolyData > vtk_tree = - vtkSmartPointer< vtkPolyData >::New( ); - vtk_tree->SetPoints( points ); - vtk_tree->SetLines( cells ); - vtk_tree->GetPointData( )->SetScalars( scalars ); - - view.AddPolyData( vtk_tree ); - view.Render( ); - view.Start( ); - - itk::ImageFileWriter< TImage >::Pointer segmentation_writer = - itk::ImageFileWriter< TImage >::New( ); - segmentation_writer->SetInput( segmentor->GetOutput( ) ); - segmentation_writer->SetFileName( output_segmentation_fn ); - - itk::ImageFileWriter< TDistanceMap >::Pointer dmap_writer = - itk::ImageFileWriter< TDistanceMap >::New( ); - dmap_writer->SetInput( distanceMap->GetOutput( ) ); - dmap_writer->SetFileName( output_distancemap_fn ); - - itk::ImageFileWriter< TDistanceMap >::Pointer dijk_writer = - itk::ImageFileWriter< TDistanceMap >::New( ); - dijk_writer->SetInput( paths->GetOutput( ) ); - dijk_writer->SetFileName( output_dijkstra_fn ); - - try - { - segmentation_writer->Update( ); - dmap_writer->Update( ); - dijk_writer->Update( ); - } - catch( itk::ExceptionObject& err ) - { - std::cerr << "Error caught: " << err << std::endl; - return( -1 ); - - } // yrt - - return( 0 ); -} - -// eof - $RCSfile$