X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2Fexamples%2Fexample_Image_Dijkstra_EndPointDetection.cxx;h=a334de3ae881a372578e7e999cc01d9315695d6f;hb=8cf92baeeefc45bc80a2d0eb2850e04b6089111f;hp=47585ad30327975306210ad1eb51fe2a5cd912c8;hpb=8480ad8f9b3e3b556aa8f3f25ff07daf1db6dc78;p=FrontAlgorithms.git diff --git a/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx b/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx index 47585ad..a334de3 100644 --- a/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx +++ b/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx @@ -20,26 +20,16 @@ #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 +#include // ------------------------------------------------------------------------- const unsigned int Dim = 3; @@ -54,6 +44,9 @@ typedef itk::ImageToVTKImageFilter< TImage > TVTKInputImage; 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 @@ -62,16 +55,30 @@ void DistanceMap( // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { - if( argc < 2 ) + if( argc < 9 ) { std::cerr - << "Usage: " << argv[ 0 ] - << " input_image" + << "Usage: " << argv[ 0 ] << std::endl + << " input_image" << 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 ]; + std::string distancemap_fn = argv[ 2 ]; + std::string output_costmap_fn = argv[ 3 ]; + std::string output_labels_fn = argv[ 4 ]; + std::string mst_output_fn = argv[ 5 ]; + std::string endpoints_output_fn = argv[ 6 ]; + std::string bifurcations_output_fn = argv[ 7 ]; + std::string branches_output_fn = argv[ 8 ]; // Read image TImage::Pointer input_image; @@ -88,6 +95,7 @@ int main( int argc, char* argv[] ) return( 1 ); } // yrt + TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( ); vtk_input_image->SetInput( input_image ); vtk_input_image->Update( ); @@ -98,6 +106,13 @@ int main( int argc, char* argv[] ) view.SetSize( 800, 800 ); view.SetImage( vtk_input_image->GetOutput( ) ); + vtkSmartPointer< vtkImageMarchingCubes > mc = + vtkSmartPointer< vtkImageMarchingCubes >::New( ); + mc->SetInputData( vtk_input_image->GetOutput( ) ); + mc->SetValue( 0, 1e-1 ); + mc->Update( ); + view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 ); + // Allow some interaction and wait for, at least, one seed view.Render( ); while( view.GetNumberOfSeeds( ) == 0 ) @@ -148,25 +163,183 @@ int main( int argc, char* argv[] ) << std::difftime( end, start ) << " s." << std::endl; - /* TODO - // Save final total cost map - itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer = - itk::ImageFileWriter< TOutputImage >::New( ); - output_image_writer->SetFileName( output_image_fn ); - output_image_writer->SetInput( filter->GetOutput( ) ); - try + // 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( ); + + // Show endpoints + vtkSmartPointer< vtkPoints > endpoints_points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > endpoints_cells = + vtkSmartPointer< vtkCellArray >::New( ); + for( + TFilter::TUniqueVertices::ConstIterator epIt = endpoints->Begin( ); + epIt != endpoints->End( ); + ++epIt + ) { - output_image_writer->Update( ); - } - catch( itk::ExceptionObject& err ) + TImage::PointType pnt; + input_image->TransformIndexToPhysicalPoint( *epIt, pnt ); + endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + endpoints_cells->InsertNextCell( 1 ); + endpoints_cells-> + InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 ); + + } // rof + vtkSmartPointer< vtkPolyData > endpoints_polydata = + vtkSmartPointer< vtkPolyData >::New( ); + endpoints_polydata->SetPoints( endpoints_points ); + endpoints_polydata->SetVerts( endpoints_cells ); + view.AddPolyData( endpoints_polydata, 0, 0, 1, 1 ); + + // Show bifurcations + vtkSmartPointer< vtkPoints > bifurcations_points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > bifurcations_cells = + vtkSmartPointer< vtkCellArray >::New( ); + for( + TFilter::TUniqueVertices::ConstIterator bfIt = bifurcations->Begin( ); + bfIt != bifurcations->End( ); + ++bfIt + ) { - std::cerr - << "Error while writing image to " << output_image_fn << ": " - << err << std::endl; - return( 1 ); + TImage::PointType pnt; + input_image->TransformIndexToPhysicalPoint( *bfIt, pnt ); + bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + bifurcations_cells->InsertNextCell( 1 ); + bifurcations_cells-> + InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 ); + + } // rof + vtkSmartPointer< vtkPolyData > bifurcations_polydata = + vtkSmartPointer< vtkPolyData >::New( ); + bifurcations_polydata->SetPoints( bifurcations_points ); + bifurcations_polydata->SetVerts( bifurcations_cells ); + view.AddPolyData( bifurcations_polydata, 0, 1, 0, 1 ); + + // Show branches (simple and detailed) + vtkSmartPointer< vtkPoints > simple_branches_points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > simple_branches_cells = + vtkSmartPointer< vtkCellArray >::New( ); + + vtkSmartPointer< vtkPoints > detailed_branches_points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > detailed_branches_cells = + vtkSmartPointer< vtkCellArray >::New( ); + vtkSmartPointer< vtkFloatArray > detailed_branches_scalars = + vtkSmartPointer< vtkFloatArray >::New( ); + + TFilter::TBranches::ConstIterator brIt = branches->Begin( ); + for( ; brIt != branches->End( ); ++brIt ) + { + // Branch's first point + TImage::PointType first_point; + input_image->TransformIndexToPhysicalPoint( brIt->first, first_point ); + unsigned long first_id = simple_branches_points->GetNumberOfPoints( ); + simple_branches_points->InsertNextPoint( + first_point[ 0 ], first_point[ 1 ], first_point[ 2 ] + ); + + TFilter::TBranches::ConstRowIterator brRowIt = branches->Begin( brIt ); + for( ; brRowIt != branches->End( brIt ); ++brRowIt ) + { + // Branch's second point + TImage::PointType second_point; + input_image-> + TransformIndexToPhysicalPoint( brRowIt->first, second_point ); + unsigned long second_id = simple_branches_points->GetNumberOfPoints( ); + simple_branches_points->InsertNextPoint( + second_point[ 0 ], second_point[ 1 ], second_point[ 2 ] + ); + simple_branches_cells->InsertNextCell( 2 ); + simple_branches_cells->InsertCellPoint( first_id ); + simple_branches_cells->InsertCellPoint( second_id ); + + // Detailed path + double pathId = double( brRowIt->second - 1 ) / double( nBranches - 1 ); + TFilter::TVertices path; + mst->GetPath( path, brIt->first, brRowIt->first ); + TFilter::TVertices::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ); ++pIt ) + { + TImage::PointType path_point; + input_image->TransformIndexToPhysicalPoint( *pIt, path_point ); + detailed_branches_points->InsertNextPoint( + path_point[ 0 ], path_point[ 1 ], path_point[ 2 ] + ); + detailed_branches_scalars->InsertNextTuple1( pathId ); + if( pIt != path.begin( ) ) + { + unsigned long nPoints = + detailed_branches_points->GetNumberOfPoints( ); + detailed_branches_cells->InsertNextCell( 2 ); + detailed_branches_cells->InsertCellPoint( nPoints - 2 ); + detailed_branches_cells->InsertCellPoint( nPoints - 1 ); + + } // fi + + } // rof + + } // rof + + } // rof + vtkSmartPointer< vtkPolyData > simple_branches_polydata = + vtkSmartPointer< vtkPolyData >::New( ); + simple_branches_polydata->SetPoints( simple_branches_points ); + simple_branches_polydata->SetLines( simple_branches_cells ); + view.AddPolyData( simple_branches_polydata, 1, 0, 1, 1 ); + + vtkSmartPointer< vtkPolyData > detailed_branches_polydata = + vtkSmartPointer< vtkPolyData >::New( ); + detailed_branches_polydata->SetPoints( detailed_branches_points ); + detailed_branches_polydata->SetLines( detailed_branches_cells ); + detailed_branches_polydata-> + GetPointData( )->SetScalars( detailed_branches_scalars ); + view.AddPolyData( detailed_branches_polydata, 1 ); + + // Let some more interaction + view.Render( ); + view.Start( ); - } // yrt - */ return( 0 ); } @@ -197,6 +370,27 @@ void ReadImage( typename I::Pointer& image, const std::string& filename ) 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(