X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2Fexamples%2Fexample_ImageAlgorithm_Skeletonization.cxx;h=22b85e773ec9b5ead755e23d76f51a3c8b44fe04;hb=b70a564ee2d7bc180b77a05c37ab431ab9c393e7;hp=156469ef0aeb414f70ec6803503a67448d3f22aa;hpb=73e3158dab6ffe30c215b2899fd80ee868f330ce;p=FrontAlgorithms.git diff --git a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx index 156469e..22b85e7 100644 --- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx +++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx @@ -2,62 +2,53 @@ #include #include +#include #include #include #include +#include #include -#include -#include +#include + #include +#include #include #include -/* - #include - #include - - #include - - - #include - #include - #include - #include - #include - #include - #include - #include - -*/ - // ------------------------------------------------------------------------- const unsigned int Dim = 3; typedef short TPixel; -typedef double TScalar; +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 itk::ImageToVTKImageFilter< TImage > TVTKImage; - typedef itk::ImageFileReader< TImage > TImageReader; - typedef itk::ImageFileWriter< TImage > TImageWriter; - typedef - fpa::Image::RegionGrowWithMultipleThresholds< TImage > - TSegmentor; - typedef - - TObserver; -*/ +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 < 6 ) + if( argc < 8 ) { std::cerr << "Usage: " << argv[ 0 ] - << " input_image thr_0 thr_1 step output_image" + << " input_image thr_0 thr_1 step" + << " output_segmentation output_distancemap output_dijkstra" << " visual_debug" << std::endl; return( 1 ); @@ -67,14 +58,15 @@ int main( int argc, char* argv[] ) 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_image_fn = argv[ 5 ]; + 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 > 6 ) - visual_debug = ( std::atoi( argv[ 6 ] ) == 1 ); + if( argc > 8 ) + visual_debug = ( std::atoi( argv[ 8 ] ) == 1 ); // Read image - itk::ImageFileReader< TImage >::Pointer input_image_reader = - itk::ImageFileReader< TImage >::New( ); + TImageReader::Pointer input_image_reader = TImageReader::New( ); input_image_reader->SetFileName( input_image_fn ); try { @@ -89,12 +81,11 @@ int main( int argc, char* argv[] ) TImage::ConstPointer input_image = input_image_reader->GetOutput( ); // Show image - itk::ImageToVTKImageFilter< TImage >::Pointer vtk_image = - itk::ImageToVTKImageFilter< TImage >::New( ); + TVTKImage::Pointer vtk_image = TVTKImage::New( ); vtk_image->SetInput( input_image ); vtk_image->Update( ); - fpa::VTK::ImageMPR view; + TMPR view; view.SetBackground( 0.3, 0.2, 0.8 ); view.SetSize( 800, 800 ); view.SetImage( vtk_image->GetOutput( ) ); @@ -114,8 +105,7 @@ int main( int argc, char* argv[] ) input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx ); // Segment input image - fpa::Image::RegionGrowWithMultipleThresholds< TImage >::Pointer segmentor = - fpa::Image::RegionGrowWithMultipleThresholds< TImage >::New( ); + TSegmentor::Pointer segmentor = TSegmentor::New( ); segmentor->AddThresholds( thr_0, thr_1, step ); segmentor->AddSeed( seed_idx, 0 ); segmentor->SetInput( input_image ); @@ -126,9 +116,8 @@ int main( int argc, char* argv[] ) if( visual_debug ) { // Configure observer - fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::Pointer obs = - fpa::VTK::Image3DObserver< fpa::Image::RegionGrowWithMultipleThresholds< TImage >, vtkRenderWindow >::New( ); - obs->SetImage( input_image, view.GetWindow( ) ); + TSegmentorObs::Pointer obs = TSegmentorObs::New( ); + obs->SetRenderWindow( view.GetWindow( ) ); segmentor->AddObserver( itk::AnyEvent( ), obs ); segmentor->ThrowEventsOn( ); } @@ -141,36 +130,40 @@ int main( int argc, char* argv[] ) std::cout << "Segmentation time = " << seconds << std::endl; // Extract distance map - itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::Pointer - distanceMap = - itk::DanielssonDistanceMapImageFilter< TImage, TDistanceMap >::New( ); + 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 - fpa::Image::Dijkstra< TDistanceMap, TScalar >::Pointer paths = - fpa::Image::Dijkstra< TDistanceMap, TScalar >::New( ); + TDijkstra::Pointer paths = TDijkstra::New( ); paths->AddSeed( segmentor->GetSeed( 0 ), TScalar( 0 ) ); paths->SetInput( distanceMap->GetOutput( ) ); - paths->SetNeighborhoodOrder( 1 ); - - // Associate an inversion cost function - fpa::Base::Functors::InvertCostFunction< TScalar >::Pointer - cost_function = - fpa::Base::Functors::InvertCostFunction< TScalar >::New( ); - paths->SetCostConversion( cost_function ); + paths->SetNeighborhoodOrder( 2 ); if( visual_debug ) { // Configure observer - fpa::VTK::Image3DObserver< fpa::Image::Dijkstra< TDistanceMap, TScalar >, vtkRenderWindow >::Pointer obs = - fpa::VTK::Image3DObserver< fpa::Image::Dijkstra< TDistanceMap, TScalar >, vtkRenderWindow >::New( ); - obs->SetImage( distanceMap->GetOutput( ), view.GetWindow( ) ); + TDijkstraObs::Pointer obs = TDijkstraObs::New( ); + obs->SetRenderWindow( view.GetWindow( ) ); paths->AddObserver( itk::AnyEvent( ), obs ); paths->ThrowEventsOn( ); } @@ -182,40 +175,82 @@ int main( int argc, char* argv[] ) seconds = double( end - start ) / double( CLOCKS_PER_SEC ); std::cout << "Paths extraction time = " << seconds << std::endl; - // Show result - /* - TVTKImage::Pointer output_image_vtk = TVTKImage::New( ); - output_image_vtk->SetInput( segmentor->GetOutput( ) ); - output_image_vtk->Update( ); + std::cout << "Final seed: " << paths->GetSeed( 0 ) << std::endl; - vtkSmartPointer< vtkImageMarchingCubes > mc = - vtkSmartPointer< vtkImageMarchingCubes >::New( ); - mc->SetInputData( output_image_vtk->GetOutput( ) ); - mc->SetValue( - 0, - double( segmentor->GetInsideValue( ) ) * double( 0.95 ) - ); - mc->Update( ); - - // Let some interaction and close program - view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 ); - view.Start( ); + // Create polydata + vtkSmartPointer< vtkPoints > points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > cells = + vtkSmartPointer< vtkCellArray >::New( ); + vtkSmartPointer< vtkFloatArray > scalars = + vtkSmartPointer< vtkFloatArray >::New( ); - // Write resulting image - TImageWriter::Pointer output_image_writer = TImageWriter::New( ); - output_image_writer->SetInput( segmentor->GetOutput( ) ); - output_image_writer->SetFileName( output_image_fn ); - try - { - output_image_writer->Update( ); - } - catch( itk::ExceptionObject& err ) + 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 ); + return( -1 ); + + } // yrt - } // yrt - */ return( 0 ); }