RegionGrow_BinaryThreshold
RegionGrow_Mori
Dijkstra_Maurer
+ SkeletonFilter
#CreateMoriInputImage
#BronchiiInitialSegmentationWithMori
#BronchiiInitialSegmentationWithBinaryThresholdRegionGrow
--- /dev/null
+
+#include <chrono>
+#include <iomanip>
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+
+#include <fpa/Image/SkeletonFilter.h>
+
+/* TODO
+ #include <itkCommand.h>
+ #include <itkImageFileWriter.h>
+ #include <itkBinaryThresholdImageFilter.h>
+ #include <fpa/Image/MoriRegionGrow.h>
+*/
+
+// -------------------------------------------------------------------------
+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$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__Graph__h__
+#define __fpa__Base__Graph__h__
+
+#include <map>
+#include <set>
+#include <vector>
+#include <itkDataObject.h>
+#include <itkObjectFactory.h>
+
+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 <fpa/Base/Graph.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __fpa__Base__Graph__h__
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @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$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__PolyLineParametricPath__h__
+#define __fpa__Base__PolyLineParametricPath__h__
+
+#include <itkPolyLineParametricPath.h>
+#include <itkImageBase.h>
+
+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 <fpa/Base/PolyLineParametricPath.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __fpa__Base__PolyLineParametricPath__h__
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__PolyLineParametricPath__hxx__
+#define __fpa__Base__PolyLineParametricPath__hxx__
+
+#include <itkMath.h>
+#include <itkNumericTraits.h>
+
+// -------------------------------------------------------------------------
+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$
SeedsInterface( itk::ProcessObject* filter );
virtual ~SeedsInterface( );
- private:
+ protected:
TSeeds m_Seeds;
itk::ProcessObject* m_Filter;
};
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__Skeleton__h__
+#define __fpa__Base__Skeleton__h__
+
+#include <vector>
+#include <fpa/Base/Graph.h>
+#include <fpa/Base/PolyLineParametricPath.h>
+
+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 <fpa/Base/Skeleton.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __fpa__Base__Skeleton__h__
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @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$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__SkeletonFilter__h__
+#define __fpa__Image__SkeletonFilter__h__
+
+#include <fpa/Image/Dijkstra.h>
+
+#include <map>
+#include <itkSignedMaurerDistanceMapImageFilter.h>
+#include <fpa/Base/Skeleton.h>
+
+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 <fpa/Image/SkeletonFilter.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __fpa__Image__SkeletonFilter__h__
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__SkeletonFilter__hxx__
+#define __fpa__Image__SkeletonFilter__hxx__
+
+#include <itkMinimumMaximumImageCalculator.h>
+
+// -------------------------------------------------------------------------
+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$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff