X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2Fexamples%2Fexample_ImageAlgorithmDijkstra_03.cxx;h=9f529f1d3966cb94a0f1865638fb71f8cab77f04;hb=67b67cafced0d039cf6ff2ccf7839088fd091395;hp=b3c04eb033305d165caf2f62a0e98e05427ca5cb;hpb=97940c7ac873a39428e8739b2d47ca8485cff70e;p=FrontAlgorithms.git diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx index b3c04eb..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,375 +16,73 @@ #include #include +#include +#include + #include #include +#include #include +#include #include +#include +#include -#include -#include -#include +#include #include #include // ------------------------------------------------------------------------- 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 fpa::Image::Dijkstra< TImage, TScalar > TBaseDijkstra; +typedef itk::ImageFileWriter< TImage > TImageWriter; +typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra; typedef fpa::VTK::ImageMPR TMPR; -typedef fpa::VTK::Image3DObserver< TBaseDijkstra, vtkRenderWindow > TDijkstraObs; +typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs; -// ------------------------------------------------------------------------- -class TDijkstra - : public TBaseDijkstra +struct TBranch { -public: - typedef TDijkstra Self; - typedef TBaseDijkstra Superclass; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - - typedef Superclass::TCost TCost; - typedef Superclass::TVertex TVertex; - typedef Superclass::InputImageType TImage; - typedef std::deque< TVertex > TVertices; - - typedef Superclass::TEndEvent TEndEvent; - typedef Superclass::TBacktrackingEvent TBacktrackingEvent; - -protected: - typedef std::pair< TCost, TVertex > _TCandidate; - typedef std::multimap< TCost, TVertex > _TCandidates; - typedef Superclass::_TNode _TNode; - -public: - itkNewMacro( Self ); - itkTypeMacro( TDijkstra, TBaseDijkstra ); - -protected: - TDijkstra( ) - : Superclass( ) - { } - virtual ~TDijkstra( ) - { } - - virtual void _BeforeMainLoop( ) - { - const TImage* img = this->GetInput( ); - - // Correct seeds - TImage::SizeType radius; - radius.Fill( 3 ); - itk::ConstNeighborhoodIterator< TImage > it( - radius, img, img->GetRequestedRegion( ) - ); - for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s ) - { - _TNode seed = this->m_Seeds[ s ]; - it.SetLocation( seed.Vertex ); - - TImage::SizeType size = it.GetSize( ); - unsigned int nN = 1; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - nN *= size[ d ]; - TImage::PixelType max_value = img->GetPixel( seed.Vertex ); - for( unsigned int i = 0; i < nN; ++i ) - { - if( it.GetPixel( i ) > max_value ) - { - seed.Vertex = it.GetIndex( i ); - seed.Parent = seed.Vertex; - max_value = it.GetPixel( i ); - - } // fi - - } // rof - this->m_Seeds[ s ] = seed; - - } // rof - - // End initialization - this->Superclass::_BeforeMainLoop( ); - this->m_Candidates.clear( ); - } + double Length; + TImage::PointType::VectorType V1; + TImage::PointType::VectorType V2; - virtual void _AfterMainLoop( ) + bool operator<( const TBranch& other ) const { - typedef unsigned char _TMark; - typedef itk::Image< _TMark, TImage::ImageDimension > _TMarkImage; - - this->Superclass::_AfterMainLoop( ); - if( this->m_Candidates.size( ) == 0 ) - return; - - this->InvokeEvent( TEndEvent( ) ); - - const TImage* input = this->GetInput( ); - TImage::SpacingType spacing = input->GetSpacing( ); - - // Prepare an object to hold marks - _TMarkImage::Pointer marks = _TMarkImage::New( ); - marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) ); - marks->SetRequestedRegion( input->GetRequestedRegion( ) ); - marks->SetBufferedRegion( input->GetBufferedRegion( ) ); - marks->SetOrigin( input->GetOrigin( ) ); - marks->SetSpacing( spacing ); - marks->SetDirection( input->GetDirection( ) ); - marks->Allocate( ); - marks->FillBuffer( _TMark( 0 ) ); - - // Iterate over the candidates, starting from the candidates that - // are near thin branches - _TCandidates::const_reverse_iterator cIt = - this->m_Candidates.rbegin( ); - for( int backId = 0; backId < 100 && cIt != this->m_Candidates.rend( ); ++cIt ) - { - // If pixel hasn't been visited, continue - TVertex v = cIt->second; - if( marks->GetPixel( v ) != _TMark( 0 ) ) - continue; - - // Compute nearest start candidate - TImage::SizeType radius; - radius.Fill( 3 ); - itk::ConstNeighborhoodIterator< TImage > iIt( - radius, input, input->GetRequestedRegion( ) - ); - iIt.SetLocation( v ); - unsigned int nN = 1; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - nN *= iIt.GetSize( )[ d ]; - TVertex max_vertex = iIt.GetIndex( 0 ); - TImage::PixelType max_value = iIt.GetPixel( 0 ); - for( unsigned int i = 1; i < nN; ++i ) - { - TImage::PixelType value = iIt.GetPixel( i ); - if( value > max_value ) - { - max_value = value; - max_vertex = iIt.GetIndex( i ); - - } // fi - - } // rof - - if( marks->GetPixel( max_vertex ) != _TMark( 0 ) ) - continue; - backId++; - std::cout << "Leaf: " << backId << " " << max_value << " " << max_vertex << std::endl; - - bool start = true; - do - { - double dist = double( 1.5 ) * std::sqrt( double( input->GetPixel( max_vertex ) ) ); - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - radius[ d ] = - ( unsigned long )( dist / spacing[ d ] ) + 1; - itk::NeighborhoodIterator< _TMarkImage > mIt( - radius, marks, marks->GetRequestedRegion( ) - ); - mIt.SetLocation( max_vertex ); - nN = 1; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - nN *= mIt.GetSize( )[ d ]; - for( unsigned int i = 0; i < nN; ++i ) - if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) ) - { - mIt.SetPixel( i, ( start )? 255: 100 ); - start = false; - } - - /* - TImage::SizeType radius; - mIt.GoToBegin( ); - mIt.SetLocation( vIt ); - - TImage::SizeType size = mIt.GetSize( ); - unsigned int nN = 1; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - nN *= size[ d ]; - for( unsigned int i = 0; i < nN; ++i ) - if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) ) - mIt.SetPixel( i, ( start )? 255: 100 ); - - start = false; - */ - // Next vertex in current path - this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) ); - max_vertex = this->_Parent( max_vertex ); - - } while( max_vertex != this->_Parent( max_vertex ) ); - - } // rof - /* - else - marks->SetPixel( v, _TMark( 1 ) ); - } // rof - */ - - /* - itk::ImageFileWriter< _TMarkImage >::Pointer w = - itk::ImageFileWriter< _TMarkImage >::New( ); - w->SetInput( marks ); - w->SetFileName( "marks.mhd" ); - w->Update( ); - - - this->m_Axes = vtkSmartPointer< vtkPolyData >::New( ); - vtkSmartPointer< vtkPoints > points = - vtkSmartPointer< vtkPoints >::New( ); - vtkSmartPointer< vtkCellArray > lines = - vtkSmartPointer< vtkCellArray >::New( ); - - { - - backId++; - std::cout << "Leaf: " << backId << " " << cIt->first << " " << vIt << std::endl; - bool start = true; - do - { - double dist = double( input->GetPixel( vIt ) ); - TImage::SizeType radius; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - radius[ d ] = - ( unsigned long )( dist / spacing[ d ] ) + 1; - itk::NeighborhoodIterator< _TMarkImage > mIt( - radius, marks, marks->GetRequestedRegion( ) - ); - mIt.GoToBegin( ); - mIt.SetLocation( vIt ); - - TImage::SizeType size = mIt.GetSize( ); - unsigned int nN = 1; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - nN *= size[ d ]; - for( unsigned int i = 0; i < nN; ++i ) - if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) ) - mIt.SetPixel( i, ( start )? 255: 100 ); - - start = false; - - // Next vertex - vIt = this->_Parent( vIt ); - - } while( vIt != this->_Parent( vIt ) ); - - // Last vertex - // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl; - */ - /* TODO - TVertices pS; - TVertex parent = this->_Parent( vIt ); - while( parent != vIt ) - { - pS.push_front( vIt ); - vIt = parent; - parent = this->_Parent( vIt ); - - TImage::PointType pnt; - this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt ); - points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - if( points->GetNumberOfPoints( ) > 1 ) - { - lines->InsertNextCell( 2 ); - lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); - lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); - - } // fi - - } // elihw - pS.push_front( vIt ); - TImage::PointType pnt; - this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt ); - points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - if( points->GetNumberOfPoints( ) > 1 ) - { - lines->InsertNextCell( 2 ); - lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); - lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); - - } // fi - - this->m_Axes->SetPoints( points ); - this->m_Axes->SetLines( lines ); - */ - - /* - } // rof - */ + return( other.Length < this->Length ); } - - virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n ) - { - TCost nc = this->_Cost( nn.Vertex, n.Vertex ); - if( TCost( 0 ) < nc ) - { - nn.Cost = n.Cost + ( TCost( 1 ) / nc ); - nn.Result = nn.Cost; - return( true ); - } - else - return( false ); - } - - virtual bool _UpdateResult( _TNode& n ) - { - bool ret = this->Superclass::_UpdateResult( n ); - if( ret ) - { - TCost nc = this->_Cost( n.Vertex, n.Parent ); - if( TCost( 0 ) < nc ) - { - TCost cc = n.Cost / nc; - this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) ); - /* TODO: debug code - this->GetOutput( )->SetPixel( n.Vertex, cc ); - */ - } // fi - - } // fi - return( ret ); - } - -private: - TDijkstra( const Self& other ); - Self& operator=( const Self& other ); - -protected: - _TCandidates m_Candidates; -public: - vtkSmartPointer< vtkPolyData > m_Axes; }; +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( ); @@ -409,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 ) { @@ -434,8 +154,143 @@ int main( int argc, char* argv[] ) double seconds = double( end - start ) / double( CLOCKS_PER_SEC ); std::cout << "Paths extraction time = " << seconds << std::endl; - view.AddPolyData( paths->m_Axes, 1, 0, 0 ); - view.Start( ); + // Create polydata + vtkSmartPointer< vtkPoints > points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > cells = + vtkSmartPointer< vtkCellArray >::New( ); + 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::TVertices& bifurcations = paths->GetBifurcationPoints( ); + const TDijkstra::TTree& tree = paths->GetFullTree( ); + const TDijkstra::TTree& reduced_tree = paths->GetReducedTree( ); + TDijkstra::TVertices::const_iterator epIt = endpoints.begin( ); + + TBranches branches; + TImage::PointType ori = input_image->GetOrigin( ); + + TDijkstra::TTree::const_iterator rIt = reduced_tree.begin( ); + for( ; rIt != reduced_tree.end( ); ++rIt ) + { + TDijkstra::TTree::const_iterator tIt = tree.find( rIt->first ); + + // 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( !start ) + { + cells->InsertNextCell( 2 ); + cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); + cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + branch.Length += prev.EuclideanDistanceTo( pnt ); + + } // fi + 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 ); + + } // 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( ); + vtk_tree->SetPoints( points ); + vtk_tree->SetLines( cells ); + vtk_tree->GetPointData( )->SetScalars( scalars ); + + view.AddPolyData( vtk_tree ); + view.Render( ); + if( visual_debug ) + view.Start( ); + + 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( ); @@ -443,40 +298,6 @@ int main( int argc, char* argv[] ) dijkstra_writer->SetFileName( "dijkstra.mhd" ); dijkstra_writer->Update( ); - // Show result - /* - TVTKImage::Pointer output_image_vtk = TVTKImage::New( ); - output_image_vtk->SetInput( segmentor->GetOutput( ) ); - output_image_vtk->Update( ); - - 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( ); - - // 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 ) - { - std::cerr << "Error caught: " << err << std::endl; - return( 1 ); - - } // yrt - */ return( 0 ); }