From b63dc485b7255d1ab70ff72096beafe13a71f1be Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Mon, 9 Mar 2015 18:21:33 -0500 Subject: [PATCH] Skeletonization debugged --- .../example_ImageAlgorithmDijkstra_03.cxx | 414 +++--------------- ...example_ImageAlgorithm_Skeletonization.cxx | 112 ++++- lib/fpa/Base/Algorithm.h | 13 +- lib/fpa/Base/Events.h | 32 +- .../Image/DijkstraWithSphereBacktracking.h | 83 ++++ .../Image/DijkstraWithSphereBacktracking.hxx | 237 ++++++++++ lib/fpa/VTK/Image3DObserver.hxx | 25 +- lib/fpa/VTK/ImageMPR.cxx | 32 ++ lib/fpa/VTK/ImageMPR.h | 2 + 9 files changed, 561 insertions(+), 389 deletions(-) create mode 100644 lib/fpa/Image/DijkstraWithSphereBacktracking.h create mode 100644 lib/fpa/Image/DijkstraWithSphereBacktracking.hxx diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx index b3c04eb..b78929b 100644 --- a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx +++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx @@ -16,12 +16,11 @@ #include #include +#include #include #include -#include -#include -#include +#include #include #include @@ -34,332 +33,10 @@ typedef itk::ImageToVTKImageFilter< TImage > TVTKImage; typedef itk::ImageFileReader< TImage > TImageReader; typedef itk::ImageFileWriter< TImage > TImageWriter; -typedef fpa::Image::Dijkstra< TImage, TScalar > TBaseDijkstra; +typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra; typedef fpa::VTK::ImageMPR TMPR; -typedef fpa::VTK::Image3DObserver< TBaseDijkstra, vtkRenderWindow > TDijkstraObs; - -// ------------------------------------------------------------------------- -class TDijkstra - : public TBaseDijkstra -{ -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( ); - } - - virtual void _AfterMainLoop( ) - { - 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 - */ - } - - 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 fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs; // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) @@ -434,49 +111,60 @@ 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( ); + + const TDijkstra::TVertices& endpoints = paths->GetEndPoints( ); + const TDijkstra::TTree& tree = paths->GetFinalTree( ); + 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 ); - TImageWriter::Pointer dijkstra_writer = - TImageWriter::New( ); - dijkstra_writer->SetInput( paths->GetOutput( ) ); - dijkstra_writer->SetFileName( "dijkstra.mhd" ); - dijkstra_writer->Update( ); + TDijkstra::TVertex idx = *epIt; + do + { + TImage::PointType pnt; + input_image->TransformIndexToPhysicalPoint( idx, pnt ); - // Show result - /* - TVTKImage::Pointer output_image_vtk = TVTKImage::New( ); - output_image_vtk->SetInput( segmentor->GetOutput( ) ); - output_image_vtk->Update( ); + 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 ); - vtkSmartPointer< vtkImageMarchingCubes > mc = - vtkSmartPointer< vtkImageMarchingCubes >::New( ); - mc->SetInputData( output_image_vtk->GetOutput( ) ); - mc->SetValue( - 0, - double( segmentor->GetInsideValue( ) ) * double( 0.95 ) - ); - mc->Update( ); + } // fi + idx = tree.find( idx )->second; - // Let some interaction and close program - view.AddPolyData( mc->GetOutput( ), 0.1, 0.6, 0.8, 0.5 ); - view.Start( ); + } while( idx != tree.find( idx )->second ); - // 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 ); + } // rof - } // yrt + 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( ); + + /* TODO + TImageWriter::Pointer dijkstra_writer = + TImageWriter::New( ); + dijkstra_writer->SetInput( paths->GetOutput( ) ); + dijkstra_writer->SetFileName( "dijkstra.mhd" ); + dijkstra_writer->Update( ); */ + return( 0 ); } diff --git a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx index e0121bd..7bc80a9 100644 --- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx +++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx @@ -9,7 +9,10 @@ #include #include +#include + #include +#include #include #include @@ -20,6 +23,7 @@ typedef double 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,9 +32,13 @@ 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[] ) @@ -131,19 +139,101 @@ int main( int argc, char* argv[] ) 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( 1 ); + + 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->GetFinalTree( ); + 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; + + } while( idx != tree.find( idx )->second ); + + } // rof - TDistanceMapWriter::Pointer distancemap_writer = - TDistanceMapWriter::New( ); - distancemap_writer->SetInput( distanceMap->GetOutput( ) ); - distancemap_writer->SetFileName( output_distancemap_fn ); - distancemap_writer->Update( ); + vtkSmartPointer< vtkPolyData > vtk_tree = + vtkSmartPointer< vtkPolyData >::New( ); + vtk_tree->SetPoints( points ); + vtk_tree->SetLines( cells ); + vtk_tree->GetPointData( )->SetScalars( scalars ); - TImageWriter::Pointer segmentation_writer = - TImageWriter::New( ); - segmentation_writer->SetInput( segmentor->GetOutput( ) ); - segmentation_writer->SetFileName( output_segmentation_fn ); - segmentation_writer->Update( ); + view.AddPolyData( vtk_tree ); + view.Render( ); + view.Start( ); + + /* TODO + TDistanceMapWriter::Pointer distancemap_writer = + TDistanceMapWriter::New( ); + distancemap_writer->SetInput( distanceMap->GetOutput( ) ); + distancemap_writer->SetFileName( output_distancemap_fn ); + distancemap_writer->Update( ); + + TImageWriter::Pointer segmentation_writer = + TImageWriter::New( ); + segmentation_writer->SetInput( segmentor->GetOutput( ) ); + segmentation_writer->SetFileName( output_segmentation_fn ); + segmentation_writer->Update( ); + */ // Show result /* diff --git a/lib/fpa/Base/Algorithm.h b/lib/fpa/Base/Algorithm.h index 11cd447..c9b2e1a 100644 --- a/lib/fpa/Base/Algorithm.h +++ b/lib/fpa/Base/Algorithm.h @@ -50,12 +50,13 @@ namespace fpa typedef std::vector< _TCollisionSitesRow > _TCollisionSites; public: - typedef BaseEvent< _TNode > TEvent; - typedef FrontEvent< _TNode > TFrontEvent; - typedef MarkEvent< _TNode > TMarkEvent; - typedef CollisionEvent< _TNode > TCollisionEvent; - typedef EndEvent< _TNode > TEndEvent; - typedef BacktrackingEvent< TVertex > TBacktrackingEvent; + typedef BaseEvent< _TNode > TEvent; + typedef FrontEvent< _TNode > TFrontEvent; + typedef MarkEvent< _TNode > TMarkEvent; + typedef CollisionEvent< _TNode > TCollisionEvent; + typedef EndEvent< _TNode > TEndEvent; + typedef BacktrackingEvent< TVertex > TBacktrackingEvent; + typedef EndBacktrackingEvent< TVertex > TEndBacktrackingEvent; public: itkTypeMacro( Algorithm, itkProcessObject ); diff --git a/lib/fpa/Base/Events.h b/lib/fpa/Base/Events.h index ad662ab..c92c4db 100644 --- a/lib/fpa/Base/Events.h +++ b/lib/fpa/Base/Events.h @@ -139,7 +139,7 @@ namespace fpa { } BacktrackingEvent( const N& n, const unsigned long& id ) : BaseEvent< N >( n ), - BackId( id ) + BackId( id ) { } virtual ~BacktrackingEvent( ) { } @@ -157,6 +157,36 @@ namespace fpa unsigned long BackId; }; + /** + */ + template< class N > + class EndBacktrackingEvent + : public BaseEvent< N > + { + public: + EndBacktrackingEvent( ) + : BaseEvent< N >( ) + { } + EndBacktrackingEvent( const unsigned long& id ) + : BaseEvent< N >( ), + BackId( id ) + { } + virtual ~EndBacktrackingEvent( ) + { } + const char* GetEventName( ) const + { return( "fpa::Base::EndBacktrackingEvent" ); } + bool CheckEvent( const itk::EventObject* e ) const + { + return( + dynamic_cast< const EndBacktrackingEvent< N >* >( e ) != NULL + ); + } + itk::EventObject* MakeObject( ) const + { return( new EndBacktrackingEvent< N >( ) ); } + + unsigned long BackId; + }; + } // ecapseman } // ecapseman diff --git a/lib/fpa/Image/DijkstraWithSphereBacktracking.h b/lib/fpa/Image/DijkstraWithSphereBacktracking.h new file mode 100644 index 0000000..0cc0c4f --- /dev/null +++ b/lib/fpa/Image/DijkstraWithSphereBacktracking.h @@ -0,0 +1,83 @@ +#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 typename Superclass::TTraits::TVertexCmp TVertexCmp; + typedef std::map< TVertex, TVertex, 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( FinalTree, TTree ); + itkGetConstMacro( EndPoints, TVertices ); + + 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_FinalTree; + TVertices m_EndPoints; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx b/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx new file mode 100644 index 0000000..d71cef8 --- /dev/null +++ b/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx @@ -0,0 +1,237 @@ +#ifndef __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__ +#define __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__ + +#include +#include +#include + +// ------------------------------------------------------------------------- +template< class I, class C > +fpa::Image::DijkstraWithSphereBacktracking< I, C >:: +DijkstraWithSphereBacktracking( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +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( 3 ); + + // 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( ) +{ + typedef itk::Image< bool, I::ImageDimension > _TMarkImage; + + // Finish base algorithm + this->Superclass::_AfterMainLoop( ); + this->m_FinalTree.clear( ); + this->m_EndPoints.clear( ); + if( this->m_Candidates.size( ) == 0 ) + return; + this->InvokeEvent( TEndEvent( ) ); + + // Get 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 = _TMarkImage::New( ); + marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) ); + marks->SetRequestedRegion( input->GetRequestedRegion( ) ); + marks->SetBufferedRegion( input->GetBufferedRegion( ) ); + marks->SetOrigin( input->GetOrigin( ) ); + marks->SetSpacing( input->GetSpacing( ) ); + marks->SetDirection( input->GetDirection( ) ); + marks->Allocate( ); + marks->FillBuffer( false ); + + // Iterate over the candidates, starting from the candidates that + // are near thin branches + typename _TCandidates::const_reverse_iterator cIt = + this->m_Candidates.rbegin( ); + for( int backId = 0; cIt != this->m_Candidates.rend( ); ++cIt ) + { + // If pixel hasn't been visited, continue + TVertex v = cIt->second; + if( marks->GetPixel( v ) ) + 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 ) ) + continue; + backId++; + this->m_EndPoints.push_back( max_vertex ); + + // Construct path (at least the part that hasn't been iterated) + do + { + if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.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.25 ) ); + itk::ImageRegionIteratorWithIndex< _TMarkImage > + mIt( marks, region ); + for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) + mIt.Set( true ); + + // Next vertex in current path + this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) ); + this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex ); + + } // fi + max_vertex = this->_Parent( max_vertex ); + + } while( max_vertex != this->_Parent( max_vertex ) ); + this->m_FinalTree[ max_vertex ] = max_vertex; + this->InvokeEvent( TEndBacktrackingEvent( backId ) ); + + } // 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 ) + { + nn.Cost = n.Cost + ( TCost( 1 ) / nc ); + 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$ diff --git a/lib/fpa/VTK/Image3DObserver.hxx b/lib/fpa/VTK/Image3DObserver.hxx index 001924a..e9f9770 100644 --- a/lib/fpa/VTK/Image3DObserver.hxx +++ b/lib/fpa/VTK/Image3DObserver.hxx @@ -28,13 +28,14 @@ template< class F, class R > void fpa::VTK::Image3DObserver< F, R >:: Execute( const itk::Object* c, const itk::EventObject& e ) { - typedef itk::ImageBase< 3 > _TImage; - typedef typename F::TEvent _TEvent; - typedef typename F::TFrontEvent _TFrontEvent; - typedef typename F::TMarkEvent _TMarkEvent; - typedef typename F::TCollisionEvent _TCollisionEvent; - typedef typename F::TEndEvent _TEndEvent; - typedef typename F::TBacktrackingEvent _TBacktrackingEvent; + typedef itk::ImageBase< 3 > _TImage; + typedef typename F::TEvent _TEvent; + typedef typename F::TFrontEvent _TFrontEvent; + typedef typename F::TMarkEvent _TMarkEvent; + typedef typename F::TCollisionEvent _TCollisionEvent; + typedef typename F::TEndEvent _TEndEvent; + typedef typename F::TBacktrackingEvent _TBacktrackingEvent; + typedef typename F::TEndBacktrackingEvent _TEndBacktrackingEvent; // Check inputs if( this->m_RenderWindow == NULL ) @@ -158,6 +159,8 @@ Execute( const itk::Object* c, const itk::EventObject& e ) { const _TBacktrackingEvent* bevt = dynamic_cast< const _TBacktrackingEvent* >( &e ); + const _TEndBacktrackingEvent* ebevt = + dynamic_cast< const _TEndBacktrackingEvent* >( &e ); if( bevt != NULL ) { static const unsigned long nColors = 10; @@ -173,12 +176,18 @@ Execute( const itk::Object* c, const itk::EventObject& e ) this->m_Data->GetPointData( )-> GetScalars( )->InsertNextTuple1( back_id ); this->m_Data->Modified( ); - this->m_RenderWindow->Render( ); return; } // fi + if( ebevt != NULL ) + { + this->m_RenderWindow->Render( ); + return; + + } // fi + } // fi } diff --git a/lib/fpa/VTK/ImageMPR.cxx b/lib/fpa/VTK/ImageMPR.cxx index 3003395..8484a03 100644 --- a/lib/fpa/VTK/ImageMPR.cxx +++ b/lib/fpa/VTK/ImageMPR.cxx @@ -48,6 +48,15 @@ public: } // fi } break; + case 'z': + case 'Z': + { + this->SeedWidget->ProcessEventsOff( ); + this->WidgetX->InteractionOff( ); + this->WidgetY->InteractionOff( ); + this->WidgetZ->InteractionOff( ); + } + break; default: break; @@ -212,6 +221,22 @@ SetSize( unsigned int w, unsigned int h ) this->m_Window->SetSize( w, h ); } +// ------------------------------------------------------------------------- +void fpa::VTK::ImageMPR:: +AddPolyData( vtkPolyData* pd, double opacity ) +{ + unsigned int i = this->m_PolyDatas.size( ); + + this->m_PolyDatas.push_back( pd ); + this->m_Mappers.push_back( vtkSmartPointer< vtkPolyDataMapper >::New( ) ); + this->m_Actors.push_back( vtkSmartPointer< vtkActor >::New( ) ); + + this->m_Mappers[ i ]->SetInputData( pd ); + this->m_Actors[ i ]->SetMapper( this->m_Mappers[ i ] ); + this->m_Actors[ i ]->GetProperty( )->SetOpacity( opacity ); + this->m_Renderer->AddActor( this->m_Actors[ i ] ); +} + // ------------------------------------------------------------------------- void fpa::VTK::ImageMPR:: AddPolyData( vtkPolyData* pd, double r, double g, double b, double opacity ) @@ -278,4 +303,11 @@ Start( ) this->m_Interactor->Start( ); } +// ------------------------------------------------------------------------- +void fpa::VTK::ImageMPR:: +Render( ) +{ + this->m_Window->Render( ); +} + // eof - $RCSfile$ diff --git a/lib/fpa/VTK/ImageMPR.h b/lib/fpa/VTK/ImageMPR.h index 0ada0bb..21a1c5c 100644 --- a/lib/fpa/VTK/ImageMPR.h +++ b/lib/fpa/VTK/ImageMPR.h @@ -37,6 +37,7 @@ namespace fpa void SetBackground( double r, double g, double b ); void SetSize( unsigned int w, unsigned int h ); + void AddPolyData( vtkPolyData* pd, double opacity = double( 1 ) ); void AddPolyData( vtkPolyData* pd, double r, double g, double b, double opacity = double( 1 ) @@ -49,6 +50,7 @@ namespace fpa vtkRenderer* GetRenderer( ) const; void Start( ); + void Render( ); protected: vtkSmartPointer< vtkImageData > m_Image; -- 2.45.1