From: Leonardo Flórez-Valencia Date: Mon, 13 Mar 2017 22:52:29 +0000 (-0500) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=49336779a8037ecd49fc0c4f6038c02636eb538b;p=FrontAlgorithms.git ... --- diff --git a/CMakeLists.txt b/CMakeLists.txt index 369cba1..aab256e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,13 +33,13 @@ SET(prj_SHORT_VERSION "${prj_MAJOR_VERSION}") ## == Functions, packages and options == ## ===================================== -#INCLUDE(cmake/BaseConfig.cmake) +INCLUDE(cmake/BaseConfig.cmake) +INCLUDE(cmake/KitwareTools.cmake) #OPTION(USE_cpPlugins "Build cpPlugins-based code" OFF) #IF(USE_cpPlugins) # FIND_PACKAGE(cpPlugins) # MARK_AS_ADVANCED(CLEAR cpPlugins_DIR) #ENDIF(USE_cpPlugins) -#INCLUDE(cmake/KitwareTools.cmake) #INCLUDE(cmake/QtTools.cmake) #INCLUDE(cmake/Functions.cmake) diff --git a/cmake/BaseConfig.cmake b/cmake/BaseConfig.cmake new file mode 100644 index 0000000..d177604 --- /dev/null +++ b/cmake/BaseConfig.cmake @@ -0,0 +1,51 @@ +## ======================================================================= +## == Force c++11 language version == +## == NOTE: It seems that by default on Visual Studio Compiler supports == +## == c++11, so it only need to be tested on other OS. == +## ======================================================================= + +IF(NOT MSVC) + INCLUDE(CheckCXXCompilerFlag) + CHECK_CXX_COMPILER_FLAG("-std=c++11" COMPILER_SUPPORTS_CXX11) + IF(COMPILER_SUPPORTS_CXX11) + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") + ELSE(COMPILER_SUPPORTS_CXX11) + CHECK_CXX_COMPILER_FLAG("-std=c++0x" COMPILER_SUPPORTS_CXX0X) + IF(COMPILER_SUPPORTS_CXX0X) + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x") + ELSE(COMPILER_SUPPORTS_CXX0X) + MESSAGE( + FATAL_ERROR + "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support." + ) + ENDIF(COMPILER_SUPPORTS_CXX0X) + ENDIF(COMPILER_SUPPORTS_CXX11) +ENDIF(NOT MSVC) + +## =================================================== +## == Prepare header generator to build shared libs == +## =================================================== + +INCLUDE(GenerateExportHeader) + +## ================================================== +## == Do not allow to build inside the source tree == +## ================================================== + +IF(PROJECT_BINARY_DIR STREQUAL ${PROJECT_SOURCE_DIR}) + MESSAGE(FATAL_ERROR "Building in the source tree is not allowed.") +ENDIF(PROJECT_BINARY_DIR STREQUAL ${PROJECT_SOURCE_DIR}) + +## ================================================= +## == Where to put targets (executables and libs) == +## ================================================= + +SET(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}) +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}) +MARK_AS_ADVANCED( + CMAKE_BACKWARDS_COMPATIBILITY + EXECUTABLE_OUTPUT_PATH + LIBRARY_OUTPUT_PATH + ) + +## eof - $RCSfile$ diff --git a/cmake/KitwareTools.cmake b/cmake/KitwareTools.cmake new file mode 100644 index 0000000..35ae868 --- /dev/null +++ b/cmake/KitwareTools.cmake @@ -0,0 +1,8 @@ +# ======================== +# == Find Kitware tools == +# ======================== + +FIND_PACKAGE(ITK REQUIRED) +INCLUDE(${ITK_USE_FILE}) + +## eof - $RCSfile$ diff --git a/examples/BronchiiInitialSegmentationWithMori.cxx b/examples/BronchiiInitialSegmentationWithMori.cxx new file mode 100644 index 0000000..ed0c294 --- /dev/null +++ b/examples/BronchiiInitialSegmentationWithMori.cxx @@ -0,0 +1,104 @@ +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +// ------------------------------------------------------------------------- +static const unsigned int Dim = 3; +typedef short TPixel; +typedef itk::Image< TPixel, Dim > TImage; + +// ------------------------------------------------------------------------- +std::string GetFullPath( const std::string& filename ) +{ + /* On windows: + #include + TCHAR full_path[MAX_PATH]; + GetFullPathName(_T("foo.dat"), MAX_PATH, full_path, NULL); + */ + + char* buffer = realpath( filename.c_str( ), NULL ); + std::string path = buffer; + std::free( buffer ); + return( path ); +} + +// ------------------------------------------------------------------------- +std::string GetFullPathToDirectory( const std::string& filename ) +{ + std::string path = GetFullPath( filename ); + std::size_t found = path.find_last_of( "/\\" ); + return( path.substr( 0, found + 1 ) ); +} + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + // Command configuration + if( argc < 3 ) + { + std::cerr + << "Usage: " << argv[ 0 ] << " input_image output_image" + << std::endl; + return( 1 ); + + } // fi + std::string input_image_filename = GetFullPath( argv[ 1 ] ); + std::string output_image_filename = argv[ 2 ]; + + // Try to guess initial seed + std::string seed_filename = GetFullPathToDirectory( argv[ 1 ] ) + "seed.txt"; + std::ifstream seed_str( seed_filename.c_str( ) ); + if( !seed_str ) + { + std::cerr + << "No \"seed.txt\" file in the same input image directory." + << std::endl; + return( 1 ); + + } // fi + TImage::IndexType seed; + seed_str >> seed[ 0 ] >> seed[ 1 ] >> seed[ 2 ]; + + // Read image + typedef itk::ImageFileReader< TImage > TReader; + TReader::Pointer reader = TReader::New( ); + reader->SetFileName( input_image_filename ); + + // Mori's algorithms + typedef fpa::Image::MoriFilter< TImage, TImage > TFilter; + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( reader->GetOutput( ) ); + filter->SetSeed( seed ); + filter->SetThresholdRange( -1024, 0 ); + filter->SetInsideValue( 1 ); + filter->SetOutsideValue( 0 ); + + // Write results + typedef itk::ImageFileWriter< TImage > TWriter; + TWriter::Pointer writer = TWriter::New( ); + writer->SetInput( filter->GetOutput( ) ); + writer->SetFileName( output_image_filename ); + + // Execute pipeline + try + { + writer->Update( ); + } + catch( std::exception& err ) + { + std::cerr << "Error caught: " << err.what( ) << std::endl; + return( 1 ); + + } // yrt + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt new file mode 100644 index 0000000..ddaadf3 --- /dev/null +++ b/examples/CMakeLists.txt @@ -0,0 +1,13 @@ +OPTION(BUILD_EXAMPLES "Build cpPlugins-free examples." OFF) +IF(BUILD_EXAMPLES) + SET( + _examples + BronchiiInitialSegmentationWithMori + ) + FOREACH(_e ${_examples}) + ADD_EXECUTABLE(fpa_example_${_e} ${_e}.cxx) + TARGET_LINK_LIBRARIES(fpa_example_${_e} ${ITK_LIBRARIES}) + ENDFOREACH(_e) +ENDIF(BUILD_EXAMPLES) + +## eof - $RCSfile$ diff --git a/libs/fpa/Base/Algorithm.h b/libs/fpa/Base/Algorithm.h new file mode 100644 index 0000000..6625467 --- /dev/null +++ b/libs/fpa/Base/Algorithm.h @@ -0,0 +1,150 @@ +#ifndef __fpa__Base__Algorithm__h__ +#define __fpa__Base__Algorithm__h__ + +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TFilter, class _TVertex, class _TOutput > + class Algorithm + : public _TFilter + { + public: + typedef Algorithm Self; + typedef _TFilter Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TVertex TVertex; + typedef _TOutput TOutput; + typedef unsigned long TFrontId; + + // Different functions + typedef std::vector< TVertex > TNeighborhood; + typedef itk::FunctionBase< TVertex, TNeighborhood > TNeighborhoodFunction; + typedef itk::FunctionBase< TOutput, TOutput > TConversionFunction; + typedef fpa::Base::Functors::VertexCostFunctionBase< TVertex, TOutput > TVertexFunction; + + // Minigraph to represent collisions + typedef std::pair< _TVertex, bool > TCollision; + typedef std::vector< TCollision > TCollisionsRow; + typedef std::vector< TCollisionsRow > TCollisions; + + protected: + struct _TQueueNode + { + TVertex Vertex; + TVertex Parent; + TOutput Result; + TFrontId FrontId; + _TQueueNode( ); + _TQueueNode( const TVertex& v ); + _TQueueNode( const TVertex& v, const _TQueueNode& n ); + }; + + public: + itkTypeMacro( Self, _TFilter ); + + fpa_Base_NewEvent( TStartEvent ); + fpa_Base_NewEvent( TEndEvent ); + fpa_Base_NewEvent( TStartLoopEvent ); + fpa_Base_NewEvent( TEndLoopEvent ); + fpa_Base_NewEventWithVertex( TPushEvent, TVertex ); + fpa_Base_NewEventWithVertex( TPopEvent, TVertex ); + fpa_Base_NewEventWithVertex( TMarkEvent, TVertex ); + + itkGetConstMacro( InitResult, _TOutput ); + itkSetMacro( InitResult, _TOutput ); + + itkBooleanMacro( StopAtOneFront ); + itkGetConstMacro( StopAtOneFront, bool ); + itkSetMacro( StopAtOneFront, bool ); + + itkGetObjectMacro( NeighborhoodFunction, TNeighborhoodFunction ); + itkGetObjectMacro( ConversionFunction, TConversionFunction ); + itkGetObjectMacro( VertexFunction, TVertexFunction ); + + itkGetConstObjectMacro( NeighborhoodFunction, TNeighborhoodFunction ); + itkGetConstObjectMacro( ConversionFunction, TConversionFunction ); + itkGetConstObjectMacro( VertexFunction, TVertexFunction ); + + itkSetObjectMacro( NeighborhoodFunction, TNeighborhoodFunction ); + itkSetObjectMacro( ConversionFunction, TConversionFunction ); + itkSetObjectMacro( VertexFunction, TVertexFunction ); + + public: + void ClearSeeds( ); + void AddSeed( const TVertex& seed, const TOutput& value ); + + protected: + Algorithm( ); + virtual ~Algorithm( ); + + virtual void GenerateData( ) override; + + // Particular methods + virtual bool _ContinueGenerateData( ); + virtual void _Loop( ); + + virtual void _BeforeGenerateData( ); + virtual void _AfterGenerateData( ); + virtual void _BeforeLoop( ); + virtual void _AfterLoop( ); + + virtual bool _ValidLoop( ) const; + virtual void _UpdateCollisions( const TVertex& a, const TVertex& b ); + + virtual _TOutput _GetInputValue( const _TQueueNode& v, const _TQueueNode& p ); + + virtual void _InitMarks( ) = 0; + virtual void _InitResults( const TOutput& init_value ) = 0; + virtual bool _IsMarked( const _TVertex& v ) const = 0; + virtual void _Mark( const _TQueueNode& n ) = 0; + virtual TFrontId _GetMark( const _TVertex& v ) const = 0; + virtual void _UpdateResult( const _TQueueNode& n ) = 0; + virtual TOutput _GetResult( const _TVertex& v ) const = 0; + virtual unsigned int _GetNumberOfDimensions( ) const = 0; + + virtual bool _UpdateValue( _TQueueNode& v, const _TQueueNode& p ) = 0; + virtual unsigned long _QueueSize( ) const = 0; + virtual void _QueueClear( ) = 0; + virtual void _QueuePush( const _TQueueNode& node ) = 0; + virtual _TQueueNode _QueuePop( ) = 0; + + private: + // Purposely not implemented. + Algorithm( const Self& other ); + Self& operator=( const Self& other ); + + protected: + _TOutput m_InitResult; + bool m_StopAtOneFront; + + typename TNeighborhoodFunction::Pointer m_NeighborhoodFunction; + typename TConversionFunction::Pointer m_ConversionFunction; + typename TVertexFunction::Pointer m_VertexFunction; + + std::vector< _TQueueNode > m_Seeds; + TCollisions m_Collisions; + unsigned int m_NumberOfFronts; + }; + + } // ecaseman + +} // ecaseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__Algorithm__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Algorithm.hxx b/libs/fpa/Base/Algorithm.hxx new file mode 100644 index 0000000..c29ff03 --- /dev/null +++ b/libs/fpa/Base/Algorithm.hxx @@ -0,0 +1,269 @@ +#ifndef __fpa__Base__Algorithm__hxx__ +#define __fpa__Base__Algorithm__hxx__ + +#include + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >::_TQueueNode:: +_TQueueNode( ) +{ + this->Vertex.Fill( 0 ); + this->Parent.Fill( 0 ); + this->Result = _TOutput( 0 ); + this->FrontId = 0; +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >::_TQueueNode:: +_TQueueNode( const _TVertex& v ) +{ + this->Vertex = v; + this->Parent = v; + this->Result = _TOutput( 0 ); + this->FrontId = 0; +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >::_TQueueNode:: +_TQueueNode( const _TVertex& v, const _TQueueNode& n ) +{ + this->Vertex = v; + this->Parent = n.Vertex; + this->Result = n.Result; + this->FrontId = n.FrontId; +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +ClearSeeds( ) +{ + this->m_Seeds.clear( ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +AddSeed( const _TVertex& seed, const TOutput& value ) +{ + _TQueueNode node( seed ); + node.FrontId = this->m_Seeds.size( ) + 1; + node.Result = value; + this->m_Seeds.push_back( node ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +Algorithm( ) + : Superclass( ), + m_InitResult( _TOutput( 0 ) ), + m_StopAtOneFront( false ) +{ +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +~Algorithm( ) +{ +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +GenerateData( ) +{ + this->InvokeEvent( TStartEvent( ) ); + this->_BeforeGenerateData( ); + this->m_NumberOfFronts = this->m_Seeds.size( ); + if( this->m_NumberOfFronts > 0 ) + { + // Prepare collisions + TCollision coll( TVertex( ), false ); + TCollisionsRow row( this->m_NumberOfFronts, coll ); + this->m_Collisions.clear( ); + this->m_Collisions.resize( this->m_NumberOfFronts, row ); + + // Put seeds on queue + this->_QueueClear( ); + for( + auto nIt = this->m_Seeds.begin( ); + nIt != this->m_Seeds.end( ); + ++nIt + ) + this->_QueuePush( *nIt ); + + // Init marks and results + this->_InitMarks( ); + this->_InitResults( this->m_InitResult ); + + // Main loop + do + this->_Loop( ); + while( this->_ContinueGenerateData( ) ); + + } // fi + this->_AfterGenerateData( ); + this->InvokeEvent( TEndEvent( ) ); +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +bool fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_ContinueGenerateData( ) +{ + return( false ); +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_Loop( ) +{ + this->InvokeEvent( TStartLoopEvent( ) ); + this->_BeforeLoop( ); + while( this->_ValidLoop( ) ) + { + _TQueueNode node = this->_QueuePop( ); + this->InvokeEvent( TPopEvent( node.Vertex, node.FrontId ) ); + + // Process actual vertex + if( this->_IsMarked( node.Vertex ) ) + continue; + this->_Mark( node ); + this->_UpdateResult( node ); + this->InvokeEvent( TMarkEvent( node.Vertex, node.FrontId ) ); + + // Add neighbors to queue + TNeighborhood neighs = this->m_NeighborhoodFunction->Evaluate( node.Vertex ); + for( auto it = neighs.begin( ); it != neighs.end( ); ++it ) + { + if( !( this->_IsMarked( *it ) ) ) + { + _TQueueNode neigh( *it, node ); + if( this->_UpdateValue( neigh, node ) ) + { + this->_QueuePush( neigh ); + this->InvokeEvent( TPushEvent( node.Vertex, node.FrontId ) ); + + } // fi + } + else + this->_UpdateCollisions( node.Vertex, *it ); + + } // rof + + } // elihw + this->_AfterLoop( ); + this->InvokeEvent( TEndLoopEvent( ) ); +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_BeforeGenerateData( ) +{ +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_AfterGenerateData( ) +{ +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_BeforeLoop( ) +{ +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_AfterLoop( ) +{ +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +bool fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_ValidLoop( ) const +{ + bool r = ( this->_QueueSize( ) > 0 ); + if( this->m_StopAtOneFront ) + r &= ( this->m_NumberOfFronts > 0 ); + return( r ); +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +void fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_UpdateCollisions( const TVertex& a, const TVertex& b ) +{ + auto ma = this->_GetMark( a ); + auto mb = this->_GetMark( b ); + if( ma == mb || ma == 0 || mb == 0 ) + return; + + // Mark collision, if it is new + ma--; mb--; + bool ret = false; + bool exists = this->m_Collisions[ ma ][ mb ].second; + exists &= this->m_Collisions[ mb ][ ma ].second; + if( !exists ) + { + this->m_Collisions[ ma ][ mb ].first = a; + this->m_Collisions[ ma ][ mb ].second = true; + this->m_Collisions[ mb ][ ma ].first = b; + this->m_Collisions[ mb ][ ma ].second = true; + + // Update number of fronts + unsigned long N = this->m_Seeds.size( ); + unsigned long count = 0; + std::vector< bool > m( N, false ); + std::queue< unsigned long > q; + q.push( 0 ); + while( !q.empty( ) ) + { + unsigned long f = q.front( ); + q.pop( ); + + if( m[ f ] ) + continue; + m[ f ] = true; + count++; + + for( unsigned int n = 0; n < N; ++n ) + if( this->m_Collisions[ f ][ n ].second && !m[ n ] ) + q.push( n ); + + } // elihw + this->m_NumberOfFronts = N - count; + + } // fi +} + +// ------------------------------------------------------------------------- +template < class _TFilter, class _TVertex, class _TOutput > +_TOutput fpa::Base::Algorithm< _TFilter, _TVertex, _TOutput >:: +_GetInputValue( const _TQueueNode& v, const _TQueueNode& p ) +{ + _TOutput res = this->m_InitResult; + if( this->m_VertexFunction.IsNotNull( ) ) + res = this->m_VertexFunction->Evaluate( v.Vertex, p.Vertex ); + if( this->m_ConversionFunction.IsNotNull( ) ) + res = this->m_ConversionFunction->Evaluate( res ); + return( res ); +} + +#endif // __fpa__Base__Algorithm__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Dijkstra.h b/libs/fpa/Base/Dijkstra.h new file mode 100644 index 0000000..08c704e --- /dev/null +++ b/libs/fpa/Base/Dijkstra.h @@ -0,0 +1,66 @@ +#ifndef __fpa__Base__Dijkstra__h__ +#define __fpa__Base__Dijkstra__h__ + +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TSuperclass, class _TMST > + class Dijkstra + : public fpa::Base::PriorityQueueAlgorithm< _TSuperclass > + { + public: + typedef Dijkstra Self; + typedef fpa::Base::PriorityQueueAlgorithm< _TSuperclass > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TMST TMST; + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TVertex TVertex; + + protected: + typedef typename Superclass::_TQueueNode _TQueueNode; + + public: + itkTypeMacro( Dijkstra, Algorithm ); + + public: + _TMST* GetMinimumSpanningTree( ); + const _TMST* GetMinimumSpanningTree( ) const; + + protected: + Dijkstra( ); + virtual ~Dijkstra( ); + + virtual void _AfterGenerateData( ) override; + + virtual void _UpdateResult( const _TQueueNode& n ) override; + virtual bool _UpdateValue( + _TQueueNode& v, const _TQueueNode& p + ) override; + + private: + // Purposely not defined + Dijkstra( const Self& other ); + Self& operator=( const Self& other ); + + protected: + unsigned long m_MSTIndex; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__Dijkstra__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Dijkstra.hxx b/libs/fpa/Base/Dijkstra.hxx new file mode 100644 index 0000000..df9c75e --- /dev/null +++ b/libs/fpa/Base/Dijkstra.hxx @@ -0,0 +1,94 @@ +#ifndef __fpa__Base__Dijkstra__hxx__ +#define __fpa__Base__Dijkstra__hxx__ + +#include +#include + +// ------------------------------------------------------------------------- +template< class _TSuperclass, class _TMST > +_TMST* fpa::Base::Dijkstra< _TSuperclass, _TMST >:: +GetMinimumSpanningTree( ) +{ + return( + dynamic_cast< _TMST* >( + this->itk::ProcessObject::GetOutput( this->m_MSTIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass, class _TMST > +const _TMST* fpa::Base::Dijkstra< _TSuperclass, _TMST >:: +GetMinimumSpanningTree( ) const +{ + return( + dynamic_cast< const _TMST* >( + this->itk::ProcessObject::GetOutput( this->m_MSTIndex ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass, class _TMST > +fpa::Base::Dijkstra< _TSuperclass, _TMST >:: +Dijkstra( ) + : Superclass( ) +{ + this->m_InitResult = TOutput( 0 ); + this->m_MSTIndex = this->GetNumberOfRequiredOutputs( ); + this->SetNumberOfRequiredOutputs( this->m_MSTIndex + 1 ); + this->itk::ProcessObject::SetNthOutput( this->m_MSTIndex, _TMST::New( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass, class _TMST > +fpa::Base::Dijkstra< _TSuperclass, _TMST >:: +~Dijkstra( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass, class _TMST > +void fpa::Base::Dijkstra< _TSuperclass, _TMST >:: +_AfterGenerateData( ) +{ + this->Superclass::_AfterGenerateData( ); + + auto mst = this->GetMinimumSpanningTree( ); + mst->ClearSeeds( ); + for( auto s = this->m_Seeds.begin( ); s != this->m_Seeds.end( ); ++s ) + mst->AddSeed( s->Vertex ); + mst->SetCollisions( this->m_Collisions ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass, class _TMST > +void fpa::Base::Dijkstra< _TSuperclass, _TMST >:: +_UpdateResult( const _TQueueNode& n ) +{ + this->Superclass::_UpdateResult( n ); + this->GetMinimumSpanningTree( )->SetParent( n.Vertex, n.Parent ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass, class _TMST > +bool fpa::Base::Dijkstra< _TSuperclass, _TMST >:: +_UpdateValue( _TQueueNode& v, const _TQueueNode& p ) +{ + v.Result = this->_GetInputValue( p.Vertex, v.Vertex ); + if( v.Result >= TOutput( 0 ) ) + { + v.Result += p.Result; + return( true ); + } + else + { + v.Result = this->m_InitResult; + return( false ); + + } // fi +} + +#endif // __fpa__Base__Dijkstra__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Events.h b/libs/fpa/Base/Events.h new file mode 100644 index 0000000..17c5564 --- /dev/null +++ b/libs/fpa/Base/Events.h @@ -0,0 +1,100 @@ +#ifndef __fpa__Base__Events__h__ +#define __fpa__Base__Events__h__ + +// ------------------------------------------------------------------------- +#define fpa_Base_NewEvent( name ) \ + class name \ + : public BaseEvent \ + { \ + public: \ + name( ) : BaseEvent( ) { } \ + virtual ~name( ) { } \ + const char* GetEventName( ) const \ + { return( "fpa::Base::##name" ); } \ + bool CheckEvent( const itk::EventObject* e ) const \ + { return( dynamic_cast< const name* >( e ) != NULL ); } \ + itk::EventObject* MakeObject( ) const \ + { return( new name( ) ); } \ + }; + +// ------------------------------------------------------------------------- +#define fpa_Base_NewEventWithVertex( name, type ) \ + class name \ + : public BaseEventWithVertex< type > \ + { \ + public: \ + name( ) : BaseEventWithVertex< type >( ) { } \ + name( const type& v, long fid ) \ + : BaseEventWithVertex< type >( v, fid ) { } \ + virtual ~name( ) { } \ + const char* GetEventName( ) const \ + { return( "fpa::Base::##name" ); } \ + bool CheckEvent( const itk::EventObject* e ) const \ + { return( dynamic_cast< const name* >( e ) != NULL ); } \ + itk::EventObject* MakeObject( ) const \ + { return( new name( ) ); } \ + }; + +namespace fpa +{ + namespace Base + { + /** + */ + class BaseEvent + : public itk::AnyEvent + { + public: + BaseEvent( ) + : itk::AnyEvent( ) + { } + virtual ~BaseEvent( ) + { } + const char* GetEventName( ) const + { return( "fpa::Base::BaseEvent" ); } + bool CheckEvent( const itk::EventObject* e ) const + { return( dynamic_cast< const BaseEvent* >( e ) != NULL ); } + itk::EventObject* MakeObject( ) const + { return( new BaseEvent( ) ); } + }; + + /** + */ + template< class _TVertex > + class BaseEventWithVertex + : public BaseEvent + { + public: + BaseEventWithVertex( ) + : BaseEvent( ) + { } + BaseEventWithVertex( const _TVertex& v, long fid ) + : BaseEvent( ), + Vertex( v ), + FrontId( fid ) + { } + virtual ~BaseEventWithVertex( ) + { } + const char* GetEventName( ) const + { return( "fpa::Base::BaseEventWithVertex" ); } + bool CheckEvent( const itk::EventObject* e ) const + { + return( + dynamic_cast< const BaseEventWithVertex< _TVertex >* >( e ) != NULL + ); + } + itk::EventObject* MakeObject( ) const + { return( new BaseEventWithVertex< _TVertex >( ) ); } + + public: + _TVertex Vertex; + long FrontId; + }; + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Base__Events__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/FastMarching.h b/libs/fpa/Base/FastMarching.h new file mode 100644 index 0000000..acea8e7 --- /dev/null +++ b/libs/fpa/Base/FastMarching.h @@ -0,0 +1,59 @@ +#ifndef __fpa__Base__FastMarching__h__ +#define __fpa__Base__FastMarching__h__ + +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TSuperclass > + class FastMarching + : public fpa::Base::PriorityQueueAlgorithm< _TSuperclass > + { + public: + typedef FastMarching Self; + typedef fpa::Base::PriorityQueueAlgorithm< _TSuperclass > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TFrontId TFrontId; + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TVertex TVertex; + + typedef std::pair< TVertex, double > TFastMarchingNeighbor; + typedef std::vector< TFastMarchingNeighbor > TFastMarchingNeighborhood; + + protected: + typedef typename Superclass::_TQueueNode _TQueueNode; + + public: + itkTypeMacro( FastMarching, Algorithm ); + + protected: + FastMarching( ); + virtual ~FastMarching( ); + + virtual TFastMarchingNeighborhood _FastMarchingNeighbors( const TVertex& v ) const = 0; + + virtual bool _UpdateValue( _TQueueNode& v, const _TQueueNode& p ) override; + + private: + // Purposely not defined + FastMarching( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__FastMarching__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/FastMarching.hxx b/libs/fpa/Base/FastMarching.hxx new file mode 100644 index 0000000..0fe1cb2 --- /dev/null +++ b/libs/fpa/Base/FastMarching.hxx @@ -0,0 +1,82 @@ +#ifndef __fpa__Base__FastMarching__hxx__ +#define __fpa__Base__FastMarching__hxx__ + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +fpa::Base::FastMarching< _TSuperclass >:: +FastMarching( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +fpa::Base::FastMarching< _TSuperclass >:: +~FastMarching( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +bool fpa::Base::FastMarching< _TSuperclass >:: +_UpdateValue( _TQueueNode& v, const _TQueueNode& p ) +{ + // Compute quadratic coefficients + double a = double( 0 ); + double b = double( 0 ); + double c = double( 0 ); + TFastMarchingNeighborhood neighs = this->_FastMarchingNeighbors( v.Vertex ); + for( unsigned int i = 0; i < neighs.size( ); i += 2 ) + { + double tn = double( this->_GetResult( neighs[ i ].first ) ); + double tp = double( this->_GetResult( neighs[ i + 1 ].first ) ); + double dd = double( 0 ); + double td = this->m_InitResult; + if( tn < tp ) + { + td = tn; + dd = neighs[ i ].second; + } + else + { + td = tp; + dd = neighs[ i + 1 ].second; + + } // fi + + if( td < double( this->m_InitResult ) ) + { + dd = double( 1 ) / ( dd * dd ); + a += dd; + b += td * dd; + c += td * td * dd; + + } // fi + + } // rof + double F = double( this->_GetInputValue( v, p ) ); + c -= double( 1 ) / ( F * F ); + + // Solve quadratic equation + double d = ( b * b ) - ( a * c ); + if( d >= double( 0 ) ) + { + d = std::sqrt( d ); + b *= double( -1 ); + double s1 = std::fabs( ( b + d ) / a ); + double s2 = std::fabs( ( b - d ) / a ); + v.Result = TOutput( ( s2 < s1 )? s1: s2 ); + if( v.Result < this->_GetResult( v.Vertex ) ) + this->_UpdateResult( v ); + } + else + { + std::cout << std::endl << "-- FM --> " << v.Vertex << " " << d << std::endl; + std::exit( 1 ); + } + return( true ); +} + +#endif // __fpa__Base__FastMarching__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Functors/GaussianModel.h b/libs/fpa/Base/Functors/GaussianModel.h new file mode 100644 index 0000000..8d32f51 --- /dev/null +++ b/libs/fpa/Base/Functors/GaussianModel.h @@ -0,0 +1,82 @@ +#ifndef __fpa__Base__Functors__GaussianModel__h__ +#define __fpa__Base__Functors__GaussianModel__h__ + +#include +#include +#include + +namespace fpa +{ + namespace Base + { + namespace Functors + { + /** + */ + template< class _TInput, class _TOutput > + class GaussianModel + : public itk::FunctionBase< _TInput, _TOutput > + { + public: + typedef GaussianModel Self; + typedef itk::FunctionBase< _TInput, _TOutput > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TInput TInput; + typedef _TOutput TOutput; + + // Model estimator + typedef + cpExtensions::Algorithms::IterativeGaussianModelEstimator< TOutput, 1 > + TModel; + + public: + itkNewMacro( Self ); + itkTypeMacro( GaussianModel, itk::FunctionBase ); + + /* TODO + itkGetConstMacro( SupportSize, unsigned int ); + itkGetConstMacro( MinimumCost, TOutput ); + itkGetObjectMacro( Model, TModel ); + itkGetConstObjectMacro( Model, TModel ); + itkSetMacro( SupportSize, unsigned int ); + itkSetMacro( MinimumCost, TOutput ); + */ + + public: + virtual TOutput Evaluate( const TInput& x ) const override; + + protected: + GaussianModel( ); + virtual ~GaussianModel( ); + + private: + // Purposely not implemented + GaussianModel( const Self& other ); + Self& operator=( const Self& other ); + + protected: + /* TODO + unsigned int m_SupportSize; + TOutput m_MinimumCost; + typename TModel::Pointer m_Model; + */ + mutable double m_S1; + mutable double m_S2; + mutable unsigned long m_N; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__Functors__GaussianModel__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Functors/GaussianModel.hxx b/libs/fpa/Base/Functors/GaussianModel.hxx new file mode 100644 index 0000000..6a44ff8 --- /dev/null +++ b/libs/fpa/Base/Functors/GaussianModel.hxx @@ -0,0 +1,100 @@ +#ifndef __fpa__Base__Functors__GaussianModel__hxx__ +#define __fpa__Base__Functors__GaussianModel__hxx__ + +// ------------------------------------------------------------------------- +template< class _TInput, class _TOutput > +typename fpa::Base::Functors::GaussianModel< _TInput, _TOutput >:: +TOutput fpa::Base::Functors::GaussianModel< _TInput, _TOutput >:: +Evaluate( const TInput& x ) const +{ + double v = double( x ); + double d = double( 0 ); + if( this->m_N > 1 ) + { + double N = double( this->m_N ); + double m = this->m_S1 / N; + double s = + ( this->m_S2 - ( ( this->m_S1 * this->m_S1 ) / N ) ) / + ( N - double( 1 ) ); + if( s > double( 0 ) ) + { + d = ( v - m ); + d *= d; + d /= s; + d = std::sqrt( d ); + if( d <= double( 1.5 ) ) // 2sigma + { + this->m_S1 += v; + this->m_S2 += v * v; + this->m_N += 1; + + } // fi + } + else + { + this->m_S1 += v; + this->m_S2 += v * v; + this->m_N += 1; + + } // fi + } + else + { + this->m_S1 += v; + this->m_S2 += v * v; + this->m_N += 1; + + } // fi + return( d ); + + /* TODO + TOutput r = TOutput( 0 ); + if( this->m_Model->GetNumberOfSamples( ) > this->m_SupportSize ) + { + r = this->m_Model->SquaredMahalanobis( TOutput( x ) ); + if( r <= TOutput( 1 ) ) + this->m_Model->AddSample( TOutput( x ) ); + } + else + { + this->m_Model->AddSample( TOutput( x ) ); + if( this->m_Model->GetNumberOfSamples( ) > 2 ) + r = this->m_Model->SquaredMahalanobis( TOutput( x ) ); + + } // fi + if( r < this->m_MinimumCost ) + return( this->m_MinimumCost ); + else + return( r ); + */ +} + +// ------------------------------------------------------------------------- +template< class _TInput, class _TOutput > +fpa::Base::Functors::GaussianModel< _TInput, _TOutput >:: +GaussianModel( ) + : Superclass( ) + /* TODO + , + m_SupportSize( 30 ), + m_MinimumCost( TOutput( 1e-10 ) ) + */ +{ + /* TODO + this->m_Model = TModel::New( ); + */ + this->m_S1 = double( 0 ); + this->m_S2 = double( 0 ); + this->m_N = 0; +} + +// ------------------------------------------------------------------------- +template< class _TInput, class _TOutput > +fpa::Base::Functors::GaussianModel< _TInput, _TOutput >:: +~GaussianModel( ) +{ +} + +#endif // __fpa__Base__Functors__GaussianModel__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Functors/Inverse.h b/libs/fpa/Base/Functors/Inverse.h new file mode 100644 index 0000000..af68022 --- /dev/null +++ b/libs/fpa/Base/Functors/Inverse.h @@ -0,0 +1,63 @@ +#ifndef __fpa__Base__Functors__Inverse__h__ +#define __fpa__Base__Functors__Inverse__h__ + +#include +#include + +namespace fpa +{ + namespace Base + { + namespace Functors + { + /** + */ + template< class _TInput, class _TOutput > + class Inverse + : public itk::FunctionBase< _TInput, _TOutput > + { + public: + typedef Inverse Self; + typedef itk::FunctionBase< _TInput, _TOutput > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TInput TInput; + typedef _TOutput TOutput; + + public: + itkNewMacro( Self ); + itkTypeMacro( Inverse, itk::FunctionBase ); + + itkGetConstMacro( NegativeValue, _TOutput ); + itkSetMacro( NegativeValue, _TOutput ); + + public: + virtual TOutput Evaluate( const TInput& x ) const override; + + protected: + Inverse( ); + virtual ~Inverse( ); + + private: + // Purposely not implemented + Inverse( const Self& other ); + Self& operator=( const Self& other ); + + protected: + _TOutput m_NegativeValue; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__Functors__Inverse__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Functors/Inverse.hxx b/libs/fpa/Base/Functors/Inverse.hxx new file mode 100644 index 0000000..4574cae --- /dev/null +++ b/libs/fpa/Base/Functors/Inverse.hxx @@ -0,0 +1,36 @@ +#ifndef __fpa__Base__Functors__Inverse__hxx__ +#define __fpa__Base__Functors__Inverse__hxx__ + +// ------------------------------------------------------------------------- +template< class _TInput, class _TOutput > +typename fpa::Base::Functors::Inverse< _TInput, _TOutput >:: +TOutput fpa::Base::Functors::Inverse< _TInput, _TOutput >:: +Evaluate( const TInput& x ) const +{ + TInput sign = TInput( ( x < TInput( 0 ) )? -1: 1 ); + TOutput y = TOutput( 1 ) / ( TOutput( 1 ) + TOutput( x * sign ) ); + if( sign < TInput( 0 ) ) + return( y * this->m_NegativeValue ); + else + return( y ); +} + +// ------------------------------------------------------------------------- +template< class _TInput, class _TOutput > +fpa::Base::Functors::Inverse< _TInput, _TOutput >:: +Inverse( ) + : Superclass( ), + m_NegativeValue( _TOutput( -1 ) ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInput, class _TOutput > +fpa::Base::Functors::Inverse< _TInput, _TOutput >:: +~Inverse( ) +{ +} + +#endif // __fpa__Base__Functors__Inverse__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Functors/RegionGrow/Base.h b/libs/fpa/Base/Functors/RegionGrow/Base.h new file mode 100644 index 0000000..d2c37aa --- /dev/null +++ b/libs/fpa/Base/Functors/RegionGrow/Base.h @@ -0,0 +1,72 @@ +#ifndef __fpa__Base__Functors__RegionGrow__Base__h__ +#define __fpa__Base__Functors__RegionGrow__Base__h__ + +#include +#include +#include + +namespace fpa +{ + namespace Base + { + namespace Functors + { + namespace RegionGrow + { + /** + */ + template< class _TVertex, class _TOutput > + class Base + : public fpa::Base::Functors::VertexCostFunctionBase< _TVertex, _TOutput > + { + public: + typedef fpa::Base::Functors::VertexCostFunctionBase< _TVertex, _TOutput > Superclass; + typedef Base Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TVertex TVertex; + typedef _TOutput TOutput; + + public: + itkTypeMacro( Base, itk::Object ); + + itkGetConstMacro( InsideValue, _TOutput ); + itkGetConstMacro( OutsideValue, _TOutput ); + + itkSetMacro( InsideValue, _TOutput ); + itkSetMacro( OutsideValue, _TOutput ); + + public: + virtual TOutput Evaluate( const TVertex& a, const TVertex& b ) const = 0; + + protected: + Base( ) + : Superclass( ), + m_InsideValue( TOutput( 1 ) ), + m_OutsideValue( TOutput( 0 ) ) + { } + virtual ~Base( ) + { } + + private: + // Purposely not defined + Base( const Self& other ); + Self& operator=( const Self& other ); + + protected: + _TOutput m_InsideValue; + _TOutput m_OutsideValue; + }; + + } // ecapseman + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Base__Functors__RegionGrow__Base__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Functors/RegionGrow/Tautology.h b/libs/fpa/Base/Functors/RegionGrow/Tautology.h new file mode 100644 index 0000000..11d016e --- /dev/null +++ b/libs/fpa/Base/Functors/RegionGrow/Tautology.h @@ -0,0 +1,62 @@ +#ifndef __fpa__Base__Functors__RegionGrow__Tautology__h__ +#define __fpa__Base__Functors__RegionGrow__Tautology__h__ + +#include + +namespace fpa +{ + namespace Base + { + namespace Functors + { + namespace RegionGrow + { + /** + */ + template< class _TVertex, class _TOutput > + class Tautology + : public Base< _TVertex, _TOutput > + { + public: + typedef Tautology Self; + typedef Base< _TVertex, _TOutput > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TVertex TVertex; + + public: + itkNewMacro( Self ); + itkTypeMacro( Tautology, Base ); + + public: + virtual TOutput Evaluate( const TVertex& a, const TVertex& b ) const override + { + return( this->m_InsideValue ); + } + + protected: + Tautology( ) + : Superclass( ) + { } + virtual ~Tautology( ) + { } + + private: + // Purposely not defined + Tautology( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Base__Functors__RegionGrow__Tautology__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/Functors/VertexCostFunctionBase.h b/libs/fpa/Base/Functors/VertexCostFunctionBase.h new file mode 100644 index 0000000..12f5c24 --- /dev/null +++ b/libs/fpa/Base/Functors/VertexCostFunctionBase.h @@ -0,0 +1,55 @@ +#ifndef __fpa__Base__Functors__VertexCostFunctionBase__h__ +#define __fpa__Base__Functors__VertexCostFunctionBase__h__ + +#include +#include + +namespace fpa +{ + namespace Base + { + namespace Functors + { + /** + */ + template< class _TVertex, class _TOutput > + class VertexCostFunctionBase + : public itk::Object + { + public: + typedef VertexCostFunctionBase Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TVertex TVertex; + typedef _TOutput TOutput; + + public: + itkTypeMacro( VertexCostFunctionBase, itk::Object ); + + public: + virtual TOutput Evaluate( const TVertex& a, const TVertex& b ) const = 0; + + protected: + VertexCostFunctionBase( ) + : Superclass( ) + { } + virtual ~VertexCostFunctionBase( ) + { } + + private: + // Purposely not defined + VertexCostFunctionBase( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Base__Functors__VertexCostFunctionBase__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/MinimumSpanningTree.h b/libs/fpa/Base/MinimumSpanningTree.h new file mode 100644 index 0000000..03940f1 --- /dev/null +++ b/libs/fpa/Base/MinimumSpanningTree.h @@ -0,0 +1,80 @@ +#ifndef __fpa__Base__MinimumSpanningTree__h__ +#define __fpa__Base__MinimumSpanningTree__h__ + +#include +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TVertex, class _TPath, class _TSuperclass > + class MinimumSpanningTree + : public _TSuperclass + { + public: + typedef MinimumSpanningTree Self; + typedef _TSuperclass Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TVertex TVertex; + typedef _TPath TPath; + typedef std::pair< _TVertex, bool > TCollision; + typedef std::vector< TCollision > TCollisionsRow; + typedef std::vector< TCollisionsRow > TCollisions; + + protected: + typedef std::vector< unsigned long > _TRow; + typedef std::vector< _TRow > _TMatrix; + + public: + itkTypeMacro( MinimumSpanningTree, _TSuperclass ); + + public: + const TCollisions& GetCollisions( ) const; + void SetCollisions( const TCollisions& collisions ); + + void ClearSeeds( ); + void AddSeed( const _TVertex& seed ); + + virtual _TVertex GetParent( const _TVertex& v ) const = 0; + virtual void SetParent( const _TVertex& v, const _TVertex& p ) = 0; + + virtual void GetPath( + typename _TPath::Pointer& path, const _TVertex& a + ) const; + virtual void GetPath( + typename _TPath::Pointer& path, const _TVertex& a, const _TVertex& b + ) const; + + protected: + MinimumSpanningTree( ); + virtual ~MinimumSpanningTree( ); + + private: + // Purposely not defined + MinimumSpanningTree( const Self& other ); + Self& operator=( const Self& other ); + + protected: + TCollisions m_Collisions; + _TMatrix m_FrontPaths; + std::vector< _TVertex > m_Seeds; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__MinimumSpanningTree__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/MinimumSpanningTree.hxx b/libs/fpa/Base/MinimumSpanningTree.hxx new file mode 100644 index 0000000..e8ed6a2 --- /dev/null +++ b/libs/fpa/Base/MinimumSpanningTree.hxx @@ -0,0 +1,244 @@ +#ifndef __fpa__Base__MinimumSpanningTree__hxx__ +#define __fpa__Base__MinimumSpanningTree__hxx__ + +#include + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TPath, class _TSuperclass > +const typename +fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +TCollisions& +fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +GetCollisions( ) const +{ + return( this->m_Collisions ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TPath, class _TSuperclass > +void fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +SetCollisions( const TCollisions& collisions ) +{ + static const unsigned long _inf = + std::numeric_limits< unsigned long >::max( ); + if( this->m_Collisions == collisions ) + return; + + this->m_Collisions = collisions; + + // Prepare a front graph + unsigned long N = this->m_Collisions.size( ); + _TMatrix dist( N, _TRow( N, _inf ) ); + this->m_FrontPaths = dist; + for( unsigned long i = 0; i < N; ++i ) + { + for( unsigned long j = 0; j < N; ++j ) + { + if( this->m_Collisions[ i ][ j ].second ) + { + dist[ i ][ j ] = 1; + dist[ j ][ i ] = 1; + this->m_FrontPaths[ i ][ j ] = j; + this->m_FrontPaths[ j ][ i ] = i; + + } // fi + + } // rof + dist[ i ][ i ] = 0; + this->m_FrontPaths[ i ][ i ] = i; + + } // rof + + // Use Floyd-Warshall to compute all possible paths between fronts + for( unsigned long k = 0; k < N; ++k ) + { + for( unsigned long i = 0; i < N; ++i ) + { + for( unsigned long j = 0; j < N; ++j ) + { + // WARNING: you don't want a numeric overflow!!! + unsigned long dik = dist[ i ][ k ]; + unsigned long dkj = dist[ k ][ j ]; + unsigned long sum = _inf; + if( dik < _inf && dkj < _inf ) + sum = dik + dkj; + + // Ok, continue Floyd-Warshall + if( sum < dist[ i ][ j ] ) + { + dist[ i ][ j ] = sum; + this->m_FrontPaths[ i ][ j ] = this->m_FrontPaths[ i ][ k ]; + + } // fi + + } // rof + + } // rof + + } // rof + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TPath, class _TSuperclass > +void fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +ClearSeeds( ) +{ + this->m_Seeds.clear( ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TPath, class _TSuperclass > +void fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +AddSeed( const _TVertex& seed ) +{ + this->m_Seeds.push_back( seed ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TPath, class _TSuperclass > +void fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +GetPath( typename _TPath::Pointer& path, const _TVertex& a ) const +{ + std::vector< _TVertex > vertices; + _TVertex it = a; + _TVertex p = this->GetParent( it ); + while( it != p ) + { + vertices.push_back( it ); + it = p; + p = this->GetParent( it ); + + } // elihw + vertices.push_back( it ); + + if( path.IsNull( ) ) + path = _TPath::New( ); + for( auto v = vertices.begin( ); v != vertices.end( ); ++v ) + path->AddVertex( *v ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TPath, class _TSuperclass > +void fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +GetPath( + typename _TPath::Pointer& path, const _TVertex& a, const _TVertex& b + ) const +{ + static const unsigned long _inf = + std::numeric_limits< unsigned long >::max( ); + + if( path.IsNull( ) ) + path = _TPath::New( ); + typename _TPath::Pointer pa, pb; + this->GetPath( pa, a ); + this->GetPath( pb, b ); + if( pa->GetSize( ) > 0 && pb->GetSize( ) > 0 ) + { + // Find front identifiers + unsigned long ia = _inf, ib = _inf; + unsigned long N = this->m_Seeds.size( ); + for( unsigned long i = 0; i < N; ++i ) + { + if( this->m_Seeds[ i ] == pa->GetVertex( pa->GetSize( ) - 1 ) ) + ia = i; + if( this->m_Seeds[ i ] == pb->GetVertex( pb->GetSize( ) - 1 ) ) + ib = i; + + } // rof + + // Check if there is a front-jump between given seeds + if( ia != ib ) + { + // Compute front path + std::vector< long > fpath; + fpath.push_back( ia ); + while( ia != ib ) + { + ia = this->m_FrontPaths[ ia ][ ib ]; + fpath.push_back( ia ); + + } // elihw + + // Continue only if both fronts are connected + unsigned int N = fpath.size( ); + if( N > 0 ) + { + // First path: from start vertex to first collision + this->GetPath( + path, a, + this->m_Collisions[ fpath[ 0 ] ][ fpath[ 1 ] ].first + ); + + // Intermediary paths + for( unsigned int i = 1; i < N - 1; ++i ) + { + typename _TPath::Pointer ipath; + this->GetPath( + ipath, + this->m_Collisions[ fpath[ i ] ][ fpath[ i - 1 ] ].first, + this->m_Collisions[ fpath[ i ] ][ fpath[ i + 1 ] ].first + ); + for( long id = 0; id < ipath->GetSize( ); ++id ) + path->AddVertex( ipath->GetVertex( id ) ); + + } // rof + + // Final path: from last collision to end point + typename _TPath::Pointer lpath; + this->GetPath( + lpath, + this->m_Collisions[ fpath[ N - 1 ] ][ fpath[ N - 2 ] ].first, b + ); + for( long id = 0; id < lpath->GetSize( ); ++id ) + path->AddVertex( lpath->GetVertex( id ) ); + + } // fi + } + else + { + // Ignore common part: find common ancestor + long aIt = pa->GetSize( ) - 1; + long bIt = pb->GetSize( ) - 1; + bool cont = true; + while( aIt >= 0 && bIt >= 0 && cont ) + { + cont = ( pa->GetVertex( aIt ) == pb->GetVertex( bIt ) ); + aIt--; + bIt--; + + } // elihw + aIt++; + bIt++; + + // Glue both parts + for( long cIt = 0; cIt <= aIt; ++cIt ) + path->AddVertex( pa->GetVertex( cIt ) ); + for( ; bIt >= 0; --bIt ) + path->AddVertex( pb->GetVertex( bIt ) ); + + } // fi + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TPath, class _TSuperclass > +fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +MinimumSpanningTree( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TPath, class _TSuperclass > +fpa::Base::MinimumSpanningTree< _TVertex, _TPath, _TSuperclass >:: +~MinimumSpanningTree( ) +{ +} + +#endif // __fpa__Base__MinimumSpanningTree__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/PriorityQueueAlgorithm.h b/libs/fpa/Base/PriorityQueueAlgorithm.h new file mode 100644 index 0000000..e34b823 --- /dev/null +++ b/libs/fpa/Base/PriorityQueueAlgorithm.h @@ -0,0 +1,64 @@ +#ifndef __fpa__Base__PriorityQueueAlgorithm__h__ +#define __fpa__Base__PriorityQueueAlgorithm__h__ + +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TSuperclass > + class PriorityQueueAlgorithm + : public _TSuperclass + { + public: + typedef PriorityQueueAlgorithm Self; + typedef _TSuperclass Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + protected: + typedef typename Superclass::_TQueueNode _TQueueNode; + struct _TQueueNodeCompare + { + bool operator( )( const _TQueueNode& a, const _TQueueNode& b ) + { + return( b.Result < a.Result ); + } + }; + typedef std::vector< _TQueueNode > _TQueue; + + public: + itkTypeMacro( PriorityQueueAlgorithm, Algorithm ); + + protected: + PriorityQueueAlgorithm( ); + virtual ~PriorityQueueAlgorithm( ); + + virtual unsigned long _QueueSize( ) const override; + virtual void _QueueClear( ) override; + virtual void _QueuePush( const _TQueueNode& node ) override; + virtual _TQueueNode _QueuePop( ) override; + + private: + // Purposely not defined + PriorityQueueAlgorithm( const Self& other ); + Self& operator=( const Self& other ); + + protected: + _TQueue m_Queue; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__PriorityQueueAlgorithm__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/PriorityQueueAlgorithm.hxx b/libs/fpa/Base/PriorityQueueAlgorithm.hxx new file mode 100644 index 0000000..11c1408 --- /dev/null +++ b/libs/fpa/Base/PriorityQueueAlgorithm.hxx @@ -0,0 +1,60 @@ +#ifndef __fpa__Base__PriorityQueueAlgorithm__hxx__ +#define __fpa__Base__PriorityQueueAlgorithm__hxx__ + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +fpa::Base::PriorityQueueAlgorithm< _TSuperclass >:: +PriorityQueueAlgorithm( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +fpa::Base::PriorityQueueAlgorithm< _TSuperclass >:: +~PriorityQueueAlgorithm( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +unsigned long fpa::Base::PriorityQueueAlgorithm< _TSuperclass >:: +_QueueSize( ) const +{ + return( this->m_Queue.size( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +void fpa::Base::PriorityQueueAlgorithm< _TSuperclass >:: +_QueueClear( ) +{ + this->m_Queue.clear( ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +void fpa::Base::PriorityQueueAlgorithm< _TSuperclass >:: +_QueuePush( const _TQueueNode& node ) +{ + static _TQueueNodeCompare cmp; + this->m_Queue.push_back( node ); + std::push_heap( this->m_Queue.begin( ), this->m_Queue.end( ), cmp ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +typename fpa::Base::PriorityQueueAlgorithm< _TSuperclass >:: +_TQueueNode fpa::Base::PriorityQueueAlgorithm< _TSuperclass >:: +_QueuePop( ) +{ + static _TQueueNodeCompare cmp; + std::pop_heap( this->m_Queue.begin( ), this->m_Queue.end( ), cmp ); + _TQueueNode f = this->m_Queue.back( ); + this->m_Queue.pop_back( ); + return( f ); +} + +#endif // __fpa__Base__PriorityQueueAlgorithm__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/QueueAlgorithm.h b/libs/fpa/Base/QueueAlgorithm.h new file mode 100644 index 0000000..dd077f7 --- /dev/null +++ b/libs/fpa/Base/QueueAlgorithm.h @@ -0,0 +1,59 @@ +#ifndef __fpa__Base__QueueAlgorithm__h__ +#define __fpa__Base__QueueAlgorithm__h__ + +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TSuperclass > + class QueueAlgorithm + : public _TSuperclass + { + public: + typedef QueueAlgorithm Self; + typedef _TSuperclass Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + protected: + typedef typename Superclass::_TQueueNode _TQueueNode; + typedef std::queue< _TQueueNode > _TQueue; + + public: + itkTypeMacro( QueueAlgorithm, _TSuperclass ); + + protected: + QueueAlgorithm( ); + virtual ~QueueAlgorithm( ); + + virtual unsigned long _QueueSize( ) const override; + virtual void _QueueClear( ) override; + virtual void _QueuePush( const _TQueueNode& node ) override; + virtual _TQueueNode _QueuePop( ) override; + + private: + // Purposely not defined + QueueAlgorithm( const Self& other ); + Self& operator=( const Self& other ); + + protected: + _TQueue m_Queue; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__QueueAlgorithm__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/QueueAlgorithm.hxx b/libs/fpa/Base/QueueAlgorithm.hxx new file mode 100644 index 0000000..db5b8bc --- /dev/null +++ b/libs/fpa/Base/QueueAlgorithm.hxx @@ -0,0 +1,57 @@ +#ifndef __fpa__Base__QueueAlgorithm__hxx__ +#define __fpa__Base__QueueAlgorithm__hxx__ + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +fpa::Base::QueueAlgorithm< _TSuperclass >:: +QueueAlgorithm( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +fpa::Base::QueueAlgorithm< _TSuperclass >:: +~QueueAlgorithm( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +unsigned long fpa::Base::QueueAlgorithm< _TSuperclass >:: +_QueueSize( ) const +{ + return( this->m_Queue.size( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +void fpa::Base::QueueAlgorithm< _TSuperclass >:: +_QueueClear( ) +{ + while( this->m_Queue.size( ) > 0 ) + this->m_Queue.pop( ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +void fpa::Base::QueueAlgorithm< _TSuperclass >:: +_QueuePush( const _TQueueNode& node ) +{ + this->m_Queue.push( node ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +typename fpa::Base::QueueAlgorithm< _TSuperclass >:: +_TQueueNode fpa::Base::QueueAlgorithm< _TSuperclass >:: +_QueuePop( ) +{ + _TQueueNode f = this->m_Queue.front( ); + this->m_Queue.pop( ); + return( f ); +} + +#endif // __fpa__Base__QueueAlgorithm__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/RegionGrow.h b/libs/fpa/Base/RegionGrow.h new file mode 100644 index 0000000..4790b85 --- /dev/null +++ b/libs/fpa/Base/RegionGrow.h @@ -0,0 +1,79 @@ +#ifndef __fpa__Base__RegionGrow__h__ +#define __fpa__Base__RegionGrow__h__ + +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TSuperclass > + class RegionGrow + : public fpa::Base::QueueAlgorithm< _TSuperclass > + { + public: + typedef RegionGrow Self; + typedef fpa::Base::QueueAlgorithm< _TSuperclass > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TVertex TVertex; + typedef fpa::Base::Functors::RegionGrow::Base< TVertex, TOutput > TGrowFunction; + + protected: + typedef typename Superclass::_TQueueNode _TQueueNode; + + public: + itkTypeMacro( RegionGrow, Algorithm ); + + public: + TGrowFunction* GetGrowFunction( ); + const TGrowFunction* GetGrowFunction( ) const; + TOutput GetInsideValue( ) const; + TOutput GetOutsideValue( ) const; + + void SetGrowFunction( TGrowFunction* f ); + void SetInsideValue( const TOutput& v ); + void SetOutsideValue( const TOutput& v ); + + protected: + RegionGrow( ); + virtual ~RegionGrow( ); + + virtual bool _UpdateValue( + _TQueueNode& v, const _TQueueNode& p + ) override; + virtual TOutput _GetInputValue( const _TQueueNode& v, const _TQueueNode& p ) override + { + TOutput res = this->m_InitResult; + if( this->m_GrowFunction.IsNotNull( ) ) + res = this->m_GrowFunction->Evaluate( v.Vertex, p.Vertex ); + return( res ); + } + + + private: + // Purposely not defined + RegionGrow( const Self& other ); + Self& operator=( const Self& other ); + + protected: + typename TGrowFunction::Pointer m_GrowFunction; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__RegionGrow__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Base/RegionGrow.hxx b/libs/fpa/Base/RegionGrow.hxx new file mode 100644 index 0000000..c3cb9f3 --- /dev/null +++ b/libs/fpa/Base/RegionGrow.hxx @@ -0,0 +1,128 @@ +#ifndef __fpa__Base__RegionGrow__hxx__ +#define __fpa__Base__RegionGrow__hxx__ + +#include + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +typename fpa::Base::RegionGrow< _TSuperclass >:: +TGrowFunction* fpa::Base::RegionGrow< _TSuperclass >:: +GetGrowFunction( ) +{ + // TODO: return( dynamic_cast< TGrowFunction* >( this->GetVertexFunction( ) ) ); + return( this->m_GrowFunction ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +const typename fpa::Base::RegionGrow< _TSuperclass >:: +TGrowFunction* fpa::Base::RegionGrow< _TSuperclass >:: +GetGrowFunction( ) const +{ + /* TODO + return( + dynamic_cast< const TGrowFunction* >( this->GetVertexFunction( ) ) + ); + */ + return( this->m_GrowFunction ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +typename fpa::Base::RegionGrow< _TSuperclass >:: +TOutput fpa::Base::RegionGrow< _TSuperclass >:: +GetInsideValue( ) const +{ + const TGrowFunction* f = this->GetGrowFunction( ); + if( f != NULL ) + return( f->GetInsideValue( ) ); + else + return( this->m_InitResult ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +typename fpa::Base::RegionGrow< _TSuperclass >:: +TOutput fpa::Base::RegionGrow< _TSuperclass >:: +GetOutsideValue( ) const +{ + const TGrowFunction* f = this->GetGrowFunction( ); + if( f != NULL ) + return( f->GetOutsideValue( ) ); + else + return( this->m_InitResult ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +void fpa::Base::RegionGrow< _TSuperclass >:: +SetGrowFunction( TGrowFunction* f ) +{ + TGrowFunction* old_f = this->GetGrowFunction( ); + if( old_f != NULL ) + { + f->SetInsideValue( old_f->GetInsideValue( ) ); + f->SetOutsideValue( old_f->GetOutsideValue( ) ); + + } // fi + this->m_GrowFunction = f; +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +void fpa::Base::RegionGrow< _TSuperclass >:: +SetInsideValue( const TOutput& v ) +{ + TGrowFunction* f = this->GetGrowFunction( ); + if( f != NULL ) + { + f->SetInsideValue( v ); + this->Modified( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +void fpa::Base::RegionGrow< _TSuperclass >:: +SetOutsideValue( const TOutput& v ) +{ + TGrowFunction* f = this->GetGrowFunction( ); + if( f != NULL ) + { + f->SetOutsideValue( v ); + this->Modified( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +fpa::Base::RegionGrow< _TSuperclass >:: +RegionGrow( ) + : Superclass( ) +{ + typedef fpa::Base::Functors::RegionGrow::Tautology< TVertex, TOutput > _TFunc; + this->SetGrowFunction( _TFunc::New( ) ); + this->m_InitResult = this->GetGrowFunction( )->GetOutsideValue( ); +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +fpa::Base::RegionGrow< _TSuperclass >:: +~RegionGrow( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSuperclass > +bool fpa::Base::RegionGrow< _TSuperclass >:: +_UpdateValue( _TQueueNode& v, const _TQueueNode& p ) +{ + v.Result = this->_GetInputValue( v, p ); + return( v.Result == this->GetInsideValue( ) ); +} + +#endif // __fpa__Base__RegionGrow__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Config.h b/libs/fpa/Config.h new file mode 100644 index 0000000..9728f0b --- /dev/null +++ b/libs/fpa/Config.h @@ -0,0 +1,12 @@ +#ifndef __fpa__Config__h__ +#define __fpa__Config__h__ + +/* + * ========================================================================= + * Language related macros + * ========================================================================= + */ + +#endif // __fpa__Config__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Algorithm.h b/libs/fpa/Image/Algorithm.h new file mode 100644 index 0000000..0105cc3 --- /dev/null +++ b/libs/fpa/Image/Algorithm.h @@ -0,0 +1,74 @@ +#ifndef __fpa__Image__Algorithm__h__ +#define __fpa__Image__Algorithm__h__ + +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TInputImage, class _TOutputImage > + class Algorithm + : public fpa::Base::Algorithm< itk::ImageToImageFilter< _TInputImage, _TOutputImage >, typename _TInputImage::IndexType, typename _TOutputImage::PixelType > + { + public: + typedef itk::ImageToImageFilter< _TInputImage, _TOutputImage > TFilter; + typedef typename _TInputImage::IndexType TVertex; + typedef typename _TOutputImage::PixelType TOutput; + + typedef Algorithm Self; + typedef fpa::Base::Algorithm< TFilter, TVertex, TOutput > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TFrontId TFrontId; + typedef typename Superclass::TNeighborhood TNeighborhood; + + typedef fpa::Image::Functors::Base< itk::ImageBase< _TInputImage::ImageDimension >, typename Superclass::TNeighborhoodFunction > TNeighborhoodFunction; + typedef fpa::Image::Functors::Base< _TInputImage, typename Superclass::TVertexFunction > TVertexFunction; + + protected: + typedef typename Superclass::_TQueueNode _TQueueNode; + + public: + itkTypeMacro( fpa::Image::Algorithm, fpa::Base::Algorithm ); + + protected: + Algorithm( ); + virtual ~Algorithm( ); + + virtual void _BeforeGenerateData( ) override; + virtual void _InitMarks( ) override; + virtual void _InitResults( const TOutput& init_value ) override; + virtual bool _IsMarked( const TVertex& v ) const override; + virtual void _Mark( const _TQueueNode& n ) override; + virtual TFrontId _GetMark( const TVertex& v ) const override; + virtual void _UpdateResult( const _TQueueNode& n ) override; + virtual TOutput _GetResult( const TVertex& v ) const override; + virtual unsigned int _GetNumberOfDimensions( ) const override; + + private: + // Purposely not defined + Algorithm( const Self& other ); + Self& operator=( const Self& other ); + + protected: + unsigned int m_MarksIdx; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__Algorithm__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Algorithm.hxx b/libs/fpa/Image/Algorithm.hxx new file mode 100644 index 0000000..3c33e4e --- /dev/null +++ b/libs/fpa/Image/Algorithm.hxx @@ -0,0 +1,149 @@ +#ifndef __fpa__Image__Algorithm__hxx__ +#define __fpa__Image__Algorithm__hxx__ + +// Send Piotr's code to Anna + +#include +#include + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +Algorithm( ) + : Superclass( ) +{ + typedef fpa::Image::Functors::SimpleNeighborhood< _TInputImage::ImageDimension > _TNeigh; + typedef itk::Image< TFrontId, _TInputImage::ImageDimension > _TMarks; + + this->m_MarksIdx = this->GetNumberOfRequiredOutputs( ); + this->itk::ProcessObject::SetNumberOfRequiredOutputs( this->m_MarksIdx + 1 ); + typename _TMarks::Pointer marks = _TMarks::New( ); + this->SetNthOutput( this->m_MarksIdx, marks ); + + typename _TNeigh::Pointer neigh = _TNeigh::New( ); + this->SetNeighborhoodFunction( neigh ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +~Algorithm( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_BeforeGenerateData( ) +{ + this->Superclass::_BeforeGenerateData( ); + this->AllocateOutputs( ); + + TNeighborhoodFunction* neighFunc = + dynamic_cast< TNeighborhoodFunction* >( + this->GetNeighborhoodFunction( ) + ); + if( neighFunc == NULL ) + itkExceptionMacro( << "NeighborhoodFunction not well defined." ); + neighFunc->SetImage( this->GetInput( ) ); + + TVertexFunction* vertexFunc = + dynamic_cast< TVertexFunction* >( + this->GetVertexFunction( ) + ); + if( vertexFunc != NULL ) + vertexFunc->SetImage( this->GetInput( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_InitMarks( ) +{ + typedef itk::Image< TFrontId, _TInputImage::ImageDimension > _TMarks; + _TMarks* marks = + dynamic_cast< _TMarks* >( + this->itk::ProcessObject::GetOutput( this->m_MarksIdx ) + ); + marks->FillBuffer( 0 ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_InitResults( const TOutput& init_value ) +{ + this->GetOutput( )->FillBuffer( init_value ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +bool fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_IsMarked( const TVertex& v ) const +{ + typedef itk::Image< TFrontId, _TInputImage::ImageDimension > _TMarks; + const _TMarks* marks = + dynamic_cast< const _TMarks* >( + this->itk::ProcessObject::GetOutput( this->m_MarksIdx ) + ); + return( marks->GetPixel( v ) != 0 ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_Mark( const _TQueueNode& n ) +{ + typedef itk::Image< TFrontId, _TInputImage::ImageDimension > _TMarks; + _TMarks* marks = + dynamic_cast< _TMarks* >( + this->itk::ProcessObject::GetOutput( this->m_MarksIdx ) + ); + marks->SetPixel( n.Vertex, n.FrontId ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +typename fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +TFrontId fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_GetMark( const TVertex& v ) const +{ + typedef itk::Image< TFrontId, _TInputImage::ImageDimension > _TMarks; + const _TMarks* marks = + dynamic_cast< const _TMarks* >( + this->itk::ProcessObject::GetOutput( this->m_MarksIdx ) + ); + return( marks->GetPixel( v ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_UpdateResult( const _TQueueNode& n ) +{ + this->GetOutput( )->SetPixel( n.Vertex, n.Result ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +typename fpa::Image::Algorithm< _TInputImage, _TOutputImage >::TOutput +fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_GetResult( const TVertex& v ) const +{ + if( this->GetOutput( )->GetLargestPossibleRegion( ).IsInside( v ) ) + return( this->GetOutput( )->GetPixel( v ) ); + else + return( this->m_InitResult ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +unsigned int fpa::Image::Algorithm< _TInputImage, _TOutputImage >:: +_GetNumberOfDimensions( ) const +{ + return( _TInputImage::ImageDimension ); +} + +#endif // __fpa__Image__Algorithm__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Dijkstra.h b/libs/fpa/Image/Dijkstra.h new file mode 100644 index 0000000..118cc27 --- /dev/null +++ b/libs/fpa/Image/Dijkstra.h @@ -0,0 +1,54 @@ +#ifndef __fpa__Image__Dijkstra__h__ +#define __fpa__Image__Dijkstra__h__ + +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TInputImage, class _TOutputImage > + class Dijkstra + : public fpa::Base::Dijkstra< fpa::Image::Algorithm< _TInputImage, _TOutputImage >, fpa::Image::MinimumSpanningTree< _TInputImage::ImageDimension > > + { + public: + typedef Dijkstra Self; + typedef fpa::Image::Algorithm< _TInputImage, _TOutputImage > TAlgorithm; + typedef fpa::Image::MinimumSpanningTree< _TInputImage::ImageDimension > TMST; + typedef fpa::Base::Dijkstra< TAlgorithm, TMST > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TVertex TVertex; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::Dijkstra, fpa::Base::Dijkstra ); + + protected: + Dijkstra( ); + virtual ~Dijkstra( ); + + private: + // Purposely not defined + Dijkstra( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__Dijkstra__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Dijkstra.hxx b/libs/fpa/Image/Dijkstra.hxx new file mode 100644 index 0000000..1f21dc4 --- /dev/null +++ b/libs/fpa/Image/Dijkstra.hxx @@ -0,0 +1,26 @@ +#ifndef __fpa__Image__Dijkstra__hxx__ +#define __fpa__Image__Dijkstra__hxx__ + +#include + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::Dijkstra< _TInputImage, _TOutputImage >:: +Dijkstra( ) + : Superclass( ) +{ + typedef fpa::Image::Functors::VertexCost< _TInputImage, TOutput > _TCost; + typename _TCost::Pointer cost = _TCost::New( ); + this->SetVertexFunction( cost ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::Dijkstra< _TInputImage, _TOutputImage >:: +~Dijkstra( ) +{ +} + +#endif // __fpa__Image__Dijkstra__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/EndPointsFilter.h b/libs/fpa/Image/EndPointsFilter.h new file mode 100644 index 0000000..afbfb10 --- /dev/null +++ b/libs/fpa/Image/EndPointsFilter.h @@ -0,0 +1,77 @@ +#ifndef __fpa__Image__EndPointsFilter__h__ +#define __fpa__Image__EndPointsFilter__h__ + +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TDistanceMap, class _TCostMap > + class EndPointsFilter + : public itk::Object + { + public: + typedef EndPointsFilter Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TDistanceMap TDistanceMap; + typedef _TCostMap TCostMap; + typedef typename TCostMap::IndexType TIndex; + typedef MinimumSpanningTree< TCostMap::ImageDimension > TMST; + + typedef itk::Functor::IndexLexicographicCompare< _TCostMap::ImageDimension > TIndexCompare; + typedef std::set< TIndex, TIndexCompare > TIndices; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::EndPointsFilter, itk::Object ); + + itkGetConstObjectMacro( DistanceMap, _TDistanceMap ); + itkGetConstObjectMacro( CostMap, _TCostMap ); + itkGetConstObjectMacro( MST, TMST ); + + itkSetConstObjectMacro( DistanceMap, _TDistanceMap ); + itkSetConstObjectMacro( CostMap, _TCostMap ); + itkSetConstObjectMacro( MST, TMST ); + + public: + const TIndices& GetBifurcations( ) const; + const TIndices& GetEndPoints( ) const; + + void Compute( ); + + protected: + EndPointsFilter( ); + virtual ~EndPointsFilter( ); + + private: + // Purposely not defined + EndPointsFilter( const Self& other ); + Self& operator=( const Self& other ); + + protected: + typename _TDistanceMap::ConstPointer m_DistanceMap; + typename _TCostMap::ConstPointer m_CostMap; + typename TMST::ConstPointer m_MST; + + TIndices m_Bifurcations; + TIndices m_EndPoints; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__EndPointsFilter__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/EndPointsFilter.hxx b/libs/fpa/Image/EndPointsFilter.hxx new file mode 100644 index 0000000..c2e1d1c --- /dev/null +++ b/libs/fpa/Image/EndPointsFilter.hxx @@ -0,0 +1,153 @@ +#ifndef __fpa__Image__EndPointsFilter__hxx__ +#define __fpa__Image__EndPointsFilter__hxx__ + +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +const typename fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >:: +TIndices& fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >:: +GetBifurcations( ) const +{ + return( this->m_Bifurcations ); +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +const typename fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >:: +TIndices& fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >:: +GetEndPoints( ) const +{ + return( this->m_EndPoints ); +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +void fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >:: +Compute( ) +{ + typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt; + typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt; + typedef std::multimap< double, TIndex, std::greater< double > > _TQueue; + typedef typename _TQueue::value_type _TQueueValue; + typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks; + + // Some values + typename _TDistanceMap::RegionType region = + this->m_DistanceMap->GetRequestedRegion( ); + typename _TMarks::Pointer marks = _TMarks::New( ); + marks->SetLargestPossibleRegion( this->m_DistanceMap->GetLargestPossibleRegion( ) ); + marks->SetRequestedRegion( this->m_DistanceMap->GetRequestedRegion( ) ); + marks->SetBufferedRegion( this->m_DistanceMap->GetBufferedRegion( ) ); + marks->SetSpacing( this->m_DistanceMap->GetSpacing( ) ); + marks->SetOrigin( this->m_DistanceMap->GetOrigin( ) ); + marks->SetDirection( this->m_DistanceMap->GetDirection( ) ); + marks->Allocate( ); + marks->FillBuffer( 0 ); + + // Create queue + _TQueue queue; + _TDistMapIt dIt( this->m_DistanceMap, region ); + _TCostMapIt cIt( this->m_CostMap, region ); + dIt.GoToBegin( ); + cIt.GoToBegin( ); + for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt ) + { + double d = double( dIt.Get( ) ); + if( d > 0 ) + { + double v = double( cIt.Get( ) ) / ( d * d ); + queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) ); + + } // fi + + } // rof + + // BFS from maximum queue + while( queue.size( ) > 0 ) + { + // Get node + auto nIt = queue.begin( ); + auto n = *nIt; + queue.erase( nIt ); + + unsigned char m = marks->GetPixel( n.second ); + if( ( m & 0x01 ) == 0x01 ) + continue; + + // Mark it + m |= 0x01; + marks->SetPixel( n.second, m ); + this->m_EndPoints.insert( n.second ); + + // Get path + typename TMST::TPath::Pointer path; + this->m_MST->GetPath( path, n.second ); + for( unsigned long i = 0; i < path->GetSize( ); ++i ) + { + TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i ); + typename _TCostMap::PointType cnt; + this->m_CostMap->TransformIndexToPhysicalPoint( idx, cnt ); + double d = double( this->m_DistanceMap->GetPixel( idx ) ); + d *= double( 2 ); + + // Mark sphere + std::queue< TIndex > q; + q.push( idx ); + while( q.size( ) > 0 ) + { + TIndex v = q.front( ); + q.pop( ); + unsigned char m = marks->GetPixel( v ); + if( ( m & 0x02 ) == 0x02 ) + continue; + m |= 0x03; + marks->SetPixel( v, m ); + + for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x ) + { + for( int y = -1; y <= 1; y += 2 ) + { + TIndex w = v; + w[ x ] += y; + if( region.IsInside( w ) ) + { + typename _TCostMap::PointType p; + this->m_CostMap->TransformIndexToPhysicalPoint( w, p ); + if( cnt.EuclideanDistanceTo( p ) <= d ) + q.push( w ); + + } // fi + + } // rof + + } // rof + + } // elihw + + } // rof + + } // elihw +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >:: +EndPointsFilter( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >:: +~EndPointsFilter( ) +{ +} + +#endif // __fpa__Image__EndPointsFilter__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/FastMarching.h b/libs/fpa/Image/FastMarching.h new file mode 100644 index 0000000..2a4056f --- /dev/null +++ b/libs/fpa/Image/FastMarching.h @@ -0,0 +1,55 @@ +#ifndef __fpa__Image__FastMarching__h__ +#define __fpa__Image__FastMarching__h__ + +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TInputImage, class _TOutputImage > + class FastMarching + : public fpa::Base::FastMarching< fpa::Image::Algorithm< _TInputImage, _TOutputImage > > + { + public: + typedef FastMarching Self; + typedef fpa::Image::Algorithm< _TInputImage, _TOutputImage > TAlgorithm; + typedef fpa::Base::FastMarching< TAlgorithm > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TFastMarchingNeighbor TFastMarchingNeighbor; + typedef typename Superclass::TFastMarchingNeighborhood TFastMarchingNeighborhood; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::FastMarching, fpa::Base::FastMarching ); + + protected: + FastMarching( ); + virtual ~FastMarching( ); + + virtual TFastMarchingNeighborhood _FastMarchingNeighbors( const TVertex& v ) const override; + + private: + // Purposely not defined + FastMarching( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__FastMarching__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/FastMarching.hxx b/libs/fpa/Image/FastMarching.hxx new file mode 100644 index 0000000..fe75b21 --- /dev/null +++ b/libs/fpa/Image/FastMarching.hxx @@ -0,0 +1,50 @@ +#ifndef __fpa__Image__FastMarching__hxx__ +#define __fpa__Image__FastMarching__hxx__ + +#include + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::FastMarching< _TInputImage, _TOutputImage >:: +FastMarching( ) + : Superclass( ) +{ + typedef fpa::Image::Functors::VertexCost< _TInputImage, TOutput > _TCost; + typename _TCost::Pointer cost = _TCost::New( ); + this->SetVertexFunction( cost ); + this->m_InitResult = std::numeric_limits< TOutput >::max( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::FastMarching< _TInputImage, _TOutputImage >:: +~FastMarching( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +typename fpa::Image::FastMarching< _TInputImage, _TOutputImage >:: +TFastMarchingNeighborhood +fpa::Image::FastMarching< _TInputImage, _TOutputImage >:: +_FastMarchingNeighbors( const TVertex& v ) const +{ + TFastMarchingNeighborhood neighs; + typename _TInputImage::SpacingType spac = this->GetInput( )->GetSpacing( ); + for( unsigned int d = 0; d < _TInputImage::ImageDimension; ++d ) + { + for( int i = -1; i <= 1; i += 2 ) + { + TVertex n = v; + n[ d ] += i; + neighs.push_back( TFastMarchingNeighbor( n, spac[ d ] ) ); + + } // rof + + } // rof + return( neighs ); +} + +#endif // __fpa__Image__FastMarching__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/Base.h b/libs/fpa/Image/Functors/Base.h new file mode 100644 index 0000000..45bd846 --- /dev/null +++ b/libs/fpa/Image/Functors/Base.h @@ -0,0 +1,56 @@ +#ifndef __fpa__Image__Functors__Base__h__ +#define __fpa__Image__Functors__Base__h__ + +#include +#include +#include + +namespace fpa +{ + namespace Image + { + namespace Functors + { + /** + */ + template< class _TImage, class _TSuperclass > + class Base + : public _TSuperclass + { + public: + typedef Base Self; + typedef _TSuperclass Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TImage TImage; + typedef itk::ImageBase< TImage::ImageDimension > TImageBase; + + public: + itkTypeMacro( Base, itk::FunctionBase ); + + itkGetConstObjectMacro( Image, TImage ); + itkSetConstObjectMacro( Image, TImage ); + + protected: + Base( ) : Superclass( ) { } + virtual ~Base( ) { } + + private: + // Purposely not implemented + Base( const Self& other ); + Self& operator=( const Self& other ); + + protected: + typename TImage::ConstPointer m_Image; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#endif // __fpa__Image__Functors__Base__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/RegionGrow/BinaryThreshold.h b/libs/fpa/Image/Functors/RegionGrow/BinaryThreshold.h new file mode 100644 index 0000000..70ac8d0 --- /dev/null +++ b/libs/fpa/Image/Functors/RegionGrow/BinaryThreshold.h @@ -0,0 +1,75 @@ +#ifndef __fpa__Image__Functors__RegionGrow__BinaryThreshold__h__ +#define __fpa__Image__Functors__RegionGrow__BinaryThreshold__h__ + +#include +#include + +namespace fpa +{ + namespace Image + { + namespace Functors + { + namespace RegionGrow + { + /** + */ + template< class _TImage, class _TOutput > + class BinaryThreshold + : public fpa::Image::Functors::Base< _TImage, fpa::Base::Functors::RegionGrow::Base< typename _TImage::IndexType, _TOutput > > + { + public: + typedef _TImage TImage; + typedef _TOutput TOutput; + typedef typename TImage::IndexType TIndex; + typedef typename TImage::PixelType TPixel; + + typedef fpa::Base::Functors::RegionGrow::Base< TIndex, TOutput > TBase; + typedef fpa::Image::Functors::Base< TImage, TBase > Superclass; + typedef BinaryThreshold Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( BinaryThreshold, Base ); + + itkGetConstMacro( Lower, TPixel ); + itkGetConstMacro( Upper, TPixel ); + itkSetMacro( Lower, TPixel ); + itkSetMacro( Upper, TPixel ); + + public: + virtual TOutput Evaluate( + const TIndex& a, const TIndex& b + ) const override; + + protected: + BinaryThreshold( ); + virtual ~BinaryThreshold( ); + + private: + // Purposely not implemented + BinaryThreshold( const Self& other ); + Self& operator=( const Self& other ); + + protected: + TPixel m_Lower; + TPixel m_Upper; + }; + + } // ecapseman + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__Functors__RegionGrow__BinaryThreshold__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/RegionGrow/BinaryThreshold.hxx b/libs/fpa/Image/Functors/RegionGrow/BinaryThreshold.hxx new file mode 100644 index 0000000..7f2b9dc --- /dev/null +++ b/libs/fpa/Image/Functors/RegionGrow/BinaryThreshold.hxx @@ -0,0 +1,49 @@ +#ifndef __fpa__Image__Functors__RegionGrow__BinaryThreshold__hxx__ +#define __fpa__Image__Functors__RegionGrow__BinaryThreshold__hxx__ + +#include + +// ------------------------------------------------------------------------- +template< class _TImage, class _TOutput > +_TOutput +fpa::Image::Functors::RegionGrow::BinaryThreshold< _TImage, _TOutput >:: +Evaluate( const TIndex& a, const TIndex& b ) const +{ + const _TImage* im = + dynamic_cast< const _TImage* >( this->m_Image.GetPointer( ) ); + if( im != NULL ) + { + TPixel v = im->GetPixel( b ); + return( + ( this->m_Lower < v && v < this->m_Upper )? + this->m_InsideValue: + this->m_OutsideValue + ); + } + else + return( this->m_OutsideValue ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TOutput > +fpa::Image::Functors::RegionGrow::BinaryThreshold< _TImage, _TOutput >:: +BinaryThreshold( ) + : Superclass( ) +{ + this->m_Upper = std::numeric_limits< TPixel >::max( ); + if( std::numeric_limits< TPixel >::is_integer ) + this->m_Lower = std::numeric_limits< TPixel >::min( ); + else + this->m_Lower = -this->m_Upper; +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TOutput > +fpa::Image::Functors::RegionGrow::BinaryThreshold< _TImage, _TOutput >:: +~BinaryThreshold( ) +{ +} + +#endif // __fpa__Image__Functors__RegionGrow__BinaryThreshold__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/SimpleNeighborhood.h b/libs/fpa/Image/Functors/SimpleNeighborhood.h new file mode 100644 index 0000000..0e4ebf7 --- /dev/null +++ b/libs/fpa/Image/Functors/SimpleNeighborhood.h @@ -0,0 +1,71 @@ +#ifndef __fpa__Image__Functors__SimpleNeighborhood__h__ +#define __fpa__Image__Functors__SimpleNeighborhood__h__ + +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + namespace Functors + { + /** + */ + template< unsigned int _VDim > + class SimpleNeighborhood + : public fpa::Image::Functors::Base< itk::ImageBase< _VDim >, itk::FunctionBase< itk::Index< _VDim >, std::vector< itk::Index< _VDim > > > > + { + public: + typedef itk::ImageBase< _VDim > TImage; + typedef typename TImage::IndexType TIndex; + typedef typename TIndex::OffsetType TOffset; + typedef std::vector< TIndex > TOutput; + typedef itk::FunctionBase< TIndex, TOutput > TBaseFunctor; + typedef fpa::Image::Functors::Base< TImage, TBaseFunctor > Superclass; + typedef SimpleNeighborhood Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( SimpleNeighborhood, Base ); + + itkGetConstMacro( Order, unsigned int ); + itkSetMacro( Order, unsigned int ); + + public: + virtual TOutput Evaluate( const TIndex& center ) const override; + + protected: + SimpleNeighborhood( ); + virtual ~SimpleNeighborhood( ); + + void _1stCombination( ) const; + void _2ndCombination( ) const; + + private: + // Purposely not implemented + SimpleNeighborhood( const Self& other ); + Self& operator=( const Self& other ); + + protected: + unsigned int m_Order; + mutable std::vector< TOffset > m_Offsets; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__Functors__SimpleNeighborhood__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/SimpleNeighborhood.hxx b/libs/fpa/Image/Functors/SimpleNeighborhood.hxx new file mode 100644 index 0000000..172425d --- /dev/null +++ b/libs/fpa/Image/Functors/SimpleNeighborhood.hxx @@ -0,0 +1,113 @@ +#ifndef __fpa__Image__Functors__SimpleNeighborhood__hxx__ +#define __fpa__Image__Functors__SimpleNeighborhood__hxx__ + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +typename fpa::Image::Functors::SimpleNeighborhood< _VDim >:: +TOutput fpa::Image::Functors::SimpleNeighborhood< _VDim >:: +Evaluate( const TIndex& center ) const +{ + if( this->m_Offsets.size( ) == 0 ) + { + if( this->m_Order == 1 ) + this->_1stCombination( ); + else + this->_2ndCombination( ); + + } // fi + + TOutput res; + typename TImage::RegionType reg = this->m_Image->GetRequestedRegion( ); + for( int i = 0; i < this->m_Offsets.size( ); ++i ) + { + TIndex idx = center + this->m_Offsets[ i ]; + if( reg.IsInside( idx ) && idx != center ) + res.push_back( idx ); + + } // rof + return( res ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +fpa::Image::Functors::SimpleNeighborhood< _VDim >:: +SimpleNeighborhood( ) + : Superclass( ), + m_Order( 1 ) +{ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > + fpa::Image::Functors::SimpleNeighborhood< _VDim >:: +~SimpleNeighborhood( ) +{ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Image::Functors::SimpleNeighborhood< _VDim >:: +_1stCombination( ) const +{ + for( int d = 0; d < TImage::ImageDimension; ++d ) + { + typename TIndex::OffsetType off; + off.Fill( 0 ); + for( int i = -1; i <= 1; i += 2 ) + { + off[ d ] = i; + this->m_Offsets.push_back( off ); + + } // rof + + } // rof +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Image::Functors::SimpleNeighborhood< _VDim >:: +_2ndCombination( ) const +{ + std::vector< std::vector< std::vector< int > > > M; + M.push_back( std::vector< std::vector< int > >( ) ); + + std::vector< int > base; + base.push_back( -1 ); + base.push_back( 0 ); + base.push_back( 1 ); + int dim = TImage::ImageDimension; + + M.push_back( std::vector< std::vector< int > >( ) ); + for( int j = 0; j < base.size( ); ++j ) + M[ 1 ].push_back( std::vector< int >( 1, base[ j ] ) ); + + for( int i = 2; i <= dim; ++i ) + { + M.push_back( std::vector< std::vector< int > >( ) ); + for( int j = 0; j < base.size( ); ++j ) + { + for( int k = 0; k < M[ i - 1 ].size( ); ++k ) + { + M[ i ].push_back( std::vector< int >( 1, base[ j ] ) ); + for( int l = 0; l < M[ i - 1 ][ k ].size( ); ++l ) + M[ i ].back( ).push_back( M[ i - 1 ][ k ][ l ] ); + + } // rof + + } // rof + + } // rof + + for( int i = 0; i < M[ dim ].size( ); ++i ) + { + TOffset off; + for( int j = 0; j < M[ dim ][ i ].size( ); ++j ) + off[ j ] = M[ dim ][ i ][ j ]; + this->m_Offsets.push_back( off ); + + } // rof +} + +#endif // __fpa__Image__Functors__SimpleNeighborhood__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/VertexCost.h b/libs/fpa/Image/Functors/VertexCost.h new file mode 100644 index 0000000..a6f3e45 --- /dev/null +++ b/libs/fpa/Image/Functors/VertexCost.h @@ -0,0 +1,68 @@ +#ifndef __fpa__Image__Functors__VertexCost__h__ +#define __fpa__Image__Functors__VertexCost__h__ + +#include +#include + +namespace fpa +{ + namespace Image + { + namespace Functors + { + /** + */ + template< class _TImage, class _TOutput > + class VertexCost + : public fpa::Image::Functors::Base< _TImage, fpa::Base::Functors::VertexCostFunctionBase< typename _TImage::IndexType, _TOutput > > + { + public: + typedef _TImage TImage; + typedef typename TImage::IndexType TIndex; + typedef _TOutput TOutput; + + typedef fpa::Base::Functors::VertexCostFunctionBase< TIndex, TOutput > TBaseFunctor; + typedef fpa::Image::Functors::Base< TImage, TBaseFunctor > Superclass; + typedef VertexCost Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( VertexCost, Base ); + + itkBooleanMacro( UseImageSpacing ); + itkGetConstMacro( UseImageSpacing, bool ); + itkSetMacro( UseImageSpacing, bool ); + + public: + virtual TOutput Evaluate( + const TIndex& a, const TIndex& b + ) const override; + + protected: + VertexCost( ); + virtual ~VertexCost( ); + + private: + // Purposely not implemented + VertexCost( const Self& other ); + Self& operator=( const Self& other ); + + protected: + bool m_UseImageSpacing; + }; + + } // ecapseman + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__Functors__VertexCost__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/Functors/VertexCost.hxx b/libs/fpa/Image/Functors/VertexCost.hxx new file mode 100644 index 0000000..3b59c52 --- /dev/null +++ b/libs/fpa/Image/Functors/VertexCost.hxx @@ -0,0 +1,42 @@ +#ifndef __fpa__Image__Functors__VertexCost__hxx__ +#define __fpa__Image__Functors__VertexCost__hxx__ + +// ------------------------------------------------------------------------- +template< class _TImage, class _TOutput > +typename fpa::Image::Functors::VertexCost< _TImage, _TOutput >:: +TOutput fpa::Image::Functors::VertexCost< _TImage, _TOutput >:: +Evaluate( const TIndex& a, const TIndex& b ) const +{ + const _TImage* im = + dynamic_cast< const _TImage* >( this->m_Image.GetPointer( ) ); + TOutput coeff = TOutput( 1 ); + if( this->m_UseImageSpacing ) + { + typename _TImage::PointType pa, pb; + im->TransformIndexToPhysicalPoint( a, pa ); + im->TransformIndexToPhysicalPoint( b, pb ); + coeff = TOutput( pa.EuclideanDistanceTo( pb ) ); + + } // fi + return( TOutput( im->GetPixel( a ) ) * coeff ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TOutput > +fpa::Image::Functors::VertexCost< _TImage, _TOutput >:: +VertexCost( ) + : Superclass( ), + m_UseImageSpacing( false ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TOutput > +fpa::Image::Functors::VertexCost< _TImage, _TOutput >:: +~VertexCost( ) +{ +} + +#endif // __fpa__Image__Functors__VertexCost__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/MinimumSpanningTree.h b/libs/fpa/Image/MinimumSpanningTree.h new file mode 100644 index 0000000..7d3941d --- /dev/null +++ b/libs/fpa/Image/MinimumSpanningTree.h @@ -0,0 +1,65 @@ +#ifndef __fpa__Image__MinimumSpanningTree__h__ +#define __fpa__Image__MinimumSpanningTree__h__ + +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< unsigned int _VDim > + class MinimumSpanningTree + : public fpa::Base::MinimumSpanningTree< itk::Index< _VDim >, cpExtensions::DataStructures::PolyLineParametricPath< _VDim >, itk::Image< itk::Offset< _VDim >, _VDim > > + { + public: + typedef itk::Index< _VDim > TVertex; + typedef cpExtensions::DataStructures::PolyLineParametricPath< _VDim > TPath; + typedef itk::Offset< _VDim > TOffset; + typedef itk::Image< TOffset, _VDim > TImage; + typedef fpa::Base::MinimumSpanningTree< TVertex, TPath, TImage > Superclass; + typedef MinimumSpanningTree Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( + fpa::Image::MinimumSpanningTree, fpa::Base::MinimumSpanningTree + ); + + public: + virtual TVertex GetParent( const TVertex& v ) const override; + virtual void SetParent( const TVertex& v, const TVertex& p ) override; + + virtual void GetPath( + typename TPath::Pointer& path, const TVertex& a + ) const override; + virtual void GetPath( + typename TPath::Pointer& path, const TVertex& a, const TVertex& b + ) const override; + + protected: + MinimumSpanningTree( ); + virtual ~MinimumSpanningTree( ); + + private: + // Purposely not defined + MinimumSpanningTree( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__MinimumSpanningTree__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/MinimumSpanningTree.hxx b/libs/fpa/Image/MinimumSpanningTree.hxx new file mode 100644 index 0000000..b99953b --- /dev/null +++ b/libs/fpa/Image/MinimumSpanningTree.hxx @@ -0,0 +1,63 @@ +#ifndef __fpa__Image__MinimumSpanningTree__hxx__ +#define __fpa__Image__MinimumSpanningTree__hxx__ + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +typename fpa::Image::MinimumSpanningTree< _VDim >:: +TVertex fpa::Image::MinimumSpanningTree< _VDim >:: +GetParent( const TVertex& v ) const +{ + TVertex p = v + this->GetPixel( v ); + return( p ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Image::MinimumSpanningTree< _VDim >:: +SetParent( const TVertex& v, const TVertex& p ) +{ + this->SetPixel( v, p - v ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Image::MinimumSpanningTree< _VDim >:: +GetPath( typename TPath::Pointer& path, const TVertex& a ) const +{ + if( path.GetPointer( ) == NULL ) + path = TPath::New( ); + path->SetReferenceImage( this ); + this->Superclass::GetPath( path, a ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Image::MinimumSpanningTree< _VDim >:: +GetPath( + typename TPath::Pointer& path, const TVertex& a, const TVertex& b + ) const +{ + if( path.GetPointer( ) == NULL ) + path = TPath::New( ); + path->SetReferenceImage( this ); + this->Superclass::GetPath( path, a, b ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +fpa::Image::MinimumSpanningTree< _VDim >:: +MinimumSpanningTree( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +fpa::Image::MinimumSpanningTree< _VDim >:: +~MinimumSpanningTree( ) +{ +} + +#endif // __fpa__Image__MinimumSpanningTree__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/MoriRegionGrow.h b/libs/fpa/Image/MoriRegionGrow.h new file mode 100644 index 0000000..41c122e --- /dev/null +++ b/libs/fpa/Image/MoriRegionGrow.h @@ -0,0 +1,94 @@ +#ifndef __fpa__Image__MoriRegionGrow__h__ +#define __fpa__Image__MoriRegionGrow__h__ + +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TInputImage, class _TOutputImage, class _TAuxPixel = unsigned short > + class MoriRegionGrow + : public itk::ImageToImageFilter< _TInputImage, _TOutputImage > + { + public: + typedef MoriRegionGrow Self; + typedef itk::ImageToImageFilter< _TInputImage, _TOutputImage > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef itk::Image< _TAuxPixel, _TInputImage::ImageDimension > TAuxImage; + typedef fpa::Image::MoriRegionGrowHelper< _TInputImage, TAuxImage > THelper; + typedef itk::BinaryThresholdImageFilter< TAuxImage, _TOutputImage > TThreshold; + + typedef typename _TInputImage::IndexType TIndex; + typedef typename _TInputImage::PixelType TInputPixel; + typedef typename _TOutputImage::PixelType TOutputPixel; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::MoriRegionGrow, itk::ImageToImageFilter ); + + itkGetConstMacro( Seed, TIndex ); + itkSetMacro( Seed, TIndex ); + + public: + TAuxImage* GetAuxiliaryImage( ); + const TAuxImage* GetAuxiliaryImage( ) const; + + TInputPixel GetLower( ) const; + TInputPixel GetUpper( ) const; + TInputPixel GetStep( ) const; + TOutputPixel GetInsideValue( ) const; + TOutputPixel GetOutsideValue( ) const; + + void SetLower( const TInputPixel& v ); + void SetUpper( const TInputPixel& v ); + void SetStep( const TInputPixel& v ); + void SetInsideValue( const TOutputPixel& v ); + void SetOutsideValue( const TOutputPixel& v ); + + + /* TODO + itkGetConstMacro( Lower, TPixel ); + itkGetConstMacro( Upper, TPixel ); + itkGetConstMacro( Step, TPixel ); + itkGetConstMacro( Sensitivity, double ); + itkSetMacro( Lower, TPixel ); + itkSetMacro( Upper, TPixel ); + itkSetMacro( Step, TPixel ); + itkSetMacro( Sensitivity, double ); + */ + + protected: + MoriRegionGrow( ); + virtual ~MoriRegionGrow( ); + + virtual void GenerateData( ) override; + + private: + // Purposely not defined + MoriRegionGrow( const Self& other ); + Self& operator=( const Self& other ); + + protected: + typename THelper::Pointer m_Helper; + typename TThreshold::Pointer m_Threshold; + TIndex m_Seed; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__MoriRegionGrow__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/MoriRegionGrow.hxx b/libs/fpa/Image/MoriRegionGrow.hxx new file mode 100644 index 0000000..737bf1f --- /dev/null +++ b/libs/fpa/Image/MoriRegionGrow.hxx @@ -0,0 +1,162 @@ +#ifndef __fpa__Image__MoriRegionGrow__hxx__ +#define __fpa__Image__MoriRegionGrow__hxx__ + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +typename +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +TAuxImage* +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +GetAuxiliaryImage( ) +{ + return( this->m_Helper->GetOutput( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +const typename +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +TAuxImage* +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +GetAuxiliaryImage( ) const +{ + return( this->m_Helper->GetOutput( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +typename +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +TInputPixel +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +GetLower( ) const +{ + return( this->m_Helper->GetLower( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +typename +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +TInputPixel +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +GetUpper( ) const +{ + return( this->m_Helper->GetUpper( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +typename +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +TInputPixel +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +GetStep( ) const +{ + return( this->m_Helper->GetStep( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +typename +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +TOutputPixel +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +GetInsideValue( ) const +{ + return( this->m_Threshold->GetInsideValue( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +typename +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +TOutputPixel +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +GetOutsideValue( ) const +{ + return( this->m_Threshold->GetInsideValue( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +SetLower( const TInputPixel& v ) +{ + this->m_Helper->SetLower( v ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +SetUpper( const TInputPixel& v ) +{ + this->m_Helper->SetUpper( v ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +SetStep( const TInputPixel& v ) +{ + this->m_Helper->SetStep( v ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +SetInsideValue( const TOutputPixel& v ) +{ + this->m_Threshold->SetInsideValue( v ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +SetOutsideValue( const TOutputPixel& v ) +{ + this->m_Threshold->SetOutsideValue( v ); + this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +MoriRegionGrow( ) + : Superclass( ) +{ + this->m_Helper = THelper::New( ); + this->m_Threshold = TThreshold::New( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +~MoriRegionGrow( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage, class _TAuxPixel > +void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxPixel >:: +GenerateData( ) +{ + this->m_Helper->ClearSeeds( ); + this->m_Helper->AddSeed( this->m_Seed, 0 ); + this->m_Helper->SetInput( this->GetInput( ) ); + this->m_Helper->Update( ); + + this->m_Threshold->SetInput( this->m_Helper->GetOutput( ) ); + this->m_Threshold->SetLowerThreshold( std::numeric_limits< _TAuxPixel >::min( ) ); + this->m_Threshold->SetUpperThreshold( this->m_Helper->GetOptimumThreshold( ) ); + this->m_Threshold->Update( ); + this->GetOutput( )->Graft( this->m_Threshold->GetOutput( ) ); +} + +#endif // __fpa__Image__MoriRegionGrow__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/MoriRegionGrowHelper.h b/libs/fpa/Image/MoriRegionGrowHelper.h new file mode 100644 index 0000000..02f51b6 --- /dev/null +++ b/libs/fpa/Image/MoriRegionGrowHelper.h @@ -0,0 +1,90 @@ +#ifndef __fpa__Image__MoriRegionGrowHelper__h__ +#define __fpa__Image__MoriRegionGrowHelper__h__ + +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TInputImage, class _TOutputImage > + class MoriRegionGrowHelper + : public fpa::Image::RegionGrow< _TInputImage, _TOutputImage > + { + public: + typedef MoriRegionGrowHelper Self; + typedef fpa::Image::RegionGrow< _TInputImage, _TOutputImage > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TGrowFunction TGrowFunction; + typedef typename _TInputImage::PixelType TPixel; + typedef fpa::Image::Functors::RegionGrow::BinaryThreshold< _TInputImage, TOutput > TBinThresholdFunction; + + protected: + typedef typename Superclass::_TQueueNode _TQueueNode; + typedef typename Superclass::_TQueue _TQueue; + + typedef std::pair< TPixel, unsigned long > TCurveData; + typedef std::vector< TCurveData > TCurve; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::MoriRegionGrowHelper, fpa::Image::RegionGrow ); + + itkGetConstMacro( Lower, TPixel ); + itkGetConstMacro( Upper, TPixel ); + itkGetConstMacro( Step, TPixel ); + itkGetConstMacro( OptimumThreshold, typename _TOutputImage::PixelType ); + + itkSetMacro( Lower, TPixel ); + itkSetMacro( Upper, TPixel ); + itkSetMacro( Step, TPixel ); + + + protected: + MoriRegionGrowHelper( ); + virtual ~MoriRegionGrowHelper( ); + + virtual bool _ContinueGenerateData( ) override; + virtual void _BeforeGenerateData( ) override; + virtual void _AfterGenerateData( ) override; + virtual void _BeforeLoop( ) override; + virtual void _AfterLoop( ) override; + virtual bool _UpdateValue( _TQueueNode& v, const _TQueueNode& p ) override; + virtual void _UpdateResult( const _TQueueNode& n ) override; + + private: + // Purposely not defined + MoriRegionGrowHelper( const Self& other ); + Self& operator=( const Self& other ); + + protected: + TPixel m_Lower; + TPixel m_Upper; + TPixel m_Step; + typename _TOutputImage::PixelType m_OptimumThreshold; + + _TQueue m_NextQueue; + unsigned long m_ActualCount; + TCurve m_Curve; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__MoriRegionGrowHelper__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/MoriRegionGrowHelper.hxx b/libs/fpa/Image/MoriRegionGrowHelper.hxx new file mode 100644 index 0000000..9dc1df3 --- /dev/null +++ b/libs/fpa/Image/MoriRegionGrowHelper.hxx @@ -0,0 +1,144 @@ +#ifndef __fpa__Image__MoriRegionGrowHelper__hxx__ +#define __fpa__Image__MoriRegionGrowHelper__hxx__ + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +MoriRegionGrowHelper( ) + : Superclass( ), + m_Step( TPixel( 1 ) ) +{ + this->m_Upper = std::numeric_limits< TPixel >::max( ); + if( std::numeric_limits< TPixel >::is_integer ) + this->m_Lower = std::numeric_limits< TPixel >::min( ); + else + this->m_Lower = -this->m_Upper; + typename TBinThresholdFunction::Pointer functor = + TBinThresholdFunction::New( ); + this->SetGrowFunction( functor ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +~MoriRegionGrowHelper( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +bool fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +_ContinueGenerateData( ) +{ + TBinThresholdFunction* functor = + dynamic_cast< TBinThresholdFunction* >( this->GetGrowFunction( ) ); + TPixel u = functor->GetUpper( ); + + // Update flooding data + this->m_Curve.push_back( TCurveData( u, this->m_ActualCount ) ); + + // Update thresholds + if( u < this->m_Upper ) + { + u += this->m_Step; + if( u > this->m_Upper ) + u = this->m_Upper; + functor->SetUpper( u ); + this->m_Queue = this->m_NextQueue; + while( this->m_NextQueue.size( ) > 0 ) + this->m_NextQueue.pop( ); + return( true ); + } + else + return( false ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +_BeforeGenerateData( ) +{ + this->Superclass::_BeforeGenerateData( ); + + while( this->m_NextQueue.size( ) > 0 ) + this->m_NextQueue.pop( ); + this->m_OptimumThreshold = ( typename _TOutputImage::PixelType )( 0 ); + this->m_ActualCount = 0; + this->m_Curve.clear( ); + TBinThresholdFunction* functor = + dynamic_cast< TBinThresholdFunction* >( this->GetGrowFunction( ) ); + functor->SetLower( this->m_Lower ); + functor->SetUpper( this->m_Lower + this->m_Step ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +_AfterGenerateData( ) +{ + typedef typename _TOutputImage::PixelType _TOut; + + this->Superclass::_AfterGenerateData( ); + while( this->m_NextQueue.size( ) > 0 ) + this->m_NextQueue.pop( ); + + // Find optimum threshold by dicotomy + unsigned long l = 0; + unsigned long r = this->m_Curve.size( ) - 1; + while( ( r - l ) > 1 ) + { + unsigned long m = ( r + l ) >> 1; + double vm = double( this->m_Curve[ m ].second ); + double dl = vm - double( this->m_Curve[ l ].second ); + double dr = double( this->m_Curve[ r ].second ) - vm; + if( dl > dr ) + r = m; + else + l = m; + + } // elihw + this->m_OptimumThreshold = _TOut( ( r + l ) >> 1 ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +_BeforeLoop( ) +{ + this->Superclass::_BeforeLoop( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +_AfterLoop( ) +{ + this->Superclass::_AfterLoop( ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +bool fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +_UpdateValue( _TQueueNode& v, const _TQueueNode& p ) +{ + typedef typename _TOutputImage::PixelType _TOut; + + bool ret = this->Superclass::_UpdateValue( v, p ); + v.Result = _TOut( this->m_Curve.size( ) + 1 ); + if( !ret ) + this->m_NextQueue.push( v ); + return( ret ); +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::MoriRegionGrowHelper< _TInputImage, _TOutputImage >:: +_UpdateResult( const _TQueueNode& n ) +{ + this->Superclass::_UpdateResult( n ); + this->m_ActualCount += 1; +} + +#endif // __fpa__Image__MoriRegionGrowHelper__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/RegionGrow.h b/libs/fpa/Image/RegionGrow.h new file mode 100644 index 0000000..dbc57f0 --- /dev/null +++ b/libs/fpa/Image/RegionGrow.h @@ -0,0 +1,55 @@ +#ifndef __fpa__Image__RegionGrow__h__ +#define __fpa__Image__RegionGrow__h__ + +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TInputImage, class _TOutputImage > + class RegionGrow + : public fpa::Base::RegionGrow< fpa::Image::Algorithm< _TInputImage, _TOutputImage > > + { + public: + typedef fpa::Image::Algorithm< _TInputImage, _TOutputImage > TAlgorithm; + typedef RegionGrow Self; + typedef fpa::Base::RegionGrow< TAlgorithm > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TVertex TVertex; + + typedef fpa::Image::Functors::Base< _TInputImage, typename Superclass::TGrowFunction > TGrowFunction; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::RegionGrow, fpa::Base::RegionGrow ); + + protected: + RegionGrow( ); + virtual ~RegionGrow( ); + + virtual void _BeforeGenerateData( ) override; + + private: + // Purposely not defined + RegionGrow( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__RegionGrow__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/RegionGrow.hxx b/libs/fpa/Image/RegionGrow.hxx new file mode 100644 index 0000000..326dc54 --- /dev/null +++ b/libs/fpa/Image/RegionGrow.hxx @@ -0,0 +1,34 @@ +#ifndef __fpa__Image__RegionGrow__hxx__ +#define __fpa__Image__RegionGrow__hxx__ + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::RegionGrow< _TInputImage, _TOutputImage >:: +RegionGrow( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +fpa::Image::RegionGrow< _TInputImage, _TOutputImage >:: +~RegionGrow( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TInputImage, class _TOutputImage > +void fpa::Image::RegionGrow< _TInputImage, _TOutputImage >:: +_BeforeGenerateData( ) +{ + this->Superclass::_BeforeGenerateData( ); + this->m_InitResult = this->GetOutsideValue( ); + TGrowFunction* grow = + dynamic_cast< TGrowFunction* >( this->GetGrowFunction( ) ); + if( grow != NULL ) + grow->SetImage( this->GetInput( ) ); +} + +#endif // __fpa__Image__RegionGrow__hxx__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/SkeletonFilter.h b/libs/fpa/Image/SkeletonFilter.h new file mode 100644 index 0000000..0fbb3d5 --- /dev/null +++ b/libs/fpa/Image/SkeletonFilter.h @@ -0,0 +1,83 @@ +#ifndef __fpa__Image__SkeletonFilter__h__ +#define __fpa__Image__SkeletonFilter__h__ + +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TImage > + class SkeletonFilter + : public Dijkstra< _TImage, _TImage > + { + public: + typedef SkeletonFilter Self; + typedef Dijkstra< _TImage, _TImage > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TImage TImage; + itkStaticConstMacro( + Dimension, unsigned int, TImage::ImageDimension + ); + + typedef typename Superclass::TMST TMST; + typedef typename TImage::IndexType TIndex; + typedef itk::Image< unsigned char, Self::Dimension > TMarks; + typedef cpExtensions::DataStructures::Skeleton< Self::Dimension > TSkeleton; + + protected: + typedef typename Superclass::_TQueueNode _TQueueNode; + typedef std::multimap< double, TIndex, std::greater< double > > _TSkeletonQueue; + typedef std::map< TIndex, TIndex, typename TIndex::LexicographicCompare > _TAdjacencies; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::SkeletonFilter, itk::Object ); + + public: + TSkeleton* GetSkeleton( ); + TMarks* GetMarks( ); + + protected: + SkeletonFilter( ); + virtual ~SkeletonFilter( ); + + virtual void _BeforeGenerateData( ) override; + virtual void _UpdateResult( const _TQueueNode& n ) override; + virtual void _AfterGenerateData( ) override; + + void _EndPoints( std::vector< TIndex >& end_points, _TAdjacencies& A ); + void _Skeleton( const std::vector< TIndex >& end_points, _TAdjacencies& A ); + void _MarkSphere( const TIndex& idx ); + + private: + // Purposely not defined + SkeletonFilter( const Self& other ); + Self& operator=( const Self& other ); + + protected: + _TSkeletonQueue m_SkeletonQueue; + + unsigned long m_SkeletonIdx; + unsigned long m_MarksIdx; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__SkeletonFilter__h__ + +// eof - $RCSfile$ diff --git a/libs/fpa/Image/SkeletonFilter.hxx b/libs/fpa/Image/SkeletonFilter.hxx new file mode 100644 index 0000000..c60e13c --- /dev/null +++ b/libs/fpa/Image/SkeletonFilter.hxx @@ -0,0 +1,274 @@ +#ifndef __fpa__Image__SkeletonFilter__hxx__ +#define __fpa__Image__SkeletonFilter__hxx__ + +#include +#include +#include + +// ------------------------------------------------------------------------- +#define fpa_Image_SkeletonFilter_OutputMacro( o_n, o_t ) \ + template< class _TImage > \ + typename fpa::Image::SkeletonFilter< _TImage >:: \ + o_t* fpa::Image::SkeletonFilter< _TImage >:: \ + Get##o_n( ) \ + { \ + return( \ + dynamic_cast< o_t* >( \ + this->itk::ProcessObject::GetOutput( this->m_##o_n##Idx ) \ + ) \ + ); \ + } + +fpa_Image_SkeletonFilter_OutputMacro( Skeleton, TSkeleton ); +fpa_Image_SkeletonFilter_OutputMacro( Marks, TMarks ); + +// ------------------------------------------------------------------------- +template< class _TImage > +fpa::Image::SkeletonFilter< _TImage >:: +SkeletonFilter( ) + : Superclass( ) +{ + typedef typename _TImage::PixelType _TPixel; + typedef fpa::Image::Functors::SimpleNeighborhood< _TImage::ImageDimension > _TNeighFunc; + typedef fpa::Base::Functors::Inverse< _TPixel, _TPixel > _TInvFunc; + + unsigned int nOutputs = this->GetNumberOfRequiredOutputs( ); + this->SetNumberOfRequiredOutputs( nOutputs + 2 ); + this->m_SkeletonIdx = nOutputs; + this->m_MarksIdx = nOutputs + 1; + + typename TSkeleton::Pointer skeleton = TSkeleton::New( ); + typename TMarks::Pointer marks = TMarks::New( ); + this->SetNthOutput( this->m_SkeletonIdx, skeleton.GetPointer( ) ); + this->SetNthOutput( this->m_MarksIdx, marks.GetPointer( ) ); + + typename _TNeighFunc::Pointer nfunc = _TNeighFunc::New( ); + nfunc->SetOrder( 2 ); + this->SetNeighborhoodFunction( nfunc ); + + typename _TInvFunc::Pointer ifunc = _TInvFunc::New( ); + this->SetConversionFunction( ifunc ); +} + +// ------------------------------------------------------------------------- +template< class _TImage > +fpa::Image::SkeletonFilter< _TImage >:: +~SkeletonFilter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void fpa::Image::SkeletonFilter< _TImage >:: +_BeforeGenerateData( ) +{ + this->Superclass::_BeforeGenerateData( ); + this->m_SkeletonQueue.clear( ); +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void fpa::Image::SkeletonFilter< _TImage >:: +_UpdateResult( const _TQueueNode& n ) +{ + typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue; + + this->Superclass::_UpdateResult( n ); + double d = double( this->GetInput( )->GetPixel( n.Vertex ) ); + if( d >= double( 0 ) ) + { + // Update skeleton candidates + d += double( 1e-5 ); + double v = double( n.Result ) / ( d * d ); + this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void fpa::Image::SkeletonFilter< _TImage >:: +_AfterGenerateData( ) +{ + this->Superclass::_AfterGenerateData( ); + + _TAdjacencies A; + std::vector< TIndex > end_points; + this->_EndPoints( end_points, A ); + this->_Skeleton( end_points, A ); +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void fpa::Image::SkeletonFilter< _TImage >:: +_EndPoints( std::vector< TIndex >& end_points, _TAdjacencies& A ) +{ + auto marks = this->GetMarks( ); + auto mst = this->GetMinimumSpanningTree( ); + auto spac = marks->GetSpacing( ); + + // Some values + marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) ); + marks->SetRequestedRegion( mst->GetRequestedRegion( ) ); + marks->SetBufferedRegion( mst->GetBufferedRegion( ) ); + marks->SetSpacing( mst->GetSpacing( ) ); + marks->SetOrigin( mst->GetOrigin( ) ); + marks->SetDirection( mst->GetDirection( ) ); + marks->Allocate( ); + marks->FillBuffer( 0 ); + + // BFS from maximum queue + while( this->m_SkeletonQueue.size( ) > 0 ) + { + // Get node + auto nIt = this->m_SkeletonQueue.begin( ); + auto n = *nIt; + this->m_SkeletonQueue.erase( nIt ); + + // Mark it and update end-points + unsigned char m = marks->GetPixel( n.second ); + if( m != 0 ) + continue; + marks->SetPixel( n.second, 1 ); + end_points.push_back( n.second ); + + // Mark path + TIndex it = n.second; + TIndex p = mst->GetParent( it ); + while( it != p ) + { + this->_MarkSphere( it ); + it = p; + p = mst->GetParent( it ); + + } // elihw + this->_MarkSphere( it ); + A[ n.second ] = it; + + } // elihw +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void fpa::Image::SkeletonFilter< _TImage >:: +_Skeleton( const std::vector< TIndex >& end_points, _TAdjacencies& A ) +{ + typedef itk::Image< unsigned long, Self::Dimension > _TTagsImage; + typedef typename TMST::TPath _TPath; + + auto mst = this->GetMinimumSpanningTree( ); + auto skeleton = this->GetSkeleton( ); + + // Tag branches + typename _TTagsImage::Pointer tags = _TTagsImage::New( ); + tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) ); + tags->SetRequestedRegion( mst->GetRequestedRegion( ) ); + tags->SetBufferedRegion( mst->GetBufferedRegion( ) ); + tags->Allocate( ); + tags->FillBuffer( 0 ); + for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt ) + { + TIndex it = *eIt; + TIndex p = mst->GetParent( it ); + while( it != p ) + { + tags->SetPixel( it, tags->GetPixel( it ) + 1 ); + it = p; + p = mst->GetParent( it ); + + } // elihw + tags->SetPixel( it, tags->GetPixel( it ) + 1 ); + + } // rof + + // Build paths (branches) + for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt ) + { + TIndex it = *eIt; + TIndex p = mst->GetParent( it ); + TIndex sIdx = *eIt; + typename _TPath::Pointer path = _TPath::New( ); + path->SetReferenceImage( mst ); + while( it != p ) + { + if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) ) + { + // Ok, a new branch has been added + path->AddVertex( it ); + skeleton->AddBranch( path ); + + // Update a new starting index + path = _TPath::New( ); + path->SetReferenceImage( mst ); + sIdx = it; + } + else + path->AddVertex( it ); + it = p; + p = mst->GetParent( it ); + + } // elihw + + // Finally, add last branch + path->AddVertex( it ); + skeleton->AddBranch( path ); + + } // rof +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void fpa::Image::SkeletonFilter< _TImage >:: +_MarkSphere( const TIndex& idx ) +{ + typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt; + + static const double _eps = std::sqrt( double( Self::Dimension + 1 ) ); + auto input = this->GetInput( ); + auto marks = this->GetMarks( ); + auto spac = input->GetSpacing( ); + auto region = input->GetRequestedRegion( ); + + typename _TImage::PointType cnt; + input->TransformIndexToPhysicalPoint( idx, cnt ); + double r = double( input->GetPixel( idx ) ) * _eps; + + TIndex i0, i1; + for( unsigned int d = 0; d < Self::Dimension; ++d ) + { + long off = long( std::ceil( r / double( spac[ d ] ) ) ); + if( off < 3 ) + off = 3; + i0[ d ] = idx[ d ] - off; + i1[ d ] = idx[ d ] + off; + + if( i0[ d ] < region.GetIndex( )[ d ] ) + i0[ d ] = region.GetIndex( )[ d ]; + + if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] ) + i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1; + + } // rof + + typename _TImage::SizeType size; + for( unsigned int d = 0; d < Self::Dimension; ++d ) + size[ d ] = i1[ d ] - i0[ d ] + 1; + + typename _TImage::RegionType neighRegion; + neighRegion.SetIndex( i0 ); + neighRegion.SetSize( size ); + + _TMarksIt mIt( marks, neighRegion ); + for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) + { + TIndex w = mIt.GetIndex( ); + typename _TImage::PointType p; + marks->TransformIndexToPhysicalPoint( w, p ); + mIt.Set( 1 ); + + } // rof +} + +#endif // __fpa__Image__SkeletonFilter__hxx__ + +// eof - $RCSfile$