X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2Fexamples%2Fexample_ImageAlgorithmDijkstra_03.cxx;h=9f529f1d3966cb94a0f1865638fb71f8cab77f04;hb=8480ad8f9b3e3b556aa8f3f25ff07daf1db6dc78;hp=b78929b5685ba0c6d85e214751ee745700c9c9c6;hpb=b63dc485b7255d1ab70ff72096beafe13a71f1be;p=FrontAlgorithms.git diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx index b78929b..9f529f1 100644 --- a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx +++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx @@ -1,7 +1,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -14,11 +16,17 @@ #include #include +#include +#include + #include #include #include #include +#include #include +#include +#include #include #include @@ -26,40 +34,55 @@ // ------------------------------------------------------------------------- const unsigned int Dim = 3; -typedef double TPixel; -typedef double TScalar; +typedef float TPixel; +typedef float TScalar; typedef itk::Image< TPixel, Dim > TImage; typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; typedef itk::ImageFileReader< TImage > TImageReader; -typedef itk::ImageFileWriter< TImage > TImageWriter; +typedef itk::ImageFileWriter< TImage > TImageWriter; typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra; typedef fpa::VTK::ImageMPR TMPR; typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs; +struct TBranch +{ + double Length; + TImage::PointType::VectorType V1; + TImage::PointType::VectorType V2; + + bool operator<( const TBranch& other ) const + { + return( other.Length < this->Length ); + } +}; +typedef std::multiset< TBranch > TBranches; + // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { - if( argc < 6 ) + if( argc < 8 ) { std::cerr << "Usage: " << argv[ 0 ] - << " input_image x y z output_image" + << " input_image x y z output_image output_imagej output_vtk" << " visual_debug" << std::endl; return( 1 ); } // fi std::string input_image_fn = argv[ 1 ]; - TImage::IndexType seed_idx; - seed_idx[ 0 ] = std::atoi( argv[ 2 ] ); - seed_idx[ 1 ] = std::atoi( argv[ 3 ] ); - seed_idx[ 2 ] = std::atoi( argv[ 4 ] ); + TImage::PointType seed_pnt; + seed_pnt[ 0 ] = std::atof( argv[ 2 ] ); + seed_pnt[ 1 ] = std::atof( argv[ 3 ] ); + seed_pnt[ 2 ] = std::atof( argv[ 4 ] ); std::string output_image_fn = argv[ 5 ]; + std::string output_imagej_fn = argv[ 6 ]; + std::string output_vtk_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 TImageReader::Pointer input_image_reader = TImageReader::New( ); @@ -86,14 +109,34 @@ int main( int argc, char* argv[] ) view.SetSize( 800, 800 ); view.SetImage( vtk_image->GetOutput( ) ); + /* TODO + vtkSmartPointer< vtkImageMarchingCubes > mc = + vtkSmartPointer< vtkImageMarchingCubes >::New( ); + mc->SetInputData( vtk_image->GetOutput( ) ); + mc->SetValue( 0, 1e-2 ); + mc->Update( ); + view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 ); + */ + // Allow some interaction - view.Start( ); + view.Render( ); + if( visual_debug ) + view.Start( ); + + TImage::IndexType seed_idx; + input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx ); + /* TODO + std::cout << seed_idx << " " << seed_pnt << std::endl; + seed_idx[ 0 ] = 256; + seed_idx[ 1 ] = 313; + seed_idx[ 2 ] = 381; + */ // Extract paths TDijkstra::Pointer paths = TDijkstra::New( ); paths->AddSeed( seed_idx, TScalar( 0 ) ); paths->SetInput( input_image ); - paths->SetNeighborhoodOrder( 1 ); + paths->SetNeighborhoodOrder( 2 ); if( visual_debug ) { @@ -119,33 +162,112 @@ int main( int argc, char* argv[] ) vtkSmartPointer< vtkFloatArray > scalars = vtkSmartPointer< vtkFloatArray >::New( ); + TDijkstra::TMarkImage::Pointer new_marks = TDijkstra::TMarkImage::New( ); + new_marks->SetLargestPossibleRegion( input_image->GetLargestPossibleRegion( ) ); + new_marks->SetRequestedRegion( input_image->GetRequestedRegion( ) ); + new_marks->SetBufferedRegion( input_image->GetBufferedRegion( ) ); + new_marks->SetOrigin( input_image->GetOrigin( ) ); + new_marks->SetSpacing( input_image->GetSpacing( ) ); + new_marks->SetDirection( input_image->GetDirection( ) ); + new_marks->Allocate( ); + new_marks->FillBuffer( 0 ); + + const TDijkstra::TMarkImage* marks = paths->GetOutputMarkImage( ); + TDijkstra::TMark max_mark = paths->GetNumberOfBranches( ); const TDijkstra::TVertices& endpoints = paths->GetEndPoints( ); - const TDijkstra::TTree& tree = paths->GetFinalTree( ); + const TDijkstra::TVertices& bifurcations = paths->GetBifurcationPoints( ); + const TDijkstra::TTree& tree = paths->GetFullTree( ); + const TDijkstra::TTree& reduced_tree = paths->GetReducedTree( ); TDijkstra::TVertices::const_iterator epIt = endpoints.begin( ); - for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId ) + + TBranches branches; + TImage::PointType ori = input_image->GetOrigin( ); + + TDijkstra::TTree::const_iterator rIt = reduced_tree.begin( ); + for( ; rIt != reduced_tree.end( ); ++rIt ) { - double pd = double( epId ) / double( endpoints.size( ) - 1 ); + TDijkstra::TTree::const_iterator tIt = tree.find( rIt->first ); - TDijkstra::TVertex idx = *epIt; + // Prepare branch "a la ImageJ" + TBranch branch; + TImage::PointType p1, p2; + input_image->TransformIndexToPhysicalPoint( rIt->second.first, p1 ); + input_image->TransformIndexToPhysicalPoint( rIt->first, p2 ); + branch.V1 = p1 - ori; + branch.V2 = p2 - ori; + branch.Length = double( 0 ); + + double pd = + ( double( tIt->second.second ) - double( 1 ) ) / + ( double( max_mark ) - double( 1 ) ); + bool start = true; + TImage::PointType prev; do { + TDijkstra::TVertex idx = tIt->first; + new_marks->SetPixel( idx, 255 ); TImage::PointType pnt; input_image->TransformIndexToPhysicalPoint( idx, pnt ); - points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); scalars->InsertNextTuple1( pd ); - if( idx != *epIt ) + if( !start ) { cells->InsertNextCell( 2 ); cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + branch.Length += prev.EuclideanDistanceTo( pnt ); } // fi - idx = tree.find( idx )->second; + start = false; + prev = pnt; + tIt = tree.find( tIt->second.first ); + + } while( tIt->first != rIt->second.first ); + + // Last point + TDijkstra::TVertex idx = tIt->first; + new_marks->SetPixel( idx, 255 ); + TImage::PointType pnt; + input_image->TransformIndexToPhysicalPoint( idx, pnt ); + points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + scalars->InsertNextTuple1( pd ); + if( !start ) + { + cells->InsertNextCell( 2 ); + cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); + cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + branch.Length += prev.EuclideanDistanceTo( pnt ); + + } // fi + branches.insert( branch ); + + } // rof + + std::ofstream bout( output_imagej_fn.c_str( ) ); + if( !bout ) + { + std::cerr << "Error: could not open " << output_imagej_fn << std::endl; + return( 1 ); - } while( idx != tree.find( idx )->second ); + } // fi + int i = 1; + for( TBranches::const_iterator bIt = branches.begin( ); bIt != branches.end( ); ++bIt, ++i ) + { + bout + << std::setiosflags( std::ios::fixed) << std::setprecision( 3 ) + << i << "\t1.000\t" + << bIt->Length << "\t" + << bIt->V1[ 0 ] << "\t" + << bIt->V1[ 1 ] << "\t" + << bIt->V1[ 2 ] << "\t" + << bIt->V2[ 0 ] << "\t" + << bIt->V2[ 1 ] << "\t" + << bIt->V2[ 2 ] << "\t" + << ( bIt->V2 - bIt->V1 ).GetNorm( ) + << std::endl; } // rof + bout.close( ); vtkSmartPointer< vtkPolyData > vtk_tree = vtkSmartPointer< vtkPolyData >::New( ); @@ -155,15 +277,26 @@ int main( int argc, char* argv[] ) view.AddPolyData( vtk_tree ); view.Render( ); - view.Start( ); + if( visual_debug ) + view.Start( ); - /* TODO - TImageWriter::Pointer dijkstra_writer = - TImageWriter::New( ); - dijkstra_writer->SetInput( paths->GetOutput( ) ); - dijkstra_writer->SetFileName( "dijkstra.mhd" ); - dijkstra_writer->Update( ); - */ + vtkSmartPointer< vtkPolyDataWriter > writer = + vtkSmartPointer< vtkPolyDataWriter >::New( ); + writer->SetInputData( vtk_tree ); + writer->SetFileName( output_vtk_fn.c_str( ) ); + writer->Update( ); + + itk::ImageFileWriter< TDijkstra::TMarkImage >::Pointer marks_writer = + itk::ImageFileWriter< TDijkstra::TMarkImage >::New( ); + marks_writer->SetInput( new_marks ); + marks_writer->SetFileName( output_image_fn ); + marks_writer->Update( ); + + TImageWriter::Pointer dijkstra_writer = + TImageWriter::New( ); + dijkstra_writer->SetInput( paths->GetOutput( ) ); + dijkstra_writer->SetFileName( "dijkstra.mhd" ); + dijkstra_writer->Update( ); return( 0 ); }