From 6fcc9fc78c44fa789bf092e2897cb6b391259b42 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Tue, 7 Apr 2015 18:59:01 -0500 Subject: [PATCH] Almost there... --- appli/examples/CMakeLists.txt | 1 + ...ample_Image_Dijkstra_EndPointDetection.cxx | 225 ++++++++++ .../example_Image_RegionGrow_AllPixels.cxx | 2 +- lib/fpa/Base/Algorithm.h | 60 ++- lib/fpa/Base/Algorithm.hxx | 221 ++++++--- lib/fpa/Base/Dijkstra.h | 52 ++- lib/fpa/Base/Dijkstra.hxx | 66 +-- lib/fpa/Base/MinimumSpanningTree.h | 30 +- lib/fpa/Base/MinimumSpanningTree.hxx | 29 +- lib/fpa/Base/RegionGrow.h | 31 +- lib/fpa/Base/RegionGrow.hxx | 39 +- lib/fpa/Image/Algorithm.h | 25 +- lib/fpa/Image/Algorithm.hxx | 71 +-- lib/fpa/Image/Dijkstra.h | 16 +- lib/fpa/Image/DijkstraWithEndPointDetection.h | 109 +++++ .../Image/DijkstraWithEndPointDetection.hxx | 420 ++++++++++++++++++ lib/fpa/Image/MinimumSpanningTree.h | 54 --- lib/fpa/Image/MinimumSpanningTree.hxx | 41 -- lib/fpa/Image/RegionGrow.h | 5 +- lib/fpa/VTK/Image2DObserver.hxx | 73 +-- lib/fpa/VTK/Image3DObserver.hxx | 225 +++++----- 21 files changed, 1244 insertions(+), 551 deletions(-) create mode 100644 appli/examples/example_Image_Dijkstra_EndPointDetection.cxx create mode 100644 lib/fpa/Image/DijkstraWithEndPointDetection.h create mode 100644 lib/fpa/Image/DijkstraWithEndPointDetection.hxx delete mode 100644 lib/fpa/Image/MinimumSpanningTree.h delete mode 100644 lib/fpa/Image/MinimumSpanningTree.hxx diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt index 5efc077..61768cb 100644 --- a/appli/examples/CMakeLists.txt +++ b/appli/examples/CMakeLists.txt @@ -8,6 +8,7 @@ IF(USE_VTK) example_Image_Dijkstra_CostFromRGBInput example_Image_Dijkstra_DanielssonCost example_Image_Dijkstra_DanielssonCost_TwoSeedsPath + example_Image_Dijkstra_EndPointDetection ) FOREACH(EX ${SIMPLE_VTK_EXAMPLES}) ADD_EXECUTABLE(${EX} ${EX}.cxx) diff --git a/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx b/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx new file mode 100644 index 0000000..070a3e8 --- /dev/null +++ b/appli/examples/example_Image_Dijkstra_EndPointDetection.cxx @@ -0,0 +1,225 @@ +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include + +/* + #include + #include + #include + #include + + #include + #include + #include + #include + #include + #include + #include + #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, class O > +void DistanceMap( + const typename I::Pointer& input, typename O::Pointer& output + ); + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 2 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image" + << std::endl; + return( 1 ); + + } // fi + std::string input_image_fn = argv[ 1 ]; + + // 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( 800, 800 ); + view.SetImage( vtk_input_image->GetOutput( ) ); + + // Allow some interaction and wait for, at least, one seed + view.Render( ); + while( view.GetNumberOfSeeds( ) == 0 ) + view.Start( ); + + // Get seed + double p[ 3 ]; + view.GetSeed( 0, p ); + TImage::PointType pnt; + TImage::IndexType seed; + pnt[ 0 ] = TImage::PointType::ValueType( p[ 0 ] ); + pnt[ 1 ] = TImage::PointType::ValueType( p[ 1 ] ); + pnt[ 2 ] = TImage::PointType::ValueType( p[ 2 ] ); + input_image->TransformPhysicalPointToIndex( pnt, seed ); + + // Compute squared distance map + TScalarImage::Pointer dmap; + DistanceMap< TImage, TScalarImage >( input_image, dmap ); + + // Prepare cost conversion function + typedef fpa::Base::Functors::InvertCostFunction< TScalar > TFunction; + TFunction::Pointer function = TFunction::New( ); + + // Prepare Dijkstra filter + typedef fpa::Image::DijkstraWithEndPointDetection< TScalarImage, TScalarImage > TFilter; + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( dmap ); + filter->SetConversionFunction( function ); + filter->SetNeighborhoodOrder( 2 ); + filter->StopAtOneFrontOff( ); + 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.0001 ); + filter->AddObserver( itk::AnyEvent( ), debugger ); + filter->ThrowEventsOn( ); + + // Go! + std::clock_t start = std::clock( ); + filter->Update( ); + std::clock_t end = std::clock( ); + std::cout + << "Extraction time = " + << ( double( end - start ) / double( CLOCKS_PER_SEC ) ) + << " s." << std::endl; + + /* TODO + // Save final total cost map + itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer = + itk::ImageFileWriter< TOutputImage >::New( ); + output_image_writer->SetFileName( output_image_fn ); + output_image_writer->SetInput( filter->GetOutput( ) ); + try + { + output_image_writer->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr + << "Error while writing image to " << output_image_fn << ": " + << err << std::endl; + return( 1 ); + + } // yrt + */ + 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( ); + + typename itk::MinimumMaximumImageCalculator< I >::Pointer minmax = + itk::MinimumMaximumImageCalculator< I >::New( ); + minmax->SetImage( reader->GetOutput( ) ); + minmax->Compute( ); + double vmin = double( minmax->GetMinimum( ) ); + double vmax = double( minmax->GetMaximum( ) ); + + typename itk::ShiftScaleImageFilter< I, I >::Pointer shift = + itk::ShiftScaleImageFilter< I, I >::New( ); + shift->SetInput( reader->GetOutput( ) ); + shift->SetScale( vmax - vmin ); + shift->SetShift( vmin / ( vmax - vmin ) ); + shift->Update( ); + + image = shift->GetOutput( ); + image->DisconnectPipeline( ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void DistanceMap( + const typename I::Pointer& input, typename O::Pointer& output + ) +{ + typename itk::SignedMaurerDistanceMapImageFilter< I, O >::Pointer dmap = + itk::SignedMaurerDistanceMapImageFilter< I, O >::New( ); + dmap->SetInput( input ); + dmap->SetBackgroundValue( ( typename I::PixelType )( 0 ) ); + dmap->InsideIsPositiveOn( ); + dmap->SquaredDistanceOn( ); + dmap->UseImageSpacingOn( ); + + std::clock_t start = std::clock( ); + dmap->Update( ); + std::clock_t end = std::clock( ); + std::cout + << "Distance map time = " + << ( double( end - start ) / double( CLOCKS_PER_SEC ) ) + << " s." << std::endl; + + output = dmap->GetOutput( ); + output->DisconnectPipeline( ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_Image_RegionGrow_AllPixels.cxx b/appli/examples/example_Image_RegionGrow_AllPixels.cxx index fffb58a..9484957 100644 --- a/appli/examples/example_Image_RegionGrow_AllPixels.cxx +++ b/appli/examples/example_Image_RegionGrow_AllPixels.cxx @@ -150,7 +150,7 @@ int main( int argc, char* argv[] ) typedef fpa::VTK::Image2DObserver< TFilter, vtkRenderWindow > TDebugger; TDebugger::Pointer debugger = TDebugger::New( ); debugger->SetRenderWindow( window ); - debugger->SetRenderPercentage( 0.001 ); + debugger->SetRenderPercentage( 0.01 ); filter->AddObserver( itk::AnyEvent( ), debugger ); filter->ThrowEventsOn( ); diff --git a/lib/fpa/Base/Algorithm.h b/lib/fpa/Base/Algorithm.h index 35c6a37..cc8497d 100644 --- a/lib/fpa/Base/Algorithm.h +++ b/lib/fpa/Base/Algorithm.h @@ -1,9 +1,11 @@ #ifndef __FPA__BASE__ALGORITHM__H__ #define __FPA__BASE__ALGORITHM__H__ +#include #include #include #include +#include namespace fpa { @@ -15,14 +17,15 @@ namespace fpa * vertex could be marked as "visited", "in the front", "not yet there" * or "freezed". * - * @param V Vertex type. - * @param C Vertex value type. - * @param R Result value type. - * @param B Base class for this algorithm. It should be any itk-based - * filter (itk::ProcessObject). + * @param V Vertex type. + * @param C Vertex value type. + * @param R Result value type. + * @param VC Vertex lexicographical compare. + * @param B Base class for this algorithm. It should be any itk-based + * filter (itk::ProcessObject). * */ - template< class V, class C, class R, class B > + template< class V, class C, class R, class VC, class B > class Algorithm : public B { @@ -32,17 +35,21 @@ namespace fpa typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; - typedef V TVertex; - typedef C TValue; - typedef R TResult; + typedef V TVertex; + typedef C TValue; + typedef R TResult; + typedef VC TVertexCompare; fpa_Base_NewEvent( TStartEvent ); fpa_Base_NewEvent( TStartLoopEvent ); + fpa_Base_NewEvent( TStartBacktrackingEvent ); fpa_Base_NewEvent( TEndEvent ); fpa_Base_NewEvent( TEndLoopEvent ); + fpa_Base_NewEvent( TEndBacktrackingEvent ); fpa_Base_NewEventWithVertex( TAliveEvent, TVertex ); fpa_Base_NewEventWithVertex( TFrontEvent, TVertex ); fpa_Base_NewEventWithVertex( TFreezeEvent, TVertex ); + fpa_Base_NewEventWithVertex( TBacktrackingEvent, TVertex ); protected: typedef std::vector< TVertex > _TVertices; @@ -68,13 +75,15 @@ namespace fpa virtual ~_TNode( ); public: - TVertex Vertex; - TVertex Parent; - TResult Result; - long FrontId; - _TNodeLabel Label; + TVertex Parent; + TResult Result; + short FrontId; + char Label; }; - typedef std::vector< _TNode > _TNodes; + typedef std::map< TVertex, _TNode, TVertexCompare > _TNodes; + + public: + typedef fpa::Base::MinimumSpanningTree< TVertex, _TCollisions, TVertexCompare > TMinimumSpanningTree; public: itkTypeMacro( Algorithm, B ); @@ -89,6 +98,10 @@ namespace fpa itkSetMacro( StopAtOneFront, bool ); public: + TMinimumSpanningTree* GetMinimumSpanningTree( ); + const TMinimumSpanningTree* GetMinimumSpanningTree( ) const; + void GraftMinimumSpanningTree( itk::DataObject* obj ); + virtual void InvokeEvent( const itk::EventObject& e ); virtual void InvokeEvent( const itk::EventObject& e ) const; @@ -134,18 +147,19 @@ namespace fpa ) const = 0; virtual void _InitResults( ) = 0; virtual const TResult& _Result( const TVertex& v ) const = 0; - virtual void _SetResult( const TVertex& v, const TResult& r ) = 0; + virtual void _SetResult( const TVertex& v, const _TNode& n ) = 0; // Marks-related abstract methods - virtual const _TNode& _Node( const TVertex& v ) const = 0; - virtual void _InitMarks( ) = 0; - virtual void _Mark( const _TNode& node ) = 0; + virtual _TNode& _Node( const TVertex& v ); + virtual const _TNode& _Node( const TVertex& v ) const; + virtual void _InitMarks( ); + virtual void _Mark( const TVertex& v, const _TNode& node ); // Queue-related abstract methods virtual void _InitQueue( ); virtual bool _IsQueueEmpty( ) const = 0; - virtual void _QueuePush( const _TNode& n ) = 0; - virtual _TNode _QueuePop( ) = 0; + virtual void _QueuePush( const TVertex& v, const _TNode& n ) = 0; + virtual void _QueuePop( TVertex& v, _TNode& n ) = 0; virtual void _QueueClear( ) = 0; private: @@ -157,7 +171,11 @@ namespace fpa bool m_ThrowEvents; bool m_StopAtOneFront; _TNodes m_Seeds; + _TVertices m_SeedVertices; + _TNodes m_Marks; _TCollisions m_Collisions; + + unsigned int m_MinimumSpanningTreeIndex; }; } // ecapseman diff --git a/lib/fpa/Base/Algorithm.hxx b/lib/fpa/Base/Algorithm.hxx index eb30758..dbcdc60 100644 --- a/lib/fpa/Base/Algorithm.hxx +++ b/lib/fpa/Base/Algorithm.hxx @@ -2,10 +2,11 @@ #define __FPA__BASE__ALGORITHM__HXX__ #include +#include // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -fpa::Base::Algorithm< V, C, R, B >::_TNode:: +template< class V, class C, class R, class VC, class B > +fpa::Base::Algorithm< V, C, R, VC, B >::_TNode:: _TNode( ) : Result( TResult( 0 ) ), FrontId( -1 ), @@ -14,15 +15,51 @@ _TNode( ) } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -fpa::Base::Algorithm< V, C, R, B >::_TNode:: +template< class V, class C, class R, class VC, class B > +fpa::Base::Algorithm< V, C, R, VC, B >::_TNode:: ~_TNode( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +typename fpa::Base::Algorithm< V, C, R, VC, B >:: +TMinimumSpanningTree* fpa::Base::Algorithm< V, C, R, VC, B >:: +GetMinimumSpanningTree( ) +{ + return( + dynamic_cast< TMinimumSpanningTree* >( + this->itk::ProcessObject::GetOutput( 1 ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class R, class VC, class B > +const typename fpa::Base::Algorithm< V, C, R, VC, B >:: +TMinimumSpanningTree* fpa::Base::Algorithm< V, C, R, VC, B >:: +GetMinimumSpanningTree( ) const +{ + return( + dynamic_cast< const TMinimumSpanningTree* >( + this->itk::ProcessObject::GetOutput( 1 ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: +GraftMinimumSpanningTree( itk::DataObject* obj ) +{ + TMinimumSpanningTree* mst = dynamic_cast< TMinimumSpanningTree* >( obj ); + if( mst != NULL ) + this->GraftNthOutput( 1, mst ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: InvokeEvent( const itk::EventObject& e ) { if( this->m_ThrowEvents ) @@ -30,8 +67,8 @@ InvokeEvent( const itk::EventObject& e ) } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: InvokeEvent( const itk::EventObject& e ) const { if( this->m_ThrowEvents ) @@ -39,66 +76,72 @@ InvokeEvent( const itk::EventObject& e ) const } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: AddSeed( const TVertex& s, const TResult& r ) { _TNode ns; - ns.Vertex = s; ns.Parent = s; ns.Result = r; - ns.FrontId = this->m_Seeds.size( ); + ns.FrontId = this->m_SeedVertices.size( ); ns.Label = Self::FrontLabel; - this->m_Seeds.push_back( ns ); + this->m_Seeds[ s ] = ns; + this->m_SeedVertices.push_back( s ); this->Modified( ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -const typename fpa::Base::Algorithm< V, C, R, B >:: -TVertex& fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +const typename fpa::Base::Algorithm< V, C, R, VC, B >:: +TVertex& fpa::Base::Algorithm< V, C, R, VC, B >:: GetSeed( const unsigned int& id ) const { - return( this->m_Seeds[ id ].Vertex ); + return( this->m_SeedVertices[ id ] ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: ClearSeeds( ) { this->m_Seeds.clear( ); + this->m_SeedVertices.clear( ); this->Modified( ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -unsigned long fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +unsigned long fpa::Base::Algorithm< V, C, R, VC, B >:: GetNumberOfSeeds( ) const { return( this->m_Seeds.size( ) ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +fpa::Base::Algorithm< V, C, R, VC, B >:: Algorithm( ) : Superclass( ), m_ThrowEvents( false ), m_StopAtOneFront( false ) { + this->m_MinimumSpanningTreeIndex = this->GetNumberOfRequiredOutputs( ); + this->SetNumberOfRequiredOutputs( this->m_MinimumSpanningTreeIndex + 1 ); + this->itk::ProcessObject::SetNthOutput( + this->m_MinimumSpanningTreeIndex, TMinimumSpanningTree::New( ) + ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +fpa::Base::Algorithm< V, C, R, VC, B >:: ~Algorithm( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: GenerateData( ) { unsigned long N = this->m_Seeds.size( ); @@ -120,8 +163,8 @@ GenerateData( ) } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: _Loop( ) { this->InvokeEvent( TStartLoopEvent( ) ); @@ -129,33 +172,34 @@ _Loop( ) while( !( this->_IsQueueEmpty( ) ) ) { // Get next candidate - _TNode candidate = this->_QueuePop( ); - if( this->_Node( candidate.Vertex ).Label == Self::AliveLabel ) + TVertex vertex; + _TNode node; + this->_QueuePop( vertex, node ); + if( this->_Node( vertex ).Label == Self::AliveLabel ) continue; // Mark it as "Alive" and update final result - candidate.Label = Self::AliveLabel; - this->_Mark( candidate ); - this->_SetResult( candidate.Vertex, candidate.Result ); - this->InvokeEvent( TAliveEvent( candidate.Vertex, candidate.FrontId ) ); + node.Label = Self::AliveLabel; + this->_Mark( vertex, node ); + this->_SetResult( vertex, node ); + this->InvokeEvent( TAliveEvent( vertex, node.FrontId ) ); // Check if a forced stop condition arises if( !( this->_NeedToStop( ) ) ) { // Compute neighborhood _TVertices neighborhood; - this->_Neighborhood( neighborhood, candidate.Vertex ); + this->_Neighborhood( neighborhood, vertex ); // Iterate over neighbors typename _TVertices::iterator nIt = neighborhood.begin( ); while( nIt != neighborhood.end( ) ) { _TNode neighbor = this->_Node( *nIt ); - neighbor.Vertex = *nIt; if( neighbor.Label == Self::AliveLabel ) { // Update collisions - if( this->_UpdateCollisions( candidate.Vertex, *nIt ) ) + if( this->_UpdateCollisions( vertex, *nIt ) ) { this->_QueueClear( ); nIt = neighborhood.end( ); @@ -166,21 +210,17 @@ _Loop( ) else { // Add new candidate to queue - if( - this->_ComputeNeighborResult( - neighbor.Result, *nIt, candidate.Vertex - ) - ) + if( this->_ComputeNeighborResult( neighbor.Result, *nIt, vertex ) ) { - neighbor.FrontId = candidate.FrontId; - neighbor.Parent = candidate.Vertex; + neighbor.FrontId = node.FrontId; + neighbor.Parent = vertex; neighbor.Label = Self::FrontLabel; - this->_QueuePush( neighbor ); - this->_Mark( neighbor ); - this->InvokeEvent( TFrontEvent( *nIt, candidate.FrontId ) ); + this->_QueuePush( *nIt, neighbor ); + this->_Mark( *nIt, neighbor ); + this->InvokeEvent( TFrontEvent( *nIt, node.FrontId ) ); } else - this->InvokeEvent( TFreezeEvent( *nIt, candidate.FrontId ) ); + this->InvokeEvent( TFreezeEvent( *nIt, node.FrontId ) ); } // fi ++nIt; @@ -196,36 +236,36 @@ _Loop( ) } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: _BeforeGenerateData( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: _AfterGenerateData( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: _BeforeLoop( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: _AfterLoop( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -bool fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +bool fpa::Base::Algorithm< V, C, R, VC, B >:: _UpdateCollisions( const TVertex& a, const TVertex& b ) { bool ret = false; @@ -242,9 +282,9 @@ _UpdateCollisions( const TVertex& a, const TVertex& b ) if( !exists ) { - this->m_Collisions[ fa ][ fb ].first = na.Vertex; + this->m_Collisions[ fa ][ fb ].first = a; this->m_Collisions[ fa ][ fb ].second = true; - this->m_Collisions[ fb ][ fa ].first = nb.Vertex; + this->m_Collisions[ fb ][ fa ].first = b; this->m_Collisions[ fb ][ fa ].second = true; // Stop if one front is desired @@ -282,16 +322,63 @@ _UpdateCollisions( const TVertex& a, const TVertex& b ) } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -bool fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +bool fpa::Base::Algorithm< V, C, R, VC, B >:: _NeedToStop( ) const { return( false ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Algorithm< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +typename fpa::Base::Algorithm< V, C, R, VC, B >:: +_TNode& fpa::Base::Algorithm< V, C, R, VC, B >:: +_Node( const TVertex& v ) +{ + typename _TNodes::iterator vIt = this->m_Marks.find( v ); + if( vIt == this->m_Marks.end( ) ) + { + _TNode new_node; + new_node.Parent = v; + new_node.Result = TResult( 0 ); + new_node.FrontId = -1; + new_node.Label = Self::FarLabel; + this->m_Marks[ v ] = new_node; + vIt = this->m_Marks.find( v ); + + } // fi + return( vIt->second ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class R, class VC, class B > +const typename fpa::Base::Algorithm< V, C, R, VC, B >:: +_TNode& fpa::Base::Algorithm< V, C, R, VC, B >:: +_Node( const TVertex& v ) const +{ + typename _TNodes::const_iterator vIt = this->m_Marks.find( v ); + return( vIt->second ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: +_InitMarks( ) +{ + this->m_Marks.clear( ); +} + +// ------------------------------------------------------------------------- +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: +_Mark( const TVertex& v, const _TNode& node ) +{ + this->m_Marks[ v ] = node; +} + +// ------------------------------------------------------------------------- +template< class V, class C, class R, class VC, class B > +void fpa::Base::Algorithm< V, C, R, VC, B >:: _InitQueue( ) { this->_QueueClear( ); @@ -301,8 +388,8 @@ _InitQueue( ) sIt++ ) { - this->_QueuePush( *sIt ); - this->_Mark( *sIt ); + this->_QueuePush( sIt->first, sIt->second ); + this->_Mark( sIt->first, sIt->second ); } // rof } diff --git a/lib/fpa/Base/Dijkstra.h b/lib/fpa/Base/Dijkstra.h index 3831526..cfa06b3 100644 --- a/lib/fpa/Base/Dijkstra.h +++ b/lib/fpa/Base/Dijkstra.h @@ -11,26 +11,40 @@ namespace fpa /** * Dijkstra is a front propagation algorithm that minimizes costs * - * @param V Vertex type. - * @param C Vertex value type. - * @param R Result value type. - * @param B Base class for this algorithm. It should be any itk-based - * filter (itk::ProcessObject). + * @param V Vertex type. + * @param C Vertex value type. + * @param R Result value type. + * @param VC Vertex lexicographical compare. + * @param B Base class for this algorithm. It should be any itk-based + * filter (itk::ProcessObject). * */ - template< class V, class C, class R, class B > + template< class V, class C, class R, class VC, class B > class Dijkstra - : public Algorithm< V, C, R, B > + : public Algorithm< V, C, R, VC, B > { public: typedef Dijkstra Self; - typedef Algorithm< V, C, R, B > Superclass; + typedef Algorithm< V, C, R, VC, B > Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; - typedef typename Superclass::TVertex TVertex; - typedef typename Superclass::TValue TValue; - typedef typename Superclass::TResult TResult; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TValue TValue; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TVertexCompare TVertexCompare; + + typedef typename Superclass::TStartEvent TStartEvent; + typedef typename Superclass::TStartLoopEvent TStartLoopEvent; + typedef typename Superclass::TEndEvent TEndEvent; + typedef typename Superclass::TEndLoopEvent TEndLoopEvent; + typedef typename Superclass::TAliveEvent TAliveEvent; + typedef typename Superclass::TFrontEvent TFrontEvent; + typedef typename Superclass::TFreezeEvent TFreezeEvent; + + typedef typename Superclass::TStartBacktrackingEvent TStartBacktrackingEvent; + typedef typename Superclass::TEndBacktrackingEvent TEndBacktrackingEvent; + typedef typename Superclass::TBacktrackingEvent TBacktrackingEvent; protected: typedef typename Superclass::_TVertices _TVertices; @@ -40,13 +54,16 @@ namespace fpa typedef typename Superclass::_TNode _TNode; typedef typename Superclass::_TNodes _TNodes; - typedef std::vector< _TNode > _TQueue; - struct _TNodeCompare + struct _TQueueNode { + TVertex Vertex; + _TNode Node; + // Make the min-heap behave as a max-heap - bool operator()( const _TNode& a, const _TNode& b ) const - { return( b.Result < a.Result ); } + bool operator<( const _TQueueNode& other ) const + { return( other.Node.Result < this->Node.Result ); } }; + typedef std::vector< _TQueueNode > _TQueue; public: itkTypeMacro( Dijkstra, Algorithm ); @@ -64,8 +81,8 @@ namespace fpa // Queue-related abstract methods virtual bool _IsQueueEmpty( ) const; - virtual void _QueuePush( const _TNode& n ); - virtual _TNode _QueuePop( ); + virtual void _QueuePush( const TVertex& v, const _TNode& n ); + virtual void _QueuePop( TVertex& v, _TNode& n ); virtual void _QueueClear( ); private: @@ -75,7 +92,6 @@ namespace fpa protected: _TQueue m_Queue; - static _TNodeCompare m_NodeCompare; }; } // ecapseman diff --git a/lib/fpa/Base/Dijkstra.hxx b/lib/fpa/Base/Dijkstra.hxx index 7d65237..2a2ccf5 100644 --- a/lib/fpa/Base/Dijkstra.hxx +++ b/lib/fpa/Base/Dijkstra.hxx @@ -4,23 +4,23 @@ #include // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -fpa::Base::Dijkstra< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +fpa::Base::Dijkstra< V, C, R, VC, B >:: Dijkstra( ) : Superclass( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -fpa::Base::Dijkstra< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +fpa::Base::Dijkstra< V, C, R, VC, B >:: ~Dijkstra( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -bool fpa::Base::Dijkstra< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +bool fpa::Base::Dijkstra< V, C, R, VC, B >:: _ComputeNeighborResult( TResult& result, const TVertex& neighbor, const TVertex& parent ) const @@ -28,49 +28,53 @@ _ComputeNeighborResult( result = this->_Cost( neighbor, parent ); result *= TResult( this->_Distance( neighbor, parent ) ); - _TNode pn = this->_Node( parent ); - if( pn.Label == Self::AliveLabel ) - result += pn.Result; - - return( true ); + // WARNING: Dijkstra does not work on negative values! + if( result >= TResult( 0 ) ) + { + _TNode pn = this->_Node( parent ); + if( pn.Label == Self::AliveLabel ) + result += pn.Result; + return( true ); + } + else + return( false ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -bool fpa::Base::Dijkstra< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +bool fpa::Base::Dijkstra< V, C, R, VC, B >:: _IsQueueEmpty( ) const { return( this->m_Queue.empty( ) ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Dijkstra< V, C, R, B >:: -_QueuePush( const _TNode& n ) +template< class V, class C, class R, class VC, class B > +void fpa::Base::Dijkstra< V, C, R, VC, B >:: +_QueuePush( const TVertex& v, const _TNode& n ) { - this->m_Queue.push_back( n ); - std::push_heap( - this->m_Queue.begin( ), this->m_Queue.end( ), Self::m_NodeCompare - ); + _TQueueNode qn; + qn.Vertex = v; + qn.Node = n; + this->m_Queue.push_back( qn ); + std::push_heap( this->m_Queue.begin( ), this->m_Queue.end( ) ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -typename fpa::Base::Dijkstra< V, C, R, B >:: -_TNode fpa::Base::Dijkstra< V, C, R, B >:: -_QueuePop( ) +template< class V, class C, class R, class VC, class B > +void fpa::Base::Dijkstra< V, C, R, VC, B >:: +_QueuePop( TVertex& v, _TNode& n ) { - _TNode n = this->m_Queue.front( ); - std::pop_heap( - this->m_Queue.begin( ), this->m_Queue.end( ), Self::m_NodeCompare - ); + _TQueueNode qn = this->m_Queue.front( ); + std::pop_heap( this->m_Queue.begin( ), this->m_Queue.end( ) ); this->m_Queue.pop_back( ); - return( n ); + v = qn.Vertex; + n = qn.Node; } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::Dijkstra< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::Dijkstra< V, C, R, VC, B >:: _QueueClear( ) { this->m_Queue.clear( ); diff --git a/lib/fpa/Base/MinimumSpanningTree.h b/lib/fpa/Base/MinimumSpanningTree.h index ec744e6..9adcd6d 100644 --- a/lib/fpa/Base/MinimumSpanningTree.h +++ b/lib/fpa/Base/MinimumSpanningTree.h @@ -1,7 +1,10 @@ #ifndef __FPA__BASE__MINIMUMSPANNINGTREE__H__ #define __FPA__BASE__MINIMUMSPANNINGTREE__H__ +#include +#include #include +#include #include namespace fpa @@ -10,30 +13,39 @@ namespace fpa { /** */ - template< class V, class C, class B > + template< class V, class C, class VC > class MinimumSpanningTree - : public B + : public itk::SimpleDataObjectDecorator< std::map< V, std::pair< V, short >, VC > > { public: - typedef MinimumSpanningTree Self; - typedef B Superclass; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; + typedef std::pair< V, short > TNodeInfo; + typedef std::map< V, TNodeInfo, VC > TDecorated; + typedef MinimumSpanningTree Self; + typedef itk::SimpleDataObjectDecorator< TDecorated > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; - typedef V TVertex; - typedef C TCollisions; + typedef V TVertex; + typedef C TCollisions; + typedef VC TVertexCompare; protected: typedef std::vector< unsigned long > _TRow; typedef std::vector< _TRow > _TMatrix; public: + itkNewMacro( Self ); itkTypeMacro( MinimumSpanningTree, B ); itkGetConstMacro( Collisions, TCollisions ); public: void SetCollisions( const TCollisions& collisions ); + void SetParent( const TVertex& v, const TVertex& p, const short& fid ) + { + this->Get( )[ v ] = TNodeInfo( p, fid ); + this->Modified( ); + } virtual void GetPath( std::vector< V >& path, const V& a, const V& b @@ -44,8 +56,6 @@ namespace fpa virtual ~MinimumSpanningTree( ); virtual void _Path( std::vector< V >& path, const V& a ) const; - virtual long _FrontId( const V& v ) const = 0; - virtual V _Parent( const V& v ) const = 0; private: // Purposely not implemented diff --git a/lib/fpa/Base/MinimumSpanningTree.hxx b/lib/fpa/Base/MinimumSpanningTree.hxx index 1140e40..d932247 100644 --- a/lib/fpa/Base/MinimumSpanningTree.hxx +++ b/lib/fpa/Base/MinimumSpanningTree.hxx @@ -56,7 +56,6 @@ SetCollisions( const TCollisions& collisions ) } // rof } // rof - this->Modified( ); } @@ -65,8 +64,14 @@ template< class V, class C, class B > void fpa::Base::MinimumSpanningTree< V, C, B >:: GetPath( std::vector< V >& path, const V& a, const V& b ) const { - long fa = this->_FrontId( a ); - long fb = this->_FrontId( b ); + typename TDecorated::const_iterator aIt = this->Get( ).find( a ); + typename TDecorated::const_iterator bIt = this->Get( ).find( b ); + + if( aIt == this->Get( ).end( ) || bIt == this->Get( ).end( ) ) + return; + + short fa = aIt->second.second; + short fb = bIt->second.second; if( fa == fb ) { @@ -164,14 +169,20 @@ template< class V, class C, class B > void fpa::Base::MinimumSpanningTree< V, C, B >:: _Path( std::vector< V >& path, const V& a ) const { - V it = a; - do + typename TDecorated::const_iterator dIt = this->Get( ).find( a ); + if( dIt != this->Get( ).end( ) ) { - path.push_back( it ); - it = this->_Parent( it ); + do + { + path.push_back( dIt->first ); + dIt = this->Get( ).find( dIt->second.first ); - } while( it != this->_Parent( it ) ); - path.push_back( it ); + } while( dIt->first != dIt->second.first && dIt != this->Get( ).end( ) ); + + if( dIt != this->Get( ).end( ) ) + path.push_back( dIt->first ); + + } // fi } #endif // __FPA__BASE__MINIMUMSPANNINGTREE__HXX__ diff --git a/lib/fpa/Base/RegionGrow.h b/lib/fpa/Base/RegionGrow.h index 6028b75..0b5b797 100644 --- a/lib/fpa/Base/RegionGrow.h +++ b/lib/fpa/Base/RegionGrow.h @@ -2,6 +2,7 @@ #define __FPA__BASE__REGIONGROW__H__ #include +#include #include namespace fpa @@ -11,26 +12,28 @@ namespace fpa /** * Region grow is a front propagation with no costs. * - * @param V Vertex type. - * @param C Vertex value type. - * @param R Result value type. - * @param B Base class for this algorithm. It should be any itk-based - * filter (itk::ProcessObject). + * @param V Vertex type. + * @param C Vertex value type. + * @param R Result value type. + * @param VC Vertex lexicographical compare. + * @param B Base class for this algorithm. It should be any itk-based + * filter (itk::ProcessObject). * */ - template< class V, class C, class R, class B > + template< class V, class C, class R, class VC, class B > class RegionGrow - : public Algorithm< V, C, R, B > + : public Algorithm< V, C, R, VC, B > { public: typedef RegionGrow Self; - typedef Algorithm< V, C, R, B > Superclass; + typedef Algorithm< V, C, R, VC, B > Superclass; typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; - typedef typename Superclass::TVertex TVertex; - typedef typename Superclass::TValue TValue; - typedef typename Superclass::TResult TResult; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TValue TValue; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TVertexCompare TVertexCompare; protected: typedef typename Superclass::_TVertices _TVertices; @@ -40,7 +43,7 @@ namespace fpa typedef typename Superclass::_TNode _TNode; typedef typename Superclass::_TNodes _TNodes; - typedef std::queue< _TNode > _TQueue; + typedef std::queue< std::pair< TVertex, _TNode > > _TQueue; public: itkTypeMacro( RegionGrow, Algorithm ); @@ -64,8 +67,8 @@ namespace fpa // Queue-related abstract methods virtual bool _IsQueueEmpty( ) const; - virtual void _QueuePush( const _TNode& n ); - virtual _TNode _QueuePop( ); + virtual void _QueuePush( const TVertex& v, const _TNode& n ); + virtual void _QueuePop( TVertex& v, _TNode& n ); virtual void _QueueClear( ); private: diff --git a/lib/fpa/Base/RegionGrow.hxx b/lib/fpa/Base/RegionGrow.hxx index 2d554c2..f0c4512 100644 --- a/lib/fpa/Base/RegionGrow.hxx +++ b/lib/fpa/Base/RegionGrow.hxx @@ -2,8 +2,8 @@ #define __FPA__BASE__REGIONGROWING__HXX__ // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -fpa::Base::RegionGrow< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +fpa::Base::RegionGrow< V, C, R, VC, B >:: RegionGrow( ) : Superclass( ), m_InsideValue( TResult( 1 ) ), @@ -12,15 +12,15 @@ RegionGrow( ) } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -fpa::Base::RegionGrow< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +fpa::Base::RegionGrow< V, C, R, VC, B >:: ~RegionGrow( ) { } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -bool fpa::Base::RegionGrow< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +bool fpa::Base::RegionGrow< V, C, R, VC, B >:: _ComputeNeighborResult( TResult& result, const TVertex& neighbor, const TVertex& parent ) const @@ -39,35 +39,34 @@ _ComputeNeighborResult( } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -bool fpa::Base::RegionGrow< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +bool fpa::Base::RegionGrow< V, C, R, VC, B >:: _IsQueueEmpty( ) const { return( this->m_Queue.empty( ) ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::RegionGrow< V, C, R, B >:: -_QueuePush( const _TNode& n ) +template< class V, class C, class R, class VC, class B > +void fpa::Base::RegionGrow< V, C, R, VC, B >:: +_QueuePush( const TVertex& v, const _TNode& n ) { - this->m_Queue.push( n ); + this->m_Queue.push( std::pair< TVertex, _TNode >( v, n ) ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -typename fpa::Base::RegionGrow< V, C, R, B >:: -_TNode fpa::Base::RegionGrow< V, C, R, B >:: -_QueuePop( ) +template< class V, class C, class R, class VC, class B > +void fpa::Base::RegionGrow< V, C, R, VC, B >:: +_QueuePop( TVertex& v, _TNode& n ) { - _TNode n = this->m_Queue.front( ); + v = this->m_Queue.front( ).first; + n = this->m_Queue.front( ).second; this->m_Queue.pop( ); - return( n ); } // ------------------------------------------------------------------------- -template< class V, class C, class R, class B > -void fpa::Base::RegionGrow< V, C, R, B >:: +template< class V, class C, class R, class VC, class B > +void fpa::Base::RegionGrow< V, C, R, VC, B >:: _QueueClear( ) { while( this->m_Queue.size( ) > 0 ) diff --git a/lib/fpa/Image/Algorithm.h b/lib/fpa/Image/Algorithm.h index e59bed2..f9ec7c4 100644 --- a/lib/fpa/Image/Algorithm.h +++ b/lib/fpa/Image/Algorithm.h @@ -2,7 +2,6 @@ #define __FPA__IMAGE__ALGORITHM__H__ #include -#include namespace fpa { @@ -30,9 +29,10 @@ namespace fpa typedef I TInputImage; typedef O TOutputImage; - typedef typename Superclass::TVertex TVertex; - typedef typename Superclass::TValue TValue; - typedef typename Superclass::TResult TResult; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TValue TValue; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TVertexCompare TVertexCompare; protected: typedef typename Superclass::_TVertices _TVertices; @@ -42,11 +42,6 @@ namespace fpa typedef typename Superclass::_TNode _TNode; typedef typename Superclass::_TNodes _TNodes; - typedef fpa::Image::MinimumSpanningTree< TVertex, _TNode, _TCollisions, I::ImageDimension, Self::AliveLabel > _TMarks; - - public: - typedef _TMarks TMinimumSpanningTree; - public: itkTypeMacro( Algorithm, TAlgorithm ); @@ -54,11 +49,6 @@ namespace fpa itkGetConstMacro( NeighborhoodOrder, unsigned int ); itkSetMacro( NeighborhoodOrder, unsigned int ); - public: - TMinimumSpanningTree* GetMinimumSpanningTree( ); - const TMinimumSpanningTree* GetMinimumSpanningTree( ) const; - void GraftMinimumSpanningTree( itk::DataObject* obj ); - protected: Algorithm( ); virtual ~Algorithm( ); @@ -80,12 +70,7 @@ namespace fpa // Results-related abstract methods virtual void _InitResults( ); virtual const TResult& _Result( const TVertex& v ) const; - virtual void _SetResult( const TVertex& v, const TResult& r ); - - // Marks-related abstract methods - virtual const _TNode& _Node( const TVertex& v ) const; - virtual void _InitMarks( ); - virtual void _Mark( const _TNode& node ); + virtual void _SetResult( const TVertex& v, const _TNode& n ); private: // Purposely not implemented diff --git a/lib/fpa/Image/Algorithm.hxx b/lib/fpa/Image/Algorithm.hxx index a5b1f85..5e8605d 100644 --- a/lib/fpa/Image/Algorithm.hxx +++ b/lib/fpa/Image/Algorithm.hxx @@ -4,42 +4,6 @@ #include #include -// ------------------------------------------------------------------------- -template< class I, class O, class A > -typename fpa::Image::Algorithm< I, O, A >:: -TMinimumSpanningTree* fpa::Image::Algorithm< I, O, A >:: -GetMinimumSpanningTree( ) -{ - return( - dynamic_cast< TMinimumSpanningTree* >( - this->itk::ProcessObject::GetOutput( 1 ) - ) - ); -} - -// ------------------------------------------------------------------------- -template< class I, class O, class A > -const typename fpa::Image::Algorithm< I, O, A >:: -TMinimumSpanningTree* fpa::Image::Algorithm< I, O, A >:: -GetMinimumSpanningTree( ) const -{ - return( - dynamic_cast< const TMinimumSpanningTree* >( - this->itk::ProcessObject::GetOutput( 1 ) - ) - ); -} - -// ------------------------------------------------------------------------- -template< class I, class O, class A > -void fpa::Image::Algorithm< I, O, A >:: -GraftMinimumSpanningTree( itk::DataObject* obj ) -{ - TMinimumSpanningTree* mst = dynamic_cast< TMinimumSpanningTree* >( obj ); - if( mst != NULL ) - this->GraftNthOutput( 1, mst ); -} - // ------------------------------------------------------------------------- template< class I, class O, class A > fpa::Image::Algorithm< I, O, A >:: @@ -47,9 +11,6 @@ Algorithm( ) : Superclass( ), m_NeighborhoodOrder( 1 ) { - this->itk::ProcessObject::SetNumberOfRequiredOutputs( 2 ); - this->itk::ProcessObject::SetNthOutput( 0, O::New( ) ); - this->itk::ProcessObject::SetNthOutput( 1, TMinimumSpanningTree::New( ) ); } // ------------------------------------------------------------------------- @@ -182,36 +143,10 @@ _Result( const TVertex& v ) const // ------------------------------------------------------------------------- template< class I, class O, class A > void fpa::Image::Algorithm< I, O, A >:: -_SetResult( const TVertex& v, const TResult& r ) -{ - this->GetOutput( )->SetPixel( v, r ); -} - -// ------------------------------------------------------------------------- -template< class I, class O, class A > -const typename fpa::Image::Algorithm< I, O, A >:: -_TNode& fpa::Image::Algorithm< I, O, A >:: -_Node( const TVertex& v ) const -{ - return( this->GetMinimumSpanningTree( )->GetPixel( v ) ); -} - -// ------------------------------------------------------------------------- -template< class I, class O, class A > -void fpa::Image::Algorithm< I, O, A >:: -_InitMarks( ) -{ - _TNode far_node; - far_node.Label = Self::FarLabel; - this->GetMinimumSpanningTree( )->FillBuffer( far_node ); -} - -// ------------------------------------------------------------------------- -template< class I, class O, class A > -void fpa::Image::Algorithm< I, O, A >:: -_Mark( const _TNode& node ) +_SetResult( const TVertex& v, const _TNode& n ) { - this->GetMinimumSpanningTree( )->SetPixel( node.Vertex, node ); + this->GetOutput( )->SetPixel( v, n.Result ); + this->GetMinimumSpanningTree( )->SetParent( v, n.Parent, n.FrontId ); } #endif // __FPA__IMAGE__ALGORITHM__HXX__ diff --git a/lib/fpa/Image/Dijkstra.h b/lib/fpa/Image/Dijkstra.h index 4e5105c..48d7a79 100644 --- a/lib/fpa/Image/Dijkstra.h +++ b/lib/fpa/Image/Dijkstra.h @@ -17,10 +17,10 @@ namespace fpa */ template< class I, class O > class Dijkstra - : public Algorithm< I, O, fpa::Base::Dijkstra< typename I::IndexType, typename I::PixelType, typename O::PixelType, itk::ImageToImageFilter< I, O > > > + : public Algorithm< I, O, fpa::Base::Dijkstra< typename I::IndexType, typename I::PixelType, typename O::PixelType, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, itk::ImageToImageFilter< I, O > > > { public: - typedef fpa::Base::Dijkstra< typename I::IndexType, typename I::PixelType, typename O::PixelType, itk::ImageToImageFilter< I, O > > TBaseAlgorithm; + typedef fpa::Base::Dijkstra< typename I::IndexType, typename I::PixelType, typename O::PixelType, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, itk::ImageToImageFilter< I, O > > TBaseAlgorithm; typedef Dijkstra Self; typedef Algorithm< I, O, TBaseAlgorithm > Superclass; @@ -33,6 +33,18 @@ namespace fpa typedef typename Superclass::TValue TValue; typedef typename Superclass::TResult TResult; + typedef typename Superclass::TStartEvent TStartEvent; + typedef typename Superclass::TStartLoopEvent TStartLoopEvent; + typedef typename Superclass::TEndEvent TEndEvent; + typedef typename Superclass::TEndLoopEvent TEndLoopEvent; + typedef typename Superclass::TAliveEvent TAliveEvent; + typedef typename Superclass::TFrontEvent TFrontEvent; + typedef typename Superclass::TFreezeEvent TFreezeEvent; + + typedef typename Superclass::TStartBacktrackingEvent TStartBacktrackingEvent; + typedef typename Superclass::TEndBacktrackingEvent TEndBacktrackingEvent; + typedef typename Superclass::TBacktrackingEvent TBacktrackingEvent; + typedef fpa::Image::Functors::ImageCostFunction< TInputImage, TResult > TCostFunction; typedef itk::FunctionBase< TResult, TResult > TConversionFunction; diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.h b/lib/fpa/Image/DijkstraWithEndPointDetection.h new file mode 100644 index 0000000..ecf197d --- /dev/null +++ b/lib/fpa/Image/DijkstraWithEndPointDetection.h @@ -0,0 +1,109 @@ +#ifndef __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__H__ +#define __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__H__ + +#include +#include + +namespace fpa +{ + namespace Image + { + /** + * @param I Input image type + * @param O Output image type + */ + template< class I, class O > + class DijkstraWithEndPointDetection + : public Dijkstra< I, O > + { + public: + typedef DijkstraWithEndPointDetection Self; + typedef Dijkstra< I, O > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TInputImage TInputImage; + typedef typename Superclass::TOutputImage TOutputImage; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TValue TValue; + typedef typename Superclass::TResult TResult; + typedef typename Superclass::TCostFunction TCostFunction; + typedef typename Superclass::TConversionFunction TConversionFunction; + typedef typename Superclass::_TVertices TVertices; + + typedef typename Superclass::TStartEvent TStartEvent; + typedef typename Superclass::TStartLoopEvent TStartLoopEvent; + typedef typename Superclass::TEndEvent TEndEvent; + typedef typename Superclass::TEndLoopEvent TEndLoopEvent; + typedef typename Superclass::TAliveEvent TAliveEvent; + typedef typename Superclass::TFrontEvent TFrontEvent; + typedef typename Superclass::TFreezeEvent TFreezeEvent; + + typedef typename Superclass::TStartBacktrackingEvent TStartBacktrackingEvent; + typedef typename Superclass::TEndBacktrackingEvent TEndBacktrackingEvent; + typedef typename Superclass::TBacktrackingEvent TBacktrackingEvent; + + typedef unsigned short TLabel; + typedef itk::Image< TLabel, I::ImageDimension > TLabelImage; + + protected: + typedef typename Superclass::_TVertices _TVertices; + typedef typename Superclass::_TCollision _TCollision; + typedef typename Superclass::_TCollisionsRow _TCollisionsRow; + typedef typename Superclass::_TCollisions _TCollisions; + typedef typename Superclass::_TNode _TNode; + typedef typename Superclass::_TNodes _TNodes; + + typedef typename I::PixelType _TPixel; + typedef typename I::RegionType _TRegion; + typedef typename I::SizeType _TSize; + + typedef std::pair< TResult, TVertex > _TCandidate; + typedef std::multimap< TResult, TVertex > _TCandidates; + + public: + itkNewMacro( Self ); + itkTypeMacro( DijkstraWithEndPointDetection, Dijkstra ); + + itkGetConstMacro( EndPoints, TVertices ); + itkGetConstMacro( BifurcationPoints, TVertices ); + itkGetConstMacro( NumberOfBranches, unsigned long ); + + public: + TLabelImage* GetLabelImage( ); + const TLabelImage* GetLabelImage( ) const; + void GraftLabelImage( itk::DataObject* obj ); + + protected: + DijkstraWithEndPointDetection( ); + virtual ~DijkstraWithEndPointDetection( ); + + virtual void _BeforeGenerateData( ); + virtual void _AfterGenerateData( ); + virtual void _SetResult( const TVertex& v, const _TNode& n ); + + _TRegion _Region( const TVertex& c, const double& r ); + + private: + // Purposely not implemented + DijkstraWithEndPointDetection( const Self& other ); + Self& operator=( const Self& other ); + + protected: + unsigned int m_LabelImageIndex; + + _TCandidates m_Candidates; + TVertices m_BifurcationPoints; + TVertices m_EndPoints; + unsigned long m_NumberOfBranches; + }; + + } // ecapseman + +} // ecapseman + +#include + +#endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__H__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/DijkstraWithEndPointDetection.hxx b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx new file mode 100644 index 0000000..2055cc4 --- /dev/null +++ b/lib/fpa/Image/DijkstraWithEndPointDetection.hxx @@ -0,0 +1,420 @@ +#ifndef __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__ +#define __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__ + +#include +#include + +// ------------------------------------------------------------------------- +template< class I, class O > +typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetLabelImage( ) +{ + return( + dynamic_cast< TLabelImage* >( + this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +const typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GetLabelImage( ) const +{ + return( + dynamic_cast< const TLabelImage* >( + this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +GraftLabelImage( itk::DataObject* obj ) +{ + TLabelImage* lbl = + dynamic_cast< TLabelImage* >( + this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex ) + ); + if( lbl != NULL ) + this->GraftNthOutput( this->m_LabelImageIndex, lbl ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +fpa::Image::DijkstraWithEndPointDetection< I, O >:: +DijkstraWithEndPointDetection( ) + : Superclass( ) +{ + this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( ); + this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 1 ); + this->itk::ProcessObject::SetNthOutput( + this->m_LabelImageIndex, TLabelImage::New( ) + ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +fpa::Image::DijkstraWithEndPointDetection< I, O >:: +~DijkstraWithEndPointDetection( ) +{ +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_BeforeGenerateData( ) +{ + // Correct seeds + /* TODO + 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::_BeforeGenerateData( ); + this->m_Candidates.clear( ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_AfterGenerateData( ) +{ + this->Superclass::_AfterGenerateData( ); + + // Finish base algorithm + /* TODO + 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( ) ); + this->InvokeEvent( TStartBacktrackingEvent( ) ); + + // 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 ); + TLabelImage* label = this->GetLabelImage( ); + label->FillBuffer( 0 ); + + // Prepare an object to hold marks + /* TODO + 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; + std::map< TLabel, std::pair< TVertex, TVertex > > branches; + for( ; cIt != this->m_Candidates.rend( ); ++cIt ) + { + // If pixel hasn't been visited, continue + TVertex v = cIt->second; + if( label->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( label->GetPixel( max_vertex ) != 0 ) + continue; + this->m_EndPoints.push_back( max_vertex ); + + // Get the path all the way to seed + std::vector< TVertex > path; + this->GetMinimumSpanningTree( )-> + GetPath( path, max_vertex, this->GetSeed( 0 ) ); + + // Mark branches + bool start = true; + bool change = false; + TVertex last_start = max_vertex; + typename std::vector< TVertex >::const_iterator pIt = path.begin( ); + for( ; pIt != path.end( ); ++pIt ) + { + this->InvokeEvent( + TBacktrackingEvent( *pIt, this->m_NumberOfBranches ) + ); + if( start ) + { + TLabel lbl = label->GetPixel( *pIt ); + if( lbl == 0 || lbl == this->m_NumberOfBranches ) + { + // Mark a sphere around current point as visited + double dist = std::sqrt( double( input->GetPixel( *pIt ) ) ); + region = this->_Region( max_vertex, dist * double( 1.1 ) ); + itk::ImageRegionIteratorWithIndex< TLabelImage > + lIt( label, region ); + for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt ) + { + if( lIt.Get( ) == 0 ) + lIt.Set( this->m_NumberOfBranches ); + + } // rof + + // Next vertex in current path + // TODO: 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! + branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt ); + std::cout << "bif " << this->m_NumberOfBranches << std::endl; + last_start = *pIt; + this->m_BifurcationPoints.push_back( *pIt ); + this->m_NumberOfBranches++; + start = false; + + } // fi + } + else + { + if( + std::find( + this->m_BifurcationPoints.begin( ), + this->m_BifurcationPoints.end( ), + *pIt + ) == 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.1 ) ); + itk::ImageRegionIteratorWithIndex< TLabelImage > + lIt( label, region ); + for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt ) + lIt.Set( this->m_NumberOfBranches ); + + // Next vertex in current path + /* TODO + 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 ); + branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt ); + std::cout << "change " << this->m_NumberOfBranches << std::endl; + + last_start = *pIt; + this->m_NumberOfBranches++; + /* TODO + this->m_FullTree[ max_vertex ] = + TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); + */ + // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl; + + } // fi + + } // fi + + } // rof + if( start || change ) + this->m_NumberOfBranches++; + this->InvokeEvent( TEndBacktrackingEvent( ) ); + + + + /* + this->InvokeEvent( TEndBacktrackingEvent( ) ); + */ + + } // rof + + std::cout << this->m_NumberOfBranches << " " << branches.size( ) << std::endl; + std::cout << this->m_EndPoints.size( ) << " " << this->m_BifurcationPoints.size( ) << std::endl; + + + // Re-enumerate labels + /* + std::map< TLabel, 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 O > +void fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_SetResult( const TVertex& v, const _TNode& n ) +{ + this->Superclass::_SetResult( v, n ); + + TResult vv = TResult( this->_VertexValue( v ) ); + if( TResult( 0 ) < vv ) + this->m_Candidates.insert( _TCandidate( n.Result / vv, v ) ); +} + +// ------------------------------------------------------------------------- +template< class I, class O > +typename fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >:: +_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__DIJKSTRAWITHENDPOINTDETECTION__HXX__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/MinimumSpanningTree.h b/lib/fpa/Image/MinimumSpanningTree.h deleted file mode 100644 index 652547d..0000000 --- a/lib/fpa/Image/MinimumSpanningTree.h +++ /dev/null @@ -1,54 +0,0 @@ -#ifndef __FPA__IMAGE__MINIMUMSPANNINGTREE__H__ -#define __FPA__IMAGE__MINIMUMSPANNINGTREE__H__ - -#include -#include - -namespace fpa -{ - namespace Image - { - /** - */ - template< class V, class N, class C, unsigned int D, unsigned long L > - class MinimumSpanningTree - : public fpa::Base::MinimumSpanningTree< V, C, itk::Image< N, D > > - { - public: - typedef fpa::Base::MinimumSpanningTree< V, C, itk::Image< N, D > > Superclass; - typedef MinimumSpanningTree Self; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - - typedef V TVertex; - typedef C TCollisions; - - protected: - typedef N _TNode; - - public: - itkNewMacro( Self ); - itkTypeMacro( MinimumSpanningTree, fpa_Base_MinimumSpanningTree ); - - protected: - MinimumSpanningTree( ); - virtual ~MinimumSpanningTree( ); - - virtual long _FrontId( const V& v ) const; - virtual V _Parent( const V& v ) const; - - private: - // Purposely not implemented - MinimumSpanningTree( const Self& other ); - Self& operator=( const Self& other ); - }; - - } // ecapseman - -} // ecapseman - -#include - -#endif // __FPA__IMAGE__MINIMUMSPANNINGTREE__H__ - -// eof - $RCSfile$ diff --git a/lib/fpa/Image/MinimumSpanningTree.hxx b/lib/fpa/Image/MinimumSpanningTree.hxx deleted file mode 100644 index 48a43cb..0000000 --- a/lib/fpa/Image/MinimumSpanningTree.hxx +++ /dev/null @@ -1,41 +0,0 @@ -#ifndef __FPA__IMAGE__MINIMUMSPANNINGTREE__HXX__ -#define __FPA__IMAGE__MINIMUMSPANNINGTREE__HXX__ - -// ------------------------------------------------------------------------- -template< class V, class N, class C, unsigned int D, unsigned long L > -fpa::Image::MinimumSpanningTree< V, N, C, D, L >:: -MinimumSpanningTree( ) - : Superclass( ) -{ -} - -// ------------------------------------------------------------------------- -template< class V, class N, class C, unsigned int D, unsigned long L > -fpa::Image::MinimumSpanningTree< V, N, C, D, L >:: -~MinimumSpanningTree( ) -{ -} - -// ------------------------------------------------------------------------- -template< class V, class N, class C, unsigned int D, unsigned long L > -long fpa::Image::MinimumSpanningTree< V, N, C, D, L >:: -_FrontId( const V& v ) const -{ - return( this->GetPixel( v ).FrontId ); -} - -// ------------------------------------------------------------------------- -template< class V, class N, class C, unsigned int D, unsigned long L > -V fpa::Image::MinimumSpanningTree< V, N, C, D, L >:: -_Parent( const V& v ) const -{ - _TNode n = this->GetPixel( v ); - if( n.Label == L ) - return( n.Parent ); - else - return( v ); -} - -#endif // __FPA__IMAGE__MINIMUMSPANNINGTREE__HXX__ - -// eof - $RCSfile$ diff --git a/lib/fpa/Image/RegionGrow.h b/lib/fpa/Image/RegionGrow.h index e20f062..71db89c 100644 --- a/lib/fpa/Image/RegionGrow.h +++ b/lib/fpa/Image/RegionGrow.h @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -16,10 +17,10 @@ namespace fpa */ template< class I, class O > class RegionGrow - : public Algorithm< I, O, fpa::Base::RegionGrow< typename I::IndexType, typename I::PixelType, typename O::PixelType, itk::ImageToImageFilter< I, O > > > + : public Algorithm< I, O, fpa::Base::RegionGrow< typename I::IndexType, typename I::PixelType, typename O::PixelType, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, itk::ImageToImageFilter< I, O > > > { public: - typedef fpa::Base::RegionGrow< typename I::IndexType, typename I::PixelType, typename O::PixelType, itk::ImageToImageFilter< I, O > > TBaseAlgorithm; + typedef fpa::Base::RegionGrow< typename I::IndexType, typename I::PixelType, typename O::PixelType, itk::Functor::IndexLexicographicCompare< I::ImageDimension >, itk::ImageToImageFilter< I, O > > TBaseAlgorithm; typedef RegionGrow Self; typedef Algorithm< I, O, TBaseAlgorithm > Superclass; diff --git a/lib/fpa/VTK/Image2DObserver.hxx b/lib/fpa/VTK/Image2DObserver.hxx index 74e5f4c..fd6247f 100644 --- a/lib/fpa/VTK/Image2DObserver.hxx +++ b/lib/fpa/VTK/Image2DObserver.hxx @@ -15,57 +15,6 @@ void fpa::VTK::Image2DObserver< F, R >:: SetRenderWindow( R* rw ) { this->m_RenderWindow = rw; - - /* TODO - this->m_Image = img; - unsigned int minD = TImage::ImageDimension; - minD = ( minD < 3 )? minD: 3; - - int e[ 6 ] = { 0 }; - typename TImage::RegionType reg = this->m_Image->GetRequestedRegion( ); - for( unsigned int i = 0; i < minD; i++ ) - { - e[ ( i << 1 ) + 0 ] = reg.GetIndex( )[ i ]; - e[ ( i << 1 ) + 1 ] = reg.GetIndex( )[ i ] + reg.GetSize( )[ i ] - 1; - - } // rof - - typename TImage::SpacingType spac = this->m_Image->GetSpacing( ); - double s[ 3 ] = { 1, 1, 1 }; - for( unsigned int i = 0; i < minD; i++ ) - s[ i ] = double( spac[ i ] ); - - typename TImage::PointType orig = this->m_Image->GetOrigin( ); - double o[ 3 ] = { 0 }; - for( unsigned int i = 0; i < minD; i++ ) - o[ i ] = double( orig[ i ] ); - - this->m_Stencil->SetExtent( e ); - this->m_Stencil->SetSpacing( s ); - this->m_Stencil->SetOrigin( o ); - this->m_Stencil->AllocateScalars( VTK_UNSIGNED_CHAR, 4 ); - for( unsigned int i = 0; i < 3; i++ ) - this->m_Stencil->GetPointData( )-> - GetScalars( )->FillComponent( i, 255 ); - this->m_Stencil->GetPointData( )->GetScalars( )->FillComponent( 3, 0 ); - - this->m_StencilActor->SetInputData( this->m_Stencil ); - this->m_StencilActor->InterpolateOff( ); - - double nPix = - double( this->m_Image->GetRequestedRegion( ).GetNumberOfPixels( ) ); - this->m_Percentage = ( unsigned long )( nPix * 0.01 ); - this->m_Number = 0; - - if( this->m_RenderWindow != NULL ) - { - vtkRenderer* ren = - this->m_RenderWindow->GetRenderers( )->GetFirstRenderer( ); - ren->AddActor( this->m_StencilActor ); - this->Render( ); - - } // fi - */ } // ------------------------------------------------------------------------- @@ -105,13 +54,13 @@ template< class F, class R > void fpa::VTK::Image2DObserver< F, R >:: Execute( const itk::Object* c, const itk::EventObject& e ) { - typedef typename F::TStartEvent TStartEvent; - typedef typename F::TStartLoopEvent TStartLoopEvent; - typedef typename F::TEndEvent TEndEvent; - typedef typename F::TEndLoopEvent TEndLoopEvent; - typedef typename F::TAliveEvent TAliveEvent; - typedef typename F::TFrontEvent TFrontEvent; - typedef typename F::TFreezeEvent TFreezeEvent; + typedef typename F::TStartEvent _TStartEvent; + typedef typename F::TStartLoopEvent _TStartLoopEvent; + typedef typename F::TEndEvent _TEndEvent; + typedef typename F::TEndLoopEvent _TEndLoopEvent; + typedef typename F::TAliveEvent _TAliveEvent; + typedef typename F::TFrontEvent _TFrontEvent; + typedef typename F::TFreezeEvent _TFreezeEvent; static unsigned char Colors[][4] = { @@ -137,7 +86,7 @@ Execute( const itk::Object* c, const itk::EventObject& e ) if( this->m_RenderWindow == NULL || filter == NULL ) return; - const TStartEvent* startEvt = dynamic_cast< const TStartEvent* >( &e ); + const _TStartEvent* startEvt = dynamic_cast< const _TStartEvent* >( &e ); if( startEvt != NULL ) { const typename F::TInputImage* img = filter->GetInput( ); @@ -185,8 +134,8 @@ Execute( const itk::Object* c, const itk::EventObject& e ) } // fi - const TAliveEvent* aliveEvt = dynamic_cast< const TAliveEvent* >( &e ); - const TFrontEvent* frontEvt = dynamic_cast< const TFrontEvent* >( &e ); + const _TAliveEvent* aliveEvt = dynamic_cast< const _TAliveEvent* >( &e ); + const _TFrontEvent* frontEvt = dynamic_cast< const _TFrontEvent* >( &e ); if( aliveEvt != NULL || frontEvt != NULL ) { if( aliveEvt != NULL ) @@ -212,7 +161,7 @@ Execute( const itk::Object* c, const itk::EventObject& e ) } // fi - const TEndEvent* endEvt = dynamic_cast< const TEndEvent* >( &e ); + const _TEndEvent* endEvt = dynamic_cast< const _TEndEvent* >( &e ); if( endEvt != NULL ) { vtkRenderer* ren = diff --git a/lib/fpa/VTK/Image3DObserver.hxx b/lib/fpa/VTK/Image3DObserver.hxx index 46ef743..721a89d 100644 --- a/lib/fpa/VTK/Image3DObserver.hxx +++ b/lib/fpa/VTK/Image3DObserver.hxx @@ -28,14 +28,18 @@ 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::TStartEvent _TStartEvent; + typedef typename F::TStartLoopEvent _TStartLoopEvent; + typedef typename F::TEndEvent _TEndEvent; + typedef typename F::TEndLoopEvent _TEndLoopEvent; + typedef typename F::TAliveEvent _TAliveEvent; + typedef typename F::TFrontEvent _TFrontEvent; + typedef typename F::TFreezeEvent _TFreezeEvent; + + typedef typename F::TStartBacktrackingEvent _TStartBacktrackingEvent; typedef typename F::TEndBacktrackingEvent _TEndBacktrackingEvent; + typedef typename F::TBacktrackingEvent _TBacktrackingEvent; // Check inputs if( this->m_RenderWindow == NULL ) @@ -51,122 +55,118 @@ Execute( const itk::Object* c, const itk::EventObject& e ) if( image == NULL ) return; - // Create actor - _TImage::RegionType reg = image->GetLargestPossibleRegion( ); - _TImage::SizeType siz = reg.GetSize( ); - if( this->m_Data.GetPointer( ) == NULL ) + const _TStartEvent* startEvt = dynamic_cast< const _TStartEvent* >( &e ); + const _TStartBacktrackingEvent* startBackEvt = + dynamic_cast< const _TStartBacktrackingEvent* >( &e ); + if( startEvt != NULL || startBackEvt != NULL ) { - this->m_Data = vtkSmartPointer< vtkPolyData >::New( ); - this->m_Mapper = vtkSmartPointer< vtkPolyDataMapper >::New( ); - this->m_Actor = vtkSmartPointer< vtkActor >::New( ); - - vtkSmartPointer< vtkPoints > points = - vtkSmartPointer< vtkPoints >::New( ); - vtkSmartPointer< vtkCellArray > vertices = - vtkSmartPointer< vtkCellArray >::New( ); - vtkSmartPointer< vtkFloatArray > scalars = - vtkSmartPointer< vtkFloatArray >::New( ); - this->m_Data->SetPoints( points ); - this->m_Data->SetVerts( vertices ); - this->m_Data->GetPointData( )->SetScalars( scalars ); - - this->m_Mapper->SetInputData( this->m_Data ); - this->m_Actor->SetMapper( this->m_Mapper ); - ren->AddActor( this->m_Actor ); - - this->m_Marks = TMarks::New( ); - this->m_Marks->SetLargestPossibleRegion( - image->GetLargestPossibleRegion( ) - ); - this->m_Marks->SetRequestedRegion( image->GetRequestedRegion( ) ); - this->m_Marks->SetBufferedRegion( image->GetBufferedRegion( ) ); - this->m_Marks->SetOrigin( image->GetOrigin( ) ); - this->m_Marks->SetSpacing( image->GetSpacing( ) ); - this->m_Marks->SetDirection( image->GetDirection( ) ); - this->m_Marks->Allocate( ); - this->m_Marks->FillBuffer( -1 ); - this->m_Count = 0; - this->m_RenderCount = reg.GetNumberOfPixels( ); + // Create actor + _TImage::RegionType reg = image->GetLargestPossibleRegion( ); + _TImage::SizeType siz = reg.GetSize( ); + if( this->m_Data.GetPointer( ) == NULL ) + { + this->m_Data = vtkSmartPointer< vtkPolyData >::New( ); + this->m_Mapper = vtkSmartPointer< vtkPolyDataMapper >::New( ); + this->m_Actor = vtkSmartPointer< vtkActor >::New( ); + + vtkSmartPointer< vtkPoints > points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > vertices = + vtkSmartPointer< vtkCellArray >::New( ); + vtkSmartPointer< vtkFloatArray > scalars = + vtkSmartPointer< vtkFloatArray >::New( ); + this->m_Data->SetPoints( points ); + this->m_Data->SetVerts( vertices ); + this->m_Data->GetPointData( )->SetScalars( scalars ); + + this->m_Mapper->SetInputData( this->m_Data ); + this->m_Actor->SetMapper( this->m_Mapper ); + ren->AddActor( this->m_Actor ); + + this->m_Marks = TMarks::New( ); + this->m_Marks->SetLargestPossibleRegion( + image->GetLargestPossibleRegion( ) + ); + this->m_Marks->SetRequestedRegion( image->GetRequestedRegion( ) ); + this->m_Marks->SetBufferedRegion( image->GetBufferedRegion( ) ); + this->m_Marks->SetOrigin( image->GetOrigin( ) ); + this->m_Marks->SetSpacing( image->GetSpacing( ) ); + this->m_Marks->SetDirection( image->GetDirection( ) ); + this->m_Marks->Allocate( ); + this->m_Marks->FillBuffer( -1 ); + this->m_Count = 0; + this->m_RenderCount = reg.GetNumberOfPixels( ); - } // fi + } // fi + return; - // Get possible events - const _TEvent* evt = dynamic_cast< const _TEvent* >( &e ); - _TFrontEvent fevt; - _TMarkEvent mevt; - _TCollisionEvent cevt; - _TEndEvent eevt; + } // fi - if( evt != NULL ) + const _TAliveEvent* aliveEvt = dynamic_cast< const _TAliveEvent* >( &e ); + const _TFrontEvent* frontEvt = dynamic_cast< const _TFrontEvent* >( &e ); + if( aliveEvt != NULL || frontEvt != NULL ) { - if( fevt.CheckEvent( evt ) ) - { - if( this->m_Marks->GetPixel( evt->Node.Vertex ) == -1 ) - { - typename _TImage::PointType pnt; - image->TransformIndexToPhysicalPoint( evt->Node.Vertex, pnt ); - - long pId = this->m_Data->GetNumberOfPoints( ); - this->m_Data->GetVerts( )->InsertNextCell( 1 ); - this->m_Data->GetVerts( )->InsertCellPoint( pId ); - this->m_Data->GetPoints( )-> - InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - this->m_Data->GetPointData( )-> - GetScalars( )->InsertNextTuple1( 0.5 ); - this->m_Data->Modified( ); - - this->m_Marks->SetPixel( evt->Node.Vertex, pId ); - this->m_Count++; - - // Render visual debug - double per = double( this->m_RenderCount ) * this->m_RenderPercentage; - if( double( this->m_Count ) >= per ) - this->m_RenderWindow->Render( ); - if( double( this->m_Count ) >= per ) - this->m_Count = 0; - - } // fi - return; - } - else if( mevt.CheckEvent( evt ) ) + _TImage::IndexType vertex; + if( aliveEvt != NULL ) + vertex = aliveEvt->Vertex; + else if( frontEvt != NULL ) + vertex = frontEvt->Vertex; + + if( this->m_Marks->GetPixel( vertex ) == -1 ) { - // TODO: std::cout << "mark" << std::endl; - return; + typename _TImage::PointType pnt; + image->TransformIndexToPhysicalPoint( vertex, pnt ); - } // fi + long pId = this->m_Data->GetNumberOfPoints( ); + this->m_Data->GetVerts( )->InsertNextCell( 1 ); + this->m_Data->GetVerts( )->InsertCellPoint( pId ); + this->m_Data->GetPoints( )-> + InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + this->m_Data->GetPointData( )-> + GetScalars( )->InsertNextTuple1( 0.5 ); + this->m_Data->Modified( ); - if( cevt.CheckEvent( evt ) ) - { - // TODO: std::cout << "collision" << std::endl; - return; + this->m_Marks->SetPixel( vertex, pId ); + this->m_Count++; + + // Render visual debug + double per = double( this->m_RenderCount ) * this->m_RenderPercentage; + if( double( this->m_Count ) >= per ) + this->m_RenderWindow->Render( ); + if( double( this->m_Count ) >= per ) + this->m_Count = 0; } // fi + return; - if( eevt.CheckEvent( evt ) ) - { - this->m_RenderWindow->Render( ); - ren->RemoveActor( this->m_Actor ); - this->m_RenderWindow->Render( ); - this->m_Marks = NULL; - this->m_Data = NULL; - this->m_Mapper = NULL; - this->m_Actor = NULL; - return; + } // fi - } // fi - } - else + const _TEndEvent* endEvt = dynamic_cast< const _TEndEvent* >( &e ); + if( endEvt != NULL ) + { + this->m_RenderWindow->Render( ); + ren->RemoveActor( this->m_Actor ); + this->m_RenderWindow->Render( ); + this->m_Marks = NULL; + this->m_Data = NULL; + this->m_Mapper = NULL; + this->m_Actor = NULL; + return; + + } // fi + + + const _TBacktrackingEvent* backEvt = + dynamic_cast< const _TBacktrackingEvent* >( &e ); + const _TEndBacktrackingEvent* endBackEvt = + dynamic_cast< const _TEndBacktrackingEvent* >( &e ); + if( backEvt != NULL ) { - 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; - double back_id = double( bevt->BackId % nColors ) / double( nColors ); + double back_id = + double( backEvt->FrontId % nColors ) / double( nColors ); typename _TImage::PointType pnt; - image->TransformIndexToPhysicalPoint( bevt->Node, pnt ); + image->TransformIndexToPhysicalPoint( backEvt->Vertex, pnt ); long pId = this->m_Data->GetNumberOfPoints( ); this->m_Data->GetVerts( )->InsertNextCell( 1 ); @@ -180,14 +180,17 @@ Execute( const itk::Object* c, const itk::EventObject& e ) } // fi - if( ebevt != NULL ) + if( endBackEvt != NULL ) { this->m_RenderWindow->Render( ); + /* TODO: DEBUG + std::cout << "Press enter: " << std::ends; + int aux; + std::cin >> aux; + */ return; } // fi - - } // fi } // ------------------------------------------------------------------------- -- 2.47.1