From f28d460c6d9ce93e55c073bd04dba8de16d55d3a Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Thu, 16 Apr 2015 12:35:03 -0500 Subject: [PATCH] Skeleton reconstruction added --- appli/examples/CMakeLists.txt | 1 + .../example_Image_Dijkstra_LabelSkeleton.cxx | 421 ++++++++++++++++++ lib/fpa/Image/DijkstraWithEndPointDetection.h | 19 +- .../Image/DijkstraWithEndPointDetection.hxx | 68 +-- .../Image/DijkstraWithSphereBacktracking.h | 97 ---- .../Image/DijkstraWithSphereBacktracking.hxx | 402 ----------------- 6 files changed, 478 insertions(+), 530 deletions(-) create mode 100644 appli/examples/example_Image_Dijkstra_LabelSkeleton.cxx delete mode 100644 lib/fpa/Image/DijkstraWithSphereBacktracking.h delete mode 100644 lib/fpa/Image/DijkstraWithSphereBacktracking.hxx diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt index cf7ebf4..122dcee 100644 --- a/appli/examples/CMakeLists.txt +++ b/appli/examples/CMakeLists.txt @@ -17,6 +17,7 @@ SET( example_Image_Dijkstra_DanielssonCost example_Image_Dijkstra_DanielssonCost_TwoSeedsPath example_Image_Dijkstra_EndPointDetection + example_Image_Dijkstra_LabelSkeleton example_ShowSkeleton ) FOREACH(EX ${SIMPLE_VTK_EXAMPLES}) diff --git a/appli/examples/example_Image_Dijkstra_LabelSkeleton.cxx b/appli/examples/example_Image_Dijkstra_LabelSkeleton.cxx new file mode 100644 index 0000000..4e65d9d --- /dev/null +++ b/appli/examples/example_Image_Dijkstra_LabelSkeleton.cxx @@ -0,0 +1,421 @@ +#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 unsigned char TPixel; +typedef float TScalar; + +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::Image< TScalar, Dim > TScalarImage; +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 F > +class SkeletonCostFunction + : public fpa::Image::Functors::ImageCostFunction< typename F::TInputImage, typename F::TResult > +{ +public: + typedef fpa::Image::Functors::ImageCostFunction< typename F::TInputImage, typename F::TResult > Superclass; + typedef SkeletonCostFunction Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TInputImage TInputImage; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TIndex TIndex; + +public: + itkNewMacro( Self ); + itkTypeMacro( SkeletonCostFunction, fpa_Image_Functors_ImageCostFunction ); + +public: + virtual void SetInputImage( const TInputImage* img ) + { + this->Superclass::SetInputImage( img ); + + typename TInputImage::SpacingType spac = img->GetSpacing( ); + typename TInputImage::SpacingType::ValueType min_spac = spac[ 0 ]; + for( unsigned int d = 1; d < TInputImage::ImageDimension; ++d ) + min_spac = ( spac[ d ] < min_spac )? spac[ d ]: min_spac; + this->m_Distance = TResult( min_spac ); + } + + virtual TResult Evaluate( const TIndex& v, const TIndex& p ) const + { + TResult res = this->Superclass::Evaluate( v, p ); + if( res > TResult( 0 ) ) + return( this->m_Distance ); + else + return( TResult( -1 ) ); + } + +protected: + SkeletonCostFunction( ) + : Superclass( ) + { } + virtual ~SkeletonCostFunction( ) + { } + +private: + // Purposely not implemented + SkeletonCostFunction( const Self& ); + void operator=( const Self& ); + +protected: + TResult m_Distance; +}; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 11 ) + { + std::cerr + << "Usage: " << argv[ 0 ] << std::endl + << " input_image" << std::endl + << " seed_x seed_y seed_z" << 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 output_costmap_fn = argv[ 5 ]; + std::string output_labels_fn = argv[ 6 ]; + std::string mst_output_fn = argv[ 7 ]; + std::string endpoints_output_fn = argv[ 8 ]; + std::string bifurcations_output_fn = argv[ 9 ]; + std::string branches_output_fn = argv[ 10 ]; + + // Get seed + TImage::PointType pnt; + pnt[ 0 ] = TImage::PointType::ValueType( std::atof( argv[ 2 ] ) ); + pnt[ 1 ] = TImage::PointType::ValueType( std::atof( argv[ 3 ] ) ); + pnt[ 2 ] = TImage::PointType::ValueType( std::atof( argv[ 4 ] ) ); + + // Read image + TImage::Pointer input_image; + try + { + ReadImage< TImage >( input_image, input_image_fn ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr + << "Error caught while reading \"" + << input_image_fn << "\": " << err + << std::endl; + return( 1 ); + + } // yrt + + TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( ); + vtk_input_image->SetInput( input_image ); + vtk_input_image->Update( ); + + // Show input image and let some interaction + fpa::VTK::ImageMPR view; + view.SetBackground( 0.3, 0.2, 0.8 ); + view.SetSize( 500, 500 ); + 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( ); + view.Start( ); + + // Transform seed + TImage::IndexType seed; + input_image->TransformPhysicalPointToIndex( pnt, seed ); + + // Prepare Dijkstra filter + typedef fpa::Image::DijkstraWithEndPointDetection< TImage, TScalarImage > TFilter; + typedef SkeletonCostFunction< TFilter > TFunction; + + TFunction::Pointer function = TFunction::New( ); + + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( input_image ); + filter->SetCostFunction( function ); + filter->SetNeighborhoodOrder( 2 ); + filter->SetSafetyNeighborhoodSize( 0 ); + filter->StopAtOneFrontOff( ); + filter->CorrectSeedsOff( ); + filter->CorrectEndPointsOff( ); + filter->AddSeed( seed, TScalar( 0 ) ); + + // Prepare graphical debugger + typedef fpa::VTK::Image3DObserver< TFilter, vtkRenderWindow > TDebugger; + TDebugger::Pointer debugger = TDebugger::New( ); + debugger->SetRenderWindow( view.GetWindow( ) ); + debugger->SetRenderPercentage( 0.0000001 ); + filter->AddObserver( itk::AnyEvent( ), debugger ); + filter->ThrowEventsOn( ); + + // Go! + std::time_t start, end; + std::time( &start ); + filter->Update( ); + std::time( &end ); + std::cout + << "Extraction time = " + << std::difftime( end, start ) + << " s." << std::endl; + + // 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( 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 + ) + { + 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 + ) + { + 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( ); + + return( 0 ); +} + +// ------------------------------------------------------------------------- +template< class I > +void ReadImage( typename I::Pointer& image, const std::string& filename ) +{ + typename itk::ImageFileReader< I >::Pointer reader = + itk::ImageFileReader< I >::New( ); + reader->SetFileName( filename ); + reader->Update( ); + image = reader->GetOutput( ); + 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 +} + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.h b/lib/fpa/Image/DijkstraWithEndPointDetection.h index 4288da3..b32d68a 100644 --- a/lib/fpa/Image/DijkstraWithEndPointDetection.h +++ b/lib/fpa/Image/DijkstraWithEndPointDetection.h @@ -76,8 +76,19 @@ namespace fpa itkNewMacro( Self ); itkTypeMacro( DijkstraWithEndPointDetection, Dijkstra ); + itkBooleanMacro( CorrectSeeds ); + itkBooleanMacro( CorrectEndPoints ); + + itkGetConstMacro( CorrectSeeds, bool ); + itkGetConstMacro( CorrectEndPoints, bool ); + itkGetConstMacro( SafetyNeighborhoodSize, unsigned int ); + itkGetConstMacro( NumberOfBranches, unsigned long ); + itkSetMacro( CorrectSeeds, bool ); + itkSetMacro( CorrectEndPoints, bool ); + itkSetMacro( SafetyNeighborhoodSize, unsigned int ); + public: TLabelImage* GetLabelImage( ); const TLabelImage* GetLabelImage( ) const; @@ -125,8 +136,12 @@ namespace fpa unsigned int m_EndPointsIndex; unsigned int m_BranchesIndex; - _TCandidates m_Candidates; - unsigned long m_NumberOfBranches; + bool m_CorrectSeeds; + bool m_CorrectEndPoints; + unsigned int m_SafetyNeighborhoodSize; + + _TCandidates m_Candidates; + unsigned long m_NumberOfBranches; }; } // ecapseman diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx index 76695c2..dc51a64 100644 --- a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx +++ b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx @@ -170,7 +170,10 @@ GraftBranches( itk::DataObject* obj ) template< class I, class O > fpa::Image::DijkstraWithEndPointDetection< I, O >:: DijkstraWithEndPointDetection( ) - : Superclass( ) + : Superclass( ), + m_CorrectSeeds( true ), + m_CorrectEndPoints( true ), + m_SafetyNeighborhoodSize( 3 ) { this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( ); this->m_EndPointsIndex = this->m_LabelImageIndex + 1; @@ -203,39 +206,43 @@ template< class I, class O > void fpa::Image::DijkstraWithEndPointDetection< I, O >:: _BeforeGenerateData( ) { - const I* input = this->GetInput( ); - typename I::SpacingType spac = input->GetSpacing( ); - double max_spac = double( spac[ 0 ] ); - for( unsigned int d = 1; d < I::ImageDimension; ++d ) - max_spac = - ( double( spac[ d ] ) > max_spac )? double( spac[ d ] ): max_spac; - max_spac *= double( 3 ); - - // Correct seeds - for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s ) + if( this->m_CorrectSeeds ) { - TVertex seed = this->m_SeedVertices[ s ]; - _TNode n = this->m_Seeds[ seed ]; - _TRegion region = this->_Region( seed, max_spac ); - itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); - - iIt.GoToBegin( ); - _TPixel max_value = iIt.Get( ); - for( ++iIt; !iIt.IsAtEnd( ); ++iIt ) + const I* input = this->GetInput( ); + typename I::SpacingType spac = input->GetSpacing( ); + double max_spac = double( spac[ 0 ] ); + for( unsigned int d = 1; d < I::ImageDimension; ++d ) + max_spac = + ( double( spac[ d ] ) > max_spac )? double( spac[ d ] ): max_spac; + max_spac *= double( 3 ); + + // Correct seeds + for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s ) { - if( iIt.Get( ) > max_value ) + TVertex seed = this->m_SeedVertices[ s ]; + _TNode n = this->m_Seeds[ seed ]; + _TRegion region = this->_Region( seed, max_spac ); + itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); + + iIt.GoToBegin( ); + _TPixel max_value = iIt.Get( ); + for( ++iIt; !iIt.IsAtEnd( ); ++iIt ) { - this->m_SeedVertices[ s ] = iIt.GetIndex( ); - max_value = iIt.Get( ); + if( iIt.Get( ) > max_value ) + { + this->m_SeedVertices[ s ] = iIt.GetIndex( ); + max_value = iIt.Get( ); - } // fi + } // fi + + } // rof + this->m_Seeds.erase( seed ); + n.Parent = this->m_SeedVertices[ s ]; + this->m_Seeds[ this->m_SeedVertices[ s ] ] = n; } // rof - this->m_Seeds.erase( seed ); - n.Parent = this->m_SeedVertices[ s ]; - this->m_Seeds[ this->m_SeedVertices[ s ] ] = n; - } // rof + } // fi // End initialization this->Superclass::_BeforeGenerateData( ); @@ -312,7 +319,8 @@ _EndPointsAndBifurcations( ) // Correct it to nearest start candidate (high distance value) double vr = std::sqrt( double( input->GetPixel( v ) ) ); - v = this->_MaxInRegion( input, v, vr * double( 1.5 ) ); + if( this->m_CorrectEndPoints ) + v = this->_MaxInRegion( input, v, vr * double( 1.5 ) ); // Now, check for real marking conditions // 1. Has it been visited by dijkstra? @@ -507,7 +515,9 @@ _Region( const TVertex& c, const double& r ) for( unsigned int d = 0; d < I::ImageDimension; ++d ) { // NOTE: 3 is a minimum neighborhood size - long s = long( std::ceil( r / double( spac[ d ] ) ) ) + 3; + long s = + long( std::ceil( r / double( spac[ d ] ) ) ) + + long( this->m_SafetyNeighborhoodSize ); i0[ d ] = c[ d ] - s; i1[ d ] = c[ d ] + s; diff --git a/lib/fpa/Image/DijkstraWithSphereBacktracking.h b/lib/fpa/Image/DijkstraWithSphereBacktracking.h deleted file mode 100644 index 137017a..0000000 --- a/lib/fpa/Image/DijkstraWithSphereBacktracking.h +++ /dev/null @@ -1,97 +0,0 @@ -#ifndef __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__H__ -#define __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__H__ - -#include -#include -#include -#include - -namespace fpa -{ - namespace Image - { - /** - * @param I Input image type - */ - template< class I, class C > - class DijkstraWithSphereBacktracking - : public fpa::Image::Dijkstra< I, C > - { - public: - typedef DijkstraWithSphereBacktracking Self; - typedef fpa::Image::Dijkstra< I, C > Superclass; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - - typedef typename Superclass::TCost TCost; - typedef typename Superclass::TVertex TVertex; - typedef typename Superclass::InputImageType TImage; - typedef std::deque< TVertex > TVertices; - - typedef unsigned short TMark; - typedef itk::Image< TMark, I::ImageDimension > TMarkImage; - - typedef typename Superclass::TTraits::TVertexCmp TVertexCmp; - typedef std::pair< TVertex, TMark > TTreeNode; - typedef std::map< TVertex, TTreeNode, TVertexCmp > TTree; - - typedef typename Superclass::TEndEvent TEndEvent; - typedef typename Superclass::TBacktrackingEvent TBacktrackingEvent; - typedef typename Superclass::TEndBacktrackingEvent TEndBacktrackingEvent; - - protected: - typedef std::pair< TCost, TVertex > _TCandidate; - typedef std::multimap< TCost, TVertex > _TCandidates; - typedef typename Superclass::_TNode _TNode; - - typedef typename I::PixelType _TPixel; - typedef typename I::RegionType _TRegion; - typedef typename I::SizeType _TSize; - - public: - itkNewMacro( Self ); - itkTypeMacro( DijkstraWithSphereBacktracking, Dijkstra ); - - itkGetConstMacro( FullTree, TTree ); - itkGetConstMacro( ReducedTree, TTree ); - itkGetConstMacro( EndPoints, TVertices ); - itkGetConstMacro( BifurcationPoints, TVertices ); - itkGetConstMacro( NumberOfBranches, TMark ); - - public: - TMarkImage* GetOutputMarkImage( ); - const TMarkImage* GetOutputMarkImage( ) const; - - protected: - DijkstraWithSphereBacktracking( ); - virtual ~DijkstraWithSphereBacktracking( ); - - virtual void _BeforeMainLoop( ); - virtual void _AfterMainLoop( ); - virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n ); - virtual bool _UpdateResult( _TNode& n ); - - _TRegion _Region( const TVertex& c, const double& r ); - - private: - DijkstraWithSphereBacktracking( const Self& other ); - Self& operator=( const Self& other ); - - protected: - _TCandidates m_Candidates; - TTree m_FullTree; - TTree m_ReducedTree; - TVertices m_BifurcationPoints; - TVertices m_EndPoints; - TMark m_NumberOfBranches; - }; - - } // ecapseman - -} // ecapseman - -#include - -#endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__H__ - -// eof - $RCSfile$ diff --git a/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx b/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx deleted file mode 100644 index e8f6c51..0000000 --- a/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx +++ /dev/null @@ -1,402 +0,0 @@ -#ifndef __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__ -#define __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__ - -#include -#include -#include - -// ------------------------------------------------------------------------- -template< class I, class C > -typename fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -TMarkImage* fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -GetOutputMarkImage( ) -{ - return( - dynamic_cast< TMarkImage* >( - this->itk::ProcessObject::GetOutput( 1 ) - ) - ); -} - -// ------------------------------------------------------------------------- -template< class I, class C > -const typename fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -TMarkImage* fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -GetOutputMarkImage( ) const -{ - return( - dynamic_cast< const TMarkImage* >( - this->itk::ProcessObject::GetOutput( 1 ) - ) - ); -} - -// ------------------------------------------------------------------------- -template< class I, class C > -fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -DijkstraWithSphereBacktracking( ) - : Superclass( ), - m_NumberOfBranches( 0 ) -{ - this->SetNumberOfRequiredOutputs( 2 ); - this->SetNthOutput( 0, I::New( ) ); - this->SetNthOutput( 1, TMarkImage::New( ) ); -} - -// ------------------------------------------------------------------------- -template< class I, class C > -fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -~DijkstraWithSphereBacktracking( ) -{ -} - -// ------------------------------------------------------------------------- -template< class I, class C > -void fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -_BeforeMainLoop( ) -{ - const I* img = this->GetInput( ); - typename I::SpacingType spac = img->GetSpacing( ); - double max_spac = spac[ 0 ]; - for( unsigned int d = 1; d < I::ImageDimension; ++d ) - max_spac = - ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac; - max_spac *= double( 30 ); - - // Correct seeds - for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s ) - { - _TNode seed = this->m_Seeds[ s ]; - _TRegion region = this->_Region( seed.Vertex, max_spac ); - itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region ); - - iIt.GoToBegin( ); - _TPixel max_value = iIt.Get( ); - for( ++iIt; !iIt.IsAtEnd( ); ++iIt ) - { - if( iIt.Get( ) > max_value ) - { - seed.Vertex = iIt.GetIndex( ); - seed.Parent = seed.Vertex; - max_value = iIt.Get( ); - - } // fi - - } // rof - this->m_Seeds[ s ] = seed; - - } // rof - - // End initialization - this->Superclass::_BeforeMainLoop( ); - this->m_Candidates.clear( ); -} - -// ------------------------------------------------------------------------- -template< class I, class C > -void fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -_AfterMainLoop( ) -{ - // Finish base algorithm - this->Superclass::_AfterMainLoop( ); - this->m_FullTree.clear( ); - this->m_ReducedTree.clear( ); - this->m_EndPoints.clear( ); - this->m_BifurcationPoints.clear( ); - if( this->m_Candidates.size( ) == 0 ) - return; - this->InvokeEvent( TEndEvent( ) ); - - // Get some input values - const I* input = this->GetInput( ); - typename I::SpacingType spac = input->GetSpacing( ); - double max_spac = spac[ 0 ]; - for( unsigned int d = 1; d < I::ImageDimension; ++d ) - max_spac = - ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac; - max_spac *= double( 3 ); - - // Prepare an object to hold marks - typename TMarkImage::Pointer marks = this->GetOutputMarkImage( ); - marks->FillBuffer( 0 ); - - // Iterate over the candidates, starting from the candidates that - // are near thin branches - typename _TCandidates::const_reverse_iterator cIt = - this->m_Candidates.rbegin( ); - this->m_NumberOfBranches = 1; - for( ; cIt != this->m_Candidates.rend( ); ++cIt ) - { - // If pixel hasn't been visited, continue - TVertex v = cIt->second; - if( marks->GetPixel( v ) != 0 ) - continue; - - // Compute nearest start candidate - _TRegion region = this->_Region( v, max_spac ); - itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region ); - iIt.GoToBegin( ); - TVertex max_vertex = iIt.GetIndex( ); - _TPixel max_value = iIt.Get( ); - for( ++iIt; !iIt.IsAtEnd( ); ++iIt ) - { - _TPixel value = iIt.Get( ); - if( value > max_value ) - { - max_value = value; - max_vertex = iIt.GetIndex( ); - - } // fi - - } // rof - - // Keep endpoint - if( marks->GetPixel( max_vertex ) != 0 ) - continue; - - this->m_EndPoints.push_back( max_vertex ); - - // Construct path (at least the part that hasn't been iterated) - bool start = true; - bool change = false; - do - { - if( start ) - { - if( this->m_FullTree.find( max_vertex ) == this->m_FullTree.end( ) ) - { - // Mark a sphere around current point as visited - double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) ); - region = this->_Region( max_vertex, dist * double( 1.5 ) ); - itk::ImageRegionIteratorWithIndex< TMarkImage > - mIt( marks, region ); - for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) - mIt.Set( this->m_NumberOfBranches ); - - // Next vertex in current path - this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) ); - this->m_FullTree[ max_vertex ] = - TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); - // std::cout << "New: " << this->m_NumberOfBranches << std::endl; - } - else - { - // A bifurcation point has been reached! - this->m_BifurcationPoints.push_back( max_vertex ); - this->m_NumberOfBranches++; - - this->m_FullTree[ max_vertex ] = - TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); - // std::cout << "Bifurcation: " << this->m_NumberOfBranches << std::endl; - - start = false; - - } // fi - } - else - { - if( - std::find( - this->m_BifurcationPoints.begin( ), - this->m_BifurcationPoints.end( ), - max_vertex - ) == this->m_BifurcationPoints.end( ) - ) - { - // Mark a sphere around current point as visited - double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) ); - region = this->_Region( max_vertex, dist * double( 1.5 ) ); - itk::ImageRegionIteratorWithIndex< TMarkImage > - mIt( marks, region ); - for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) - mIt.Set( this->m_NumberOfBranches ); - - // Next vertex in current path - this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) ); - this->m_FullTree[ max_vertex ] = - TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); - change = true; - // std::cout << "Change: " << this->m_NumberOfBranches << std::endl; - } - else - { - // A bifurcation point has been reached! - // TODO: this->m_BifurcationPoints.push_back( max_vertex ); - this->m_NumberOfBranches++; - this->m_FullTree[ max_vertex ] = - TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); - // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl; - - } // fi - - } // fi - max_vertex = this->_Parent( max_vertex ); - - } while( max_vertex != this->_Parent( max_vertex ) ); - if( start || change ) - this->m_NumberOfBranches++; - this->m_FullTree[ max_vertex ] = TTreeNode( max_vertex, this->m_NumberOfBranches ); - - this->InvokeEvent( TEndBacktrackingEvent( ) ); - - } // rof - - std::map< TMark, unsigned long > histo; - for( - typename TTree::iterator treeIt = this->m_FullTree.begin( ); - treeIt != this->m_FullTree.end( ); - ++treeIt - ) - histo[ treeIt->second.second ]++; - - std::map< TMark, TMark > changes; - TMark last_change = 1; - for( TMark i = 1; i <= this->m_NumberOfBranches; ++i ) - { - if( histo[ i ] != 0 ) - changes[ i ] = last_change++; - - } // rof - this->m_NumberOfBranches = changes.size( ); - - for( - typename TTree::iterator treeIt = this->m_FullTree.begin( ); - treeIt != this->m_FullTree.end( ); - ++treeIt - ) - { - TMark old = treeIt->second.second; - treeIt->second.second = changes[ old ]; - - } // fi - - // Construct reduced tree - for( - typename TVertices::const_iterator eIt = this->m_EndPoints.begin( ); - eIt != this->m_EndPoints.end( ); - ++eIt - ) - { - typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt ); - if( tIt != this->m_FullTree.end( ) ) - { - TMark id = tIt->second.second; - do - { - tIt = this->m_FullTree.find( tIt->second.first ); - if( tIt == this->m_FullTree.end( ) ) - break; - - } while( tIt->second.second == id ); - this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id ); - - } // fi - - } // rof - - for( - typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( ); - bIt != this->m_BifurcationPoints.end( ); - ++bIt - ) - { - typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt ); - if( tIt != this->m_FullTree.end( ) ) - { - TMark id = tIt->second.second; - do - { - tIt = this->m_FullTree.find( tIt->second.first ); - if( tIt == this->m_FullTree.end( ) ) - break; - - } while( tIt->second.second == id ); - this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id ); - - } // fi - - } // rof -} - -// ------------------------------------------------------------------------- -template< class I, class C > -bool fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -_UpdateNeigh( _TNode& nn, const _TNode& n ) -{ - C nc = this->_Cost( nn.Vertex, n.Vertex ); - if( TCost( 0 ) < nc ) - { - typename I::PointType pnn, pn; - this->GetInput( )->TransformIndexToPhysicalPoint( nn.Vertex, pnn ); - this->GetInput( )->TransformIndexToPhysicalPoint( n.Vertex, pn ); - - - nc += TCost( 1 ); - nn.Cost = n.Cost + ( TCost( pnn.EuclideanDistanceTo( pn ) ) / std::pow( nc, 4 ) ); - nn.Result = nn.Cost; - return( true ); - } - else - return( false ); -} - -// ------------------------------------------------------------------------- -template< class I, class C > -bool fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -_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 ); -} - -// ------------------------------------------------------------------------- -template< class I, class C > -typename fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -_TRegion fpa::Image::DijkstraWithSphereBacktracking< I, C >:: -_Region( const TVertex& c, const double& r ) -{ - typename I::ConstPointer input = this->GetInput( ); - typename I::SpacingType spac = input->GetSpacing( ); - _TRegion region = input->GetLargestPossibleRegion( ); - typename I::IndexType idx0 = region.GetIndex( ); - typename I::IndexType idx1 = idx0 + region.GetSize( ); - - // Compute region size and index - typename I::IndexType i0, i1; - _TSize size; - for( unsigned int d = 0; d < I::ImageDimension; ++d ) - { - long s = long( std::ceil( r / double( spac[ d ] ) ) ); - i0[ d ] = c[ d ] - s; - i1[ d ] = c[ d ] + s; - - if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ]; - if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ]; - if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ]; - if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ]; - size[ d ] = i1[ d ] - i0[ d ]; - - } // rof - - // Prepare region and return it - region.SetIndex( i0 ); - region.SetSize( size ); - return( region ); -} - -#endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__ - -// eof - $RCSfile$ -- 2.47.1