X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2Fexamples%2Fexample_ImageAlgorithm_Skeletonization.cxx;h=22b85e773ec9b5ead755e23d76f51a3c8b44fe04;hb=67b67cafced0d039cf6ff2ccf7839088fd091395;hp=e0121bd6ba9c173d41a532a56e6d0cb0c507e1b8;hpb=97940c7ac873a39428e8739b2d47ca8485cff70e;p=FrontAlgorithms.git diff --git a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx index e0121bd..22b85e7 100644 --- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx +++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx @@ -9,17 +9,21 @@ #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; @@ -28,19 +32,23 @@ 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 < 7 ) + if( argc < 8 ) { std::cerr << "Usage: " << argv[ 0 ] << " input_image thr_0 thr_1 step" - << " output_segmentation output_distancemap" + << " output_segmentation output_distancemap output_dijkstra" << " visual_debug" << std::endl; return( 1 ); @@ -52,9 +60,10 @@ int main( int argc, char* argv[] ) 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 > 7 ) - visual_debug = ( std::atoi( argv[ 7 ] ) == 1 ); + if( argc > 8 ) + visual_debug = ( std::atoi( argv[ 8 ] ) == 1 ); // Read image TImageReader::Pointer input_image_reader = TImageReader::New( ); @@ -125,60 +134,123 @@ int main( int argc, char* argv[] ) 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; - std::cout << "Final seed: " << seed_idx << 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 ); - TDistanceMapWriter::Pointer distancemap_writer = - TDistanceMapWriter::New( ); - distancemap_writer->SetInput( distanceMap->GetOutput( ) ); - distancemap_writer->SetFileName( output_distancemap_fn ); - distancemap_writer->Update( ); + } // rof - TImageWriter::Pointer segmentation_writer = - TImageWriter::New( ); + 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 ); - segmentation_writer->Update( ); - // Show result - /* - TVTKImage::Pointer output_image_vtk = TVTKImage::New( ); - output_image_vtk->SetInput( segmentor->GetOutput( ) ); - output_image_vtk->Update( ); + itk::ImageFileWriter< TDistanceMap >::Pointer dmap_writer = + itk::ImageFileWriter< TDistanceMap >::New( ); + dmap_writer->SetInput( distanceMap->GetOutput( ) ); + dmap_writer->SetFileName( output_distancemap_fn ); - 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( ); + itk::ImageFileWriter< TDistanceMap >::Pointer dijk_writer = + itk::ImageFileWriter< TDistanceMap >::New( ); + dijk_writer->SetInput( paths->GetOutput( ) ); + dijk_writer->SetFileName( output_dijkstra_fn ); - // 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 ) - { + 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 ); }