From cd6b3f8433deaf2ba7c6bb0ece8e5912c760db17 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Sun, 14 May 2017 20:51:22 -0500 Subject: [PATCH] ... --- examples/CMakeLists.txt | 1 + examples/SkeletonFilter.cxx | 125 +++++++++ lib/fpa/Base/Graph.h | 170 ++++++++++++ lib/fpa/Base/Graph.hxx | 272 +++++++++++++++++++ lib/fpa/Base/PolyLineParametricPath.h | 97 +++++++ lib/fpa/Base/PolyLineParametricPath.hxx | 194 ++++++++++++++ lib/fpa/Base/SeedsInterface.h | 2 +- lib/fpa/Base/Skeleton.h | 67 +++++ lib/fpa/Base/Skeleton.hxx | 99 +++++++ lib/fpa/Image/SkeletonFilter.h | 96 +++++++ lib/fpa/Image/SkeletonFilter.hxx | 336 ++++++++++++++++++++++++ 11 files changed, 1458 insertions(+), 1 deletion(-) create mode 100644 examples/SkeletonFilter.cxx create mode 100644 lib/fpa/Base/Graph.h create mode 100644 lib/fpa/Base/Graph.hxx create mode 100644 lib/fpa/Base/PolyLineParametricPath.h create mode 100644 lib/fpa/Base/PolyLineParametricPath.hxx create mode 100644 lib/fpa/Base/Skeleton.h create mode 100644 lib/fpa/Base/Skeleton.hxx create mode 100644 lib/fpa/Image/SkeletonFilter.h create mode 100644 lib/fpa/Image/SkeletonFilter.hxx diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 4feb556..72d9ad2 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -6,6 +6,7 @@ if(BUILD_EXAMPLES) RegionGrow_BinaryThreshold RegionGrow_Mori Dijkstra_Maurer + SkeletonFilter #CreateMoriInputImage #BronchiiInitialSegmentationWithMori #BronchiiInitialSegmentationWithBinaryThresholdRegionGrow diff --git a/examples/SkeletonFilter.cxx b/examples/SkeletonFilter.cxx new file mode 100644 index 0000000..25d54c4 --- /dev/null +++ b/examples/SkeletonFilter.cxx @@ -0,0 +1,125 @@ + +#include +#include +#include +#include + +#include + +/* TODO + #include + #include + #include + #include +*/ + +// ------------------------------------------------------------------------- +static const unsigned int VDim = 3; +typedef unsigned char TPixel; +typedef double TScalar; +typedef itk::Image< TPixel, VDim > TImage; +typedef itk::ImageFileReader< TImage > TReader; +typedef fpa::Image::SkeletonFilter< TImage, TScalar > TFilter; + +/* TODO + typedef itk::ImageFileWriter< TImage > TWriter; + typedef fpa::Image::MoriRegionGrow< TImage, TImage > TFilter; + typedef itk::BinaryThresholdImageFilter< TImage, TImage > TThreshold; +*/ + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + // Get arguments + if( argc != 3 && argc != 3 + VDim ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " input_image output_skeleton ["; + for( unsigned int i = 0; i < VDim; ++i ) + std::cerr << " s_" << i; + std::cerr << " ]" << std::endl; + return( 1 ); + + } // fi + std::string input_image_filename = argv[ 1 ]; + std::string output_skeleton_filename = argv[ 2 ]; + + TReader::Pointer reader = TReader::New( ); + reader->SetFileName( input_image_filename ); + try + { + reader->Update( ); + } + catch( std::exception& err ) + { + std::cerr << "ERROR: " << err.what( ) << std::endl; + return( 1 ); + + } // yrt + + TFilter::Pointer filter = TFilter::New( ); + filter->SetInput( reader->GetOutput( ) ); + + /* TODO + filter->SetThresholdRange( lower, upper, delta ); + TImage::PointType pnt; + for( int i = 0; i < VDim; ++i ) + pnt[ i ] = std::atof( argv[ i + 6 ] ); + + TImage::IndexType seed; + if( !( reader->GetOutput( )->TransformPhysicalPointToIndex( pnt, seed ) ) ) + { + std::cerr << "ERROR: seed outside image." << std::endl; + return( 1 ); + + } // fi + filter->AddSeed( seed ); + + std::chrono::time_point< std::chrono::system_clock > start, end; + start = std::chrono::system_clock::now( ); + filter->Update( ); + end = std::chrono::system_clock::now( ); + std::chrono::duration< double > elapsed_seconds = end - start; + + TThreshold::Pointer threshold = TThreshold::New( ); + threshold->SetInput( filter->GetOutput( ) ); + threshold->SetInsideValue( 255 ); + threshold->SetOutsideValue( 0 ); + threshold->SetLowerThreshold( 0 ); + threshold->SetUpperThreshold( filter->GetOptimumThreshold( ) ); + + TWriter::Pointer writer = TWriter::New( ); + writer->SetInput( threshold->GetOutput( ) ); + writer->SetFileName( output_image_filename ); + try + { + writer->Update( ); + } + catch( std::exception& err ) + { + std::cerr << "ERROR: " << err.what( ) << std::endl; + return( 1 ); + + } // yrt + + // Show data + TFilter::TCurve curve = filter->GetCurve( ); + for( TFilter::TCurveData data: curve ) + { + std::cout << data.XValue << " " << data.YValue << " " << data.Diff1 << std::endl; + } + std::cout + << std::endl + << "# Opt: " + << curve[ filter->GetOptimumThreshold( ) ].XValue + << "(" + << filter->GetOptimumThreshold( ) + << ")" + << std::endl; + std::cout << "Time: " << elapsed_seconds.count( ) << "s" << std::endl; + */ + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Graph.h b/lib/fpa/Base/Graph.h new file mode 100644 index 0000000..41b2fb9 --- /dev/null +++ b/lib/fpa/Base/Graph.h @@ -0,0 +1,170 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__Graph__h__ +#define __fpa__Base__Graph__h__ + +#include +#include +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** \brief A generic graph with templated index types. + * + * @param _TVertex Vertex type. + * @param _TCost Cost type. + * @param _TIndex Index type (it should be a strict weak ordering type). + */ + template< class _TVertex, class _TCost, class _TIndex = unsigned long, class _TIndexCompare = std::less< _TIndex > > + class Graph + : public itk::DataObject + { + public: + typedef Graph Self; + typedef itk::DataObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TVertex TVertex; + typedef _TCost TCost; + typedef _TIndex TIndex; + typedef _TIndexCompare TIndexCompare; + + // Base types + typedef std::map< TIndex, TVertex, TIndexCompare > TVertices; + typedef std::vector< TCost > TEdges; + typedef std::map< TIndex, TEdges, TIndexCompare > TMatrixRow; + typedef std::map< TIndex, TMatrixRow, TIndexCompare > TMatrix; + + public: + itkNewMacro( Self ); + itkTypeMacro( Graph, itk::DataObject ); + + public: + /*! \brief Iterators over vertices. + * These allow you to iterate over all of graph's vertices. + * + * Typical iteration should be done as: + * + * TGraph g; + * ... + * TGraph::TVertices::[const_]iterator vIt = g.BeginVertices( ); + * for( ; vIt != g.EndVertices( ); ++vIt ) + * { + * vIt->first; --> this is the vertex's index <-- + * vIt->second; --> this is the vertex's value <-- + * } + */ + inline typename TVertices::iterator BeginVertices( ) + { return( this->m_Vertices.begin( ) ); } + inline typename TVertices::iterator EndVertices( ) + { return( this->m_Vertices.end( ) ); } + inline typename TVertices::const_iterator BeginVertices( ) const + { return( this->m_Vertices.begin( ) ); } + inline typename TVertices::const_iterator EndVertices( ) const + { return( this->m_Vertices.end( ) ); } + + /*! \brief Iterators over edges. + * These allow you to iterate over all of graph's edges. + * + * Typical iteration should be done as: + * + * TGraph g; + * ... + * TGraph::TMatrix::[const_]iterator mIt = g.BeginEdgesRows( ); + * for( ; mIt != g.EndEdgesRows( ); ++mIt ) + * { + * mIt->first; --> this is the row index. <-- + * TGraph::TMatrixRow::[const_]iterator rIt = mIt->second.begin( ); + * for( ; rIt != mIt->second.end( ); ++rIt ) + * { + * rIt->first; --> this is the column index. + * TGraph::TEdges::[const_]iterator eIt = rIt->second.begin( ); + * for( ; eIt != rIt->second.end( ); ++eIt ) + * *eIt; --> this is the cost between mIt->first and rIt->first + * } + * } + */ + inline typename TMatrix::iterator BeginEdgesRows( ) + { return( this->m_Matrix.begin( ) ); } + inline typename TMatrix::iterator EndEdgetsRows( ) + { return( this->m_Matrix.end( ) ); } + inline typename TMatrix::const_iterator BeginEdgesRows( ) const + { return( this->m_Matrix.begin( ) ); } + inline typename TMatrix::const_iterator EndEdgesRows( ) const + { return( this->m_Matrix.end( ) ); } + + /*! \brief Clear all vertices and edges. + */ + void Clear( ); + + /*! \brief Clear all edges. + */ + inline void ClearEdges( ) + { this->m_Matrix.clear( ); } + + /*! \brief Vertex manipulation methods. + * Names are self-explanatory. + */ + inline bool HasVertexIndex( const TIndex& i ) const + { return( this->m_Vertices.find( i ) != this->m_Vertices.end( ) ); } + inline void SetVertex( const TIndex& index, TVertex& vertex ) + { this->m_Vertices[ index ] = vertex; } + inline TVertex& GetVertex( const TIndex& index ) + { return( this->m_Vertices[ index ] ); } + inline const TVertex& GetVertex( const TIndex& index ) const + { return( this->m_Vertices[ index ] ); } + bool RenameVertex( const TIndex& old_index, const TIndex& new_index ); + void RemoveVertex( const TIndex& index ); + + /*! \brief Edge manipulation methods. + * Names are self-explanatory. + */ + inline void AddEdge( const TIndex& orig, const TIndex& dest, const TCost& cost ) + { this->m_Matrix[ orig ][ dest ].push_back( cost ); } + TEdges& GetEdges( const TIndex& orig, const TIndex& dest ); + const TEdges& GetEdges( const TIndex& orig, const TIndex& dest ) const; + bool HasEdge( const TIndex& orig, const TIndex& dest ) const; + void RemoveEdge( const TIndex& orig, const TIndex& dest, const TCost& cost ); + void RemoveEdges( const TIndex& orig, const TIndex& dest ); + + /*! \brief Returns graph's sinks. + * + * A sink is a special vertex which does not have any "exiting" edges. + * + * @return Sinks ordered by their index. + */ + std::set< TIndex, TIndexCompare > GetSinks( ) const; + + protected: + Graph( ); + virtual ~Graph( ); + + private: + // Purposely not implemented + Graph( const Self& other ); + Self& operator=( const Self& other ); + + protected: + TVertices m_Vertices; + TMatrix m_Matrix; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__Graph__h__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Graph.hxx b/lib/fpa/Base/Graph.hxx new file mode 100644 index 0000000..89fec7b --- /dev/null +++ b/lib/fpa/Base/Graph.hxx @@ -0,0 +1,272 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__Graph__hxx__ +#define __fpa__Base__Graph__hxx__ + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +void fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +Clear( ) +{ + this->m_Vertices.clear( ); + this->m_Matrix.clear( ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +bool fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +RenameVertex( const TIndex& old_index, const TIndex& new_index ) +{ + auto old_v = this->m_Vertices.find( old_index ); + auto new_v = this->m_Vertices.find( new_index ); + if( old_v != this->m_Vertices.end( ) && new_v == this->m_Vertices.end( ) ) + { + // Replace vertex + this->m_Vertices[ new_index ] = old_v->second; + this->m_Vertices.erase( old_index ); + + // Duplicate edges + auto mIt = this->m_Matrix.begin( ); + auto found_row = this->m_Matrix.end( ); + for( ; mIt != this->m_Matrix.end( ); ++mIt ) + { + if( mIt->first == old_index ) + found_row = mIt; + + auto rIt = mIt->second.begin( ); + for( ; rIt != mIt->second.end( ); ++rIt ) + { + if( mIt->first == old_index ) + this->m_Matrix[ new_index ][ rIt->first ] = rIt->second; + else if( rIt->first == old_index ) + this->m_Matrix[ mIt->first ][ new_index ] = rIt->second; + + } // rof + + } // rof + + // Delete old edges + if( found_row != this->m_Matrix.end( ) ) + this->m_Matrix.erase( found_row ); + + mIt = this->m_Matrix.begin( ); + for( ; mIt != this->m_Matrix.end( ); ++mIt ) + { + auto rIt = mIt->second.begin( ); + while( rIt != mIt->second.end( ) ) + { + if( rIt->first == old_index ) + { + mIt->second.erase( rIt ); + rIt = mIt->second.begin( ); + } + else + ++rIt; + + } // elihw + + } // rof + return( true ); + } + else + return( false ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +void fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +RemoveVertex( const TIndex& index ) +{ + auto i = this->m_Vertices.find( index ); + if( i != this->m_Vertices.end( ) ) + { + // Delete vertex + this->m_Vertices.erase( i ); + + // Delete edges starting from given vertex + auto mIt = this->m_Matrix.find( index ); + if( mIt != this->m_Matrix.end( ) ) + this->m_Matrix.erase( mIt ); + + // Delete edges arriving to given vertex + mIt = this->m_Matrix.begin( ); + for( ; mIt != this->m_Matrix.end( ); ++mIt ) + { + auto rIt = mIt->second.begin( ); + while( rIt != mIt->second.end( ) ) + { + if( rIt->first == index ) + { + mIt->second.erase( rIt ); + rIt = mIt->second.begin( ); + } + else + ++rIt; + + } // elihw + + } // rof + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +typename +fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +TEdges& +fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +GetEdges( const TIndex& orig, const TIndex& dest ) +{ + static TEdges null_edges; + auto o = this->m_Matrix.find( orig ); + if( o != this->m_Matrix.find( orig ) ) + { + auto d = o->second.find( dest ); + if( d == o->second.end( ) ) + { + null_edges.clear( ); + return( null_edges ); + } + else + return( d->second ); + } + else + { + null_edges.clear( ); + return( null_edges ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +const typename +fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +TEdges& +fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +GetEdges( const TIndex& orig, const TIndex& dest ) const +{ + static const TEdges null_edges; + auto o = this->m_Matrix.find( orig ); + if( o != this->m_Matrix.find( orig ) ) + { + auto d = o->second.find( dest ); + if( d == o->second.end( ) ) + return( null_edges ); + else + return( d->second ); + } + else + return( null_edges ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +bool fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +HasEdge( const TIndex& orig, const TIndex& dest ) const +{ + auto mIt = this->m_Matrix.find( orig ); + if( mIt != this->m_Matrix.end( ) ) + return( mIt->second.find( dest ) != mIt->second.end( ) ); + else + return( false ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +void fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +RemoveEdge( const TIndex& orig, const TIndex& dest, const TCost& cost ) +{ + auto m = this->m_Matrix.find( orig ); + if( m != this->m_Matrix.end( ) ) + { + auto r = m->second.find( dest ); + if( r != m->second.end( ) ) + { + auto e = r->second.end( ); + for( + auto i = r->second.begin( ); + i != r->second.end( ) && e == r->second.end( ); + ++i + ) + if( *i == cost ) + e = i; + if( e != r->second.end( ) ) + { + r->second.erase( e ); + if( r->second.size( ) == 0 ) + { + m->second.erase( r ); + if( m->second.size( ) == 0 ) + this->m_Matrix.erase( m ); + + } // fi + + } // fi + + } // fi + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +void fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +RemoveEdges( const TIndex& orig, const TIndex& dest ) +{ + auto m = this->m_Matrix.find( orig ); + if( m != this->m_Matrix.end( ) ) + { + auto r = m->second.find( dest ); + if( r != m->second.end( ) ) + { + m->second.erase( r ); + if( m->second.size( ) == 0 ) + this->m_Matrix.erase( m ); + + } // fi + + } // fi + +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +std::set< _TIndex, _TIndexCompare > +fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +GetSinks( ) const +{ + std::set< _TIndex, _TIndexCompare > sinks; + + auto vIt = this->m_Vertices.begin( ); + for( ; vIt != this->m_Vertices.end( ); ++vIt ) + sinks.insert( vIt->first ); + auto mIt = this->m_Matrix.begin( ); + for( ; mIt != this->m_Matrix.end( ); ++mIt ) + sinks.erase( mIt->first ); + + return( sinks ); +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +Graph( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TVertex, class _TCost, class _TIndex, class _TIndexCompare > +fpa::Base::Graph< _TVertex, _TCost, _TIndex, _TIndexCompare >:: +~Graph( ) +{ +} + +#endif // __fpa__Base__Graph__hxx__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/PolyLineParametricPath.h b/lib/fpa/Base/PolyLineParametricPath.h new file mode 100644 index 0000000..79eb530 --- /dev/null +++ b/lib/fpa/Base/PolyLineParametricPath.h @@ -0,0 +1,97 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__PolyLineParametricPath__h__ +#define __fpa__Base__PolyLineParametricPath__h__ + +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< unsigned int _VDim > + class PolyLineParametricPath + : public itk::PolyLineParametricPath< _VDim > + { + public: + typedef PolyLineParametricPath Self; + typedef itk::PolyLineParametricPath< _VDim > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef itk::ImageBase< _VDim > TImageBase; + typedef typename TImageBase::SpacingType TSpacing; + typedef typename TImageBase::PointType TPoint; + typedef typename TImageBase::DirectionType TDirection; + typedef typename Superclass::ContinuousIndexType TContinuousIndex; + typedef typename TContinuousIndex::IndexType TIndex; + + public: + itkNewMacro( Self ); + itkTypeMacro( PolyLineParametricPath, itk::PolyLineParametricPath ); + + itkGetConstReferenceMacro( Spacing, TSpacing ); + itkGetConstReferenceMacro( Origin, TPoint ); + itkGetConstReferenceMacro( Direction, TDirection ); + itkGetConstReferenceMacro( InverseDirection, TDirection ); + + itkSetMacro( Origin, TPoint ); + + public: + unsigned long GetSize( ) const; + TContinuousIndex GetContinuousVertex( unsigned long i ) const; + TIndex GetVertex( unsigned long i ) const; + TPoint GetPoint( unsigned long i ) const; + + virtual void SetSpacing( const TSpacing& spac ); + virtual void SetSpacing( const double spac[ _VDim ] ); + virtual void SetSpacing( const float spac[ _VDim ] ); + virtual void SetOrigin( const double ori[ _VDim ] ); + virtual void SetOrigin( const float ori[ _VDim ] ); + virtual void SetDirection( const TDirection& dir ); + + template< class _TRefImage > + inline void SetReferenceImage( const _TRefImage* image ) + { + this->SetSpacing( image->GetSpacing( ) ); + this->SetOrigin( image->GetOrigin( ) ); + this->SetDirection( image->GetDirection( ) ); + } + + protected: + PolyLineParametricPath( ); + virtual ~PolyLineParametricPath( ); + + virtual void _ComputeIndexToPhysicalPointMatrices( ); + + private: + // Purposely not implemented + PolyLineParametricPath( const Self& other ); + Self& operator=( const Self& other ); + + protected: + TSpacing m_Spacing; + TPoint m_Origin; + TDirection m_Direction; + TDirection m_InverseDirection; + TDirection m_IndexToPhysicalPoint; + TDirection m_PhysicalPointToIndex; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__PolyLineParametricPath__h__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/PolyLineParametricPath.hxx b/lib/fpa/Base/PolyLineParametricPath.hxx new file mode 100644 index 0000000..18732d5 --- /dev/null +++ b/lib/fpa/Base/PolyLineParametricPath.hxx @@ -0,0 +1,194 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__PolyLineParametricPath__hxx__ +#define __fpa__Base__PolyLineParametricPath__hxx__ + +#include +#include + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +unsigned long fpa::Base::PolyLineParametricPath< _VDim >:: +GetSize( ) const +{ + return( this->GetVertexList( )->Size( ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +typename fpa::Base::PolyLineParametricPath< _VDim >:: +TContinuousIndex +fpa::Base::PolyLineParametricPath< _VDim >:: +GetContinuousVertex( unsigned long i ) const +{ + return( this->GetVertexList( )->GetElement( i ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +typename fpa::Base::PolyLineParametricPath< _VDim >:: +TIndex fpa::Base::PolyLineParametricPath< _VDim >:: +GetVertex( unsigned long i ) const +{ + TContinuousIndex cidx = this->GetContinuousVertex( i ); + TIndex idx; + for( unsigned int d = 0; d < _VDim; ++d ) + idx[ d ] = cidx[ d ]; + return( idx ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +typename fpa::Base::PolyLineParametricPath< _VDim >:: +TPoint fpa::Base::PolyLineParametricPath< _VDim >:: +GetPoint( unsigned long i ) const +{ + typedef typename TPoint::CoordRepType _TCoordRep; + TPoint pnt; + TContinuousIndex idx = this->GetVertex( i ); + for( unsigned int r = 0; r < _VDim; ++r ) + { + _TCoordRep sum = itk::NumericTraits< _TCoordRep >::ZeroValue( ); + for( unsigned int c = 0; c < _VDim; ++c ) + sum += this->m_IndexToPhysicalPoint( r, c ) * idx[ c ]; + pnt[ r ] = sum + this->m_Origin[ r ]; + + } // rof + return( pnt ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Base::PolyLineParametricPath< _VDim >:: +SetSpacing( const TSpacing& spac ) +{ + if( this->m_Spacing != spac ) + { + this->m_Spacing = spac; + this->_ComputeIndexToPhysicalPointMatrices( ); + this->Modified( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Base::PolyLineParametricPath< _VDim >:: +SetSpacing( const double spac[ _VDim ] ) +{ + this->SetSpacing( TSpacing( spac ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Base::PolyLineParametricPath< _VDim >:: +SetSpacing( const float spac[ _VDim ] ) +{ + TSpacing s; + for( unsigned int d = 0; d < _VDim; ++d ) + s[ d ] = spac[ d ]; + this->SetSpacing( s ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Base::PolyLineParametricPath< _VDim >:: +SetOrigin( const double ori[ _VDim ] ) +{ + this->SetOrigin( TPoint( ori ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Base::PolyLineParametricPath< _VDim >:: +SetOrigin( const float ori[ _VDim ] ) +{ + this->SetOrigin( TPoint( ori ) ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Base::PolyLineParametricPath< _VDim >:: +SetDirection( const TDirection& dir ) +{ + bool modified = false; + for( unsigned int r = 0; r < _VDim; r++ ) + { + for( unsigned int c = 0; c < _VDim; c++ ) + { + if( + itk::Math::NotExactlyEquals( + this->m_Direction[ r ][ c ], dir[ r ][ c ] + ) + ) + { + this->m_Direction[ r ][ c ] = dir[ r ][ c ]; + modified = true; + } // fi + + } // rof + + } // rof + if( modified ) + { + this->_ComputeIndexToPhysicalPointMatrices( ); + this->m_InverseDirection = this->m_Direction.GetInverse( ); + this->Modified( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +fpa::Base::PolyLineParametricPath< _VDim >:: +PolyLineParametricPath( ) + : Superclass( ) +{ + this->m_Spacing.Fill( 1.0 ); + this->m_Origin.Fill( 0.0 ); + this->m_Direction.SetIdentity( ); + this->m_InverseDirection.SetIdentity( ); + this->m_IndexToPhysicalPoint.SetIdentity( ); + this->m_PhysicalPointToIndex.SetIdentity( ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +fpa::Base::PolyLineParametricPath< _VDim >:: +~PolyLineParametricPath( ) +{ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Base::PolyLineParametricPath< _VDim >:: +_ComputeIndexToPhysicalPointMatrices( ) +{ + TDirection scale; + scale.Fill( 0.0 ); + for( unsigned int i = 0; i < _VDim; i++ ) + { + if( this->m_Spacing[ i ] == 0.0 ) + itkExceptionMacro( + "A spacing of 0 is not allowed: Spacing is " << this->m_Spacing + ); + scale[ i ][ i ] = this->m_Spacing[ i ]; + + } // rof + + if( vnl_determinant( this->m_Direction.GetVnlMatrix( ) ) == 0.0 ) + itkExceptionMacro( + << "Bad direction, determinant is 0. Direction is " + << this->m_Direction + ); + this->m_IndexToPhysicalPoint = this->m_Direction * scale; + this->m_PhysicalPointToIndex = this->m_IndexToPhysicalPoint.GetInverse( ); + this->Modified( ); +} + +#endif // __fpa__Base__PolyLineParametricPath__hxx__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/SeedsInterface.h b/lib/fpa/Base/SeedsInterface.h index a994125..cfc1bd5 100644 --- a/lib/fpa/Base/SeedsInterface.h +++ b/lib/fpa/Base/SeedsInterface.h @@ -38,7 +38,7 @@ namespace fpa SeedsInterface( itk::ProcessObject* filter ); virtual ~SeedsInterface( ); - private: + protected: TSeeds m_Seeds; itk::ProcessObject* m_Filter; }; diff --git a/lib/fpa/Base/Skeleton.h b/lib/fpa/Base/Skeleton.h new file mode 100644 index 0000000..223d74e --- /dev/null +++ b/lib/fpa/Base/Skeleton.h @@ -0,0 +1,67 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__Skeleton__h__ +#define __fpa__Base__Skeleton__h__ + +#include +#include +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< unsigned int _VDim > + class Skeleton + : public fpa::Base::Graph< typename fpa::Base::PolyLineParametricPath< _VDim >::TIndex, typename fpa::Base::PolyLineParametricPath< _VDim >::Pointer, typename fpa::Base::PolyLineParametricPath< _VDim >::TIndex, typename fpa::Base::PolyLineParametricPath< _VDim >::TIndex::LexicographicCompare > + { + public: + typedef fpa::Base::PolyLineParametricPath< _VDim > TPath; + typedef typename TPath::TIndex TIndex; + typedef typename TIndex::LexicographicCompare TIndexCompare; + typedef typename TPath::Pointer TPathPointer; + + itkStaticConstMacro( Dimension, unsigned int, _VDim ); + + typedef fpa::Base::Graph< TIndex, TPathPointer, TIndex, TIndexCompare > Superclass; + typedef Skeleton Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( Skeleton, fpa::Base::Graph ); + + public: + void AddBranch( TPath* path ); + const TPath* GetBranch( const TIndex& a, const TIndex& b ) const; + + std::vector< TIndex > GetEndPoints( ) const; + std::vector< TIndex > GetBifurcations( ) const; + + protected: + Skeleton( ); + virtual ~Skeleton( ); + + private: + // Purposely not implemented + Skeleton( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__Skeleton__h__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/Skeleton.hxx b/lib/fpa/Base/Skeleton.hxx new file mode 100644 index 0000000..91ec404 --- /dev/null +++ b/lib/fpa/Base/Skeleton.hxx @@ -0,0 +1,99 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__Skeleton__hxx__ +#define __fpa__Base__Skeleton__hxx__ + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +void fpa::Base::Skeleton< _VDim >:: +AddBranch( TPath* path ) +{ + // Check inputs + if( path == NULL ) + return; + unsigned long size = path->GetSize( ); + if( size == 0 ) + return; + TIndex a = path->GetVertex( 0 ); + TIndex b = path->GetVertex( size - 1 ); + if( this->HasEdge( a, b ) ) + return; + + // Add path + this->SetVertex( a, a ); + this->SetVertex( b, b ); + this->AddEdge( a, b, path ); + this->AddEdge( b, a, path ); + // TODO: this->Modified( ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +const typename fpa::Base::Skeleton< _VDim >:: +TPath* fpa::Base::Skeleton< _VDim >:: +GetBranch( const TIndex& a, const TIndex& b ) const +{ + static const TPath* null_path = NULL; + if( this->HasEdge( a, b ) ) + return( this->GetEdges( a, b ).front( ) ); + else + return( null_path ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +std::vector< typename fpa::Base::Skeleton< _VDim >::TIndex > +fpa::Base::Skeleton< _VDim >:: +GetEndPoints( ) const +{ + std::vector< TIndex > res; + auto mIt = this->BeginEdgesRows( ); + for( ; mIt != this->EndEdgesRows( ); ++mIt ) + { + unsigned long count = mIt->second.size( ); + if( count == 1 ) + res.push_back( mIt->first ); + + } // rof + return( res ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +std::vector< typename fpa::Base::Skeleton< _VDim >::TIndex > +fpa::Base::Skeleton< _VDim >:: +GetBifurcations( ) const +{ + std::vector< TIndex > res; + auto mIt = this->BeginEdgesRows( ); + for( ; mIt != this->EndEdgesRows( ); ++mIt ) + { + unsigned long count = mIt->second.size( ); + if( count > 1 ) + res.push_back( mIt->first ); + + } // rof + return( res ); +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +fpa::Base::Skeleton< _VDim >:: +Skeleton( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< unsigned int _VDim > +fpa::Base::Skeleton< _VDim >:: +~Skeleton( ) +{ +} + +#endif // __fpa__Base__Skeleton__hxx__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/SkeletonFilter.h b/lib/fpa/Image/SkeletonFilter.h new file mode 100644 index 0000000..9f0542f --- /dev/null +++ b/lib/fpa/Image/SkeletonFilter.h @@ -0,0 +1,96 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__SkeletonFilter__h__ +#define __fpa__Image__SkeletonFilter__h__ + +#include + +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TImage, class _TScalar = double > + class SkeletonFilter + : public fpa::Image::Dijkstra< itk::Image< _TScalar, _TImage::ImageDimension >, itk::Image< _TScalar, _TImage::ImageDimension > > + { + public: + // Smart pointers + typedef _TImage TImage; + typedef _TScalar TScalar; + typedef itk::Image< TScalar, TImage::ImageDimension > TScalarImage; + typedef SkeletonFilter Self; + typedef fpa::Image::Dijkstra< TScalarImage, TScalarImage > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TOutputValue TOutputValue; + + typedef itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage > TDistanceMap; + typedef itk::Image< unsigned char, TImage::ImageDimension > TMarks; + typedef fpa::Base::Skeleton< TImage::ImageDimension > TSkeleton; + + protected: + typedef std::multimap< TScalar, TVertex, std::greater< TScalar > > _TSkeletonQueue; + typedef std::map< TVertex, TVertex, typename TVertex::LexicographicCompare > _TAdjacencies; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::SkeletonFilter, fpa::Image::Dijkstra ); + + itkBooleanMacro( SeedFromMaximumDistance ); + itkGetConstMacro( SeedFromMaximumDistance, bool ); + itkSetMacro( SeedFromMaximumDistance, bool ); + + public: + virtual void SetInput( const TScalarImage* image ) override; + void SetInput( const TImage* image ); + TImage* GetInput( ); + const TImage* GetInput( ) const; + + TSkeleton* GetSkeleton( ); + TMarks* GetMarks( ); + + protected: + SkeletonFilter( ); + virtual ~SkeletonFilter( ); + + virtual void GenerateData( ) override; + virtual void _SetOutputValue( const TVertex& vertex, const TOutputValue& value ) override; + + void _EndPoints( std::vector< TVertex >& end_points, _TAdjacencies& A ); + void _Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A ); + void _MarkSphere( const TVertex& idx ); + + private: + SkeletonFilter( const Self& other ); + Self& operator=( const Self& other ); + + protected: + bool m_SeedFromMaximumDistance; + typename TDistanceMap::Pointer m_DistanceMap; + _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/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx new file mode 100644 index 0000000..2b1cb65 --- /dev/null +++ b/lib/fpa/Image/SkeletonFilter.hxx @@ -0,0 +1,336 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__SkeletonFilter__hxx__ +#define __fpa__Image__SkeletonFilter__hxx__ + +#include + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +SetInput( const TScalarImage* image ) +{ + // Do nothing +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +SetInput( const TImage* image ) +{ + this->m_DistanceMap->SetInput( image ); + this->Superclass::SetInput( this->m_DistanceMap->GetOutput( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +_TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +GetInput( ) +{ + return( this->m_DistanceMap->GetInput( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +const _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +GetInput( ) const +{ + return( this->m_DistanceMap->GetInput( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +typename fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +TSkeleton* fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +GetSkeleton( ) +{ + return( + dynamic_cast< TSkeleton* >( + this->itk::ProcessObject::GetOutput( this->m_SkeletonIdx ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +typename fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +TMarks* fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +GetMarks( ) +{ + return( + dynamic_cast< TMarks* >( + this->itk::ProcessObject::GetOutput( this->m_MarksIdx ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +SkeletonFilter( ) + : Superclass( ), + m_SeedFromMaximumDistance( false ) +{ + this->m_DistanceMap = TDistanceMap::New( ); + this->m_DistanceMap->InsideIsPositiveOn( ); + this->m_DistanceMap->SquaredDistanceOff( ); + this->m_DistanceMap->UseImageSpacingOn( ); + + 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( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +~SkeletonFilter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +GenerateData( ) +{ + // Prepare distance map + this->m_DistanceMap->Update( ); + + // Update start seed + TVertex seed; + if( this->m_SeedFromMaximumDistance ) + { + typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax; + typename _TMinMax::Pointer minmax = _TMinMax::New( ); + minmax->SetImage( this->m_DistanceMap->GetOutput( ) ); + minmax->Compute( ); + seed = minmax->GetIndexOfMaximum( ); + } + else + seed = *( this->BeginSeeds( ) ); + this->m_Seeds.clear( ); + this->m_Seeds.insert( seed ); + + // Prepare skeleton candidates queue + this->m_SkeletonQueue.clear( ); + + // Go! + this->Superclass::GenerateData( ); + + // Backtracking + _TAdjacencies A; + std::vector< TVertex > end_points; + this->_EndPoints( end_points, A ); + this->_Skeleton( end_points, A ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TScalar > +void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +_SetOutputValue( const TVertex& vertex, const TOutputValue& value ) +{ + /* TODO + 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, class _TScalar > +void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +_EndPoints( std::vector< TVertex >& end_points, _TAdjacencies& A ) +{ + /* TODO + 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 + TVertex it = n.second; + TVertex 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, class _TScalar > +void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +_Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A ) +{ + /* TODO + 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 ) + { + TVertex it = *eIt; + TVertex 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 ) + { + TVertex it = *eIt; + TVertex p = mst->GetParent( it ); + TVertex 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, class _TScalar > +void fpa::Image::SkeletonFilter< _TImage, _TScalar >:: +_MarkSphere( const TVertex& idx ) +{ + /* TODO + 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; + + TVertex 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 ) + { + TVertex w = mIt.GetIndex( ); + typename _TImage::PointType p; + marks->TransformIndexToPhysicalPoint( w, p ); + mIt.Set( 1 ); + + } // rof + */ +} + +#endif // __fpa__Image__SkeletonFilter__hxx__ + +// eof - $RCSfile$ -- 2.47.1