From: Leonardo Flórez-Valencia Date: Mon, 15 May 2017 12:34:13 +0000 (-0500) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=FrontAlgorithms.git;a=commitdiff_plain;h=e8548d1529b25ab74bca646333612fbe3e16e73f ... --- diff --git a/examples/SkeletonFilter.cxx b/examples/SkeletonFilter.cxx index 25d54c4..4e1b316 100644 --- a/examples/SkeletonFilter.cxx +++ b/examples/SkeletonFilter.cxx @@ -4,28 +4,17 @@ #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; -*/ +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; +typedef fpa::Base::SkeletonWriter< TFilter::TSkeleton > TWriter; // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) @@ -61,64 +50,46 @@ int main( int argc, char* argv[] ) 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 ); + if( argc > 3 ) + { + filter->SeedFromMaximumDistanceOff( ); + TImage::PointType pnt; + for( int i = 0; i < VDim; ++i ) + pnt[ i ] = std::atof( argv[ i + 3 ] ); + TImage::IndexType seed; + if( !( reader->GetOutput( )->TransformPhysicalPointToIndex( pnt, seed ) ) ) + { + std::cerr << "ERROR: seed outside image." << std::endl; + return( 1 ); + + } // fi + filter->AddSeed( seed ); + } + else + filter->SeedFromMaximumDistanceOn( ); + + 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; + + TWriter::Pointer writer = TWriter::New( ); + writer->SetInput( filter->GetSkeleton( ) ); + writer->SetFileName( output_skeleton_filename ); + try + { + writer->Update( ); + } + catch( std::exception& err ) + { + std::cerr << "ERROR: " << err.what( ) << std::endl; + return( 1 ); - } // yrt + } // 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; - */ + // Show data + std::cout << "Time: " << elapsed_seconds.count( ) << "s" << std::endl; return( 0 ); } diff --git a/lib/fpa/Base/SkeletonReader.h b/lib/fpa/Base/SkeletonReader.h new file mode 100644 index 0000000..ea1eaa3 --- /dev/null +++ b/lib/fpa/Base/SkeletonReader.h @@ -0,0 +1,83 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__SkeletonReader__h__ +#define __fpa__Base__SkeletonReader__h__ + +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TSkeleton > + class SkeletonReader + : public itk::ProcessObject + { + public: + // Basic types + typedef SkeletonReader Self; + typedef itk::ProcessObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TSkeleton TSkeleton; + + public: + itkNewMacro( Self ); + itkTypeMacro( SkeletonReader, itk::ProcessObject ); + + itkGetConstMacro( FileName, std::string ); + itkSetMacro( FileName, std::string ); + + public: + TSkeleton* GetOutput( ); + TSkeleton* GetOutput( unsigned int i ); + + virtual void GraftOutput( itk::DataObject* out ); + virtual void GraftOutput( + const typename Superclass::DataObjectIdentifierType& key, + itk::DataObject* out + ); + virtual void GraftNthOutput( unsigned int i, itk::DataObject* out ); + virtual itk::DataObject::Pointer MakeOutput( + itk::ProcessObject::DataObjectPointerArraySizeType i + ) override; + + virtual void Update( ) override + { this->GenerateData( ); } + + protected: + SkeletonReader( ); + virtual ~SkeletonReader( ); + + virtual void GenerateData( ) override; + + // Do nothing + virtual void GenerateOutputInformation( ) override + { } + + private: + // Purposely not implemented + SkeletonReader( const Self& ); + void operator=( const Self& ); + + protected: + std::string m_FileName; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__SkeletonReader__h__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/SkeletonReader.hxx b/lib/fpa/Base/SkeletonReader.hxx new file mode 100644 index 0000000..1a59958 --- /dev/null +++ b/lib/fpa/Base/SkeletonReader.hxx @@ -0,0 +1,176 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__SkeletonReader__hxx__ +#define __fpa__Base__SkeletonReader__hxx__ + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +_TSkeleton* fpa::Base::SkeletonReader< _TSkeleton >:: +GetOutput( ) +{ + return( + itkDynamicCastInDebugMode< TSkeleton* >( this->GetPrimaryOutput( ) ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +_TSkeleton* fpa::Base::SkeletonReader< _TSkeleton >:: +GetOutput( unsigned int i ) +{ + return( + itkDynamicCastInDebugMode< TSkeleton* >( + this->itk::ProcessObject::GetOutput( i ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +void fpa::Base::SkeletonReader< _TSkeleton >:: +GraftOutput( itk::DataObject* out ) +{ + this->GraftNthOutput( 0, out ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +void fpa::Base::SkeletonReader< _TSkeleton >:: +GraftOutput( + const typename Superclass::DataObjectIdentifierType& key, + itk::DataObject* out + ) +{ + if( out == NULL ) + { + itkExceptionMacro( + << "Requested to graft output that is a NULL pointer" + ); + + } // fi + itk::DataObject* output = this->itk::ProcessObject::GetOutput( key ); + output->Graft( out ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +void fpa::Base::SkeletonReader< _TSkeleton >:: +GraftNthOutput( unsigned int i, itk::DataObject* out ) +{ + if( i >= this->GetNumberOfIndexedOutputs( ) ) + { + itkExceptionMacro( + << "Requested to graft output " << i + << " but this filter only has " + << this->GetNumberOfIndexedOutputs( ) + << " indexed Outputs." + ); + + } // fi + this->GraftOutput( this->MakeNameFromOutputIndex( i ), out ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +itk::DataObject::Pointer +fpa::Base::SkeletonReader< _TSkeleton >:: +MakeOutput( itk::ProcessObject::DataObjectPointerArraySizeType i ) +{ + return( TSkeleton::New( ).GetPointer( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +fpa::Base::SkeletonReader< _TSkeleton >:: +SkeletonReader( ) + : Superclass( ) +{ + typename TSkeleton::Pointer out = + static_cast< TSkeleton* >( this->MakeOutput( 0 ).GetPointer( ) ); + this->itk::ProcessObject::SetNumberOfRequiredInputs( 0 ); + this->itk::ProcessObject::SetNumberOfRequiredOutputs( 1 ); + this->itk::ProcessObject::SetNthOutput( 0, out.GetPointer( ) ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +fpa::Base::SkeletonReader< _TSkeleton >:: +~SkeletonReader( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +void fpa::Base::SkeletonReader< _TSkeleton >:: +GenerateData( ) +{ + typedef typename TSkeleton::TPath _TPath; + typedef typename _TPath::TSpacing _TSpacing; + typedef typename _TPath::TPoint _TPoint; + typedef typename _TPath::TDirection _TDirection; + typedef typename _TPath::TContinuousIndex _TContinuousIndex; + + std::string buffer; + /* TODO + if( !( fpa::Read( buffer, this->m_FileName ) ) ) + { + itkExceptionMacro( + << "Error reading skeleton from \"" << this->m_FileName << "\"" + ); + return; + } // fi + */ + + std::istringstream in( buffer ); + unsigned int dim; + in >> dim; + if( dim != TSkeleton::Dimension ) + { + itkExceptionMacro( + << "Mismatched skeletons dimension: " << dim + << " != " << TSkeleton::Dimension + ); + return; + + } // fi + + TSkeleton* out = this->GetOutput( ); + unsigned long size; + in >> size; + while( size > 0 ) + { + _TSpacing spa; + _TPoint ori; + _TDirection dir; + for( unsigned int d = 0; d < dim; ++d ) + in >> spa[ d ]; + for( unsigned int d = 0; d < dim; ++d ) + in >> ori[ d ]; + for( unsigned int d = 0; d < dim; ++d ) + for( unsigned int e = 0; e < dim; ++e ) + in >> dir[ d ][ e ]; + + typename _TPath::Pointer path = _TPath::New( ); + path->SetSpacing( spa ); + path->SetOrigin( ori ); + path->SetDirection( dir ); + for( unsigned long s = 0; s < size; ++s ) + { + _TContinuousIndex idx; + for( unsigned int d = 0; d < dim; ++d ) + in >> idx[ d ]; + path->AddVertex( idx ); + + } // rof + out->AddBranch( path ); + in >> size; + + } // elihw +} + +#endif // __fpa__Base__SkeletonReader__hxx__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/SkeletonWriter.h b/lib/fpa/Base/SkeletonWriter.h new file mode 100644 index 0000000..ab07bba --- /dev/null +++ b/lib/fpa/Base/SkeletonWriter.h @@ -0,0 +1,67 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__SkeletonWriter__h__ +#define __fpa__Base__SkeletonWriter__h__ + +#include + +namespace fpa +{ + namespace Base + { + /** + */ + template< class _TSkeleton > + class SkeletonWriter + : public itk::ProcessObject + { + public: + // Basic types + typedef SkeletonWriter Self; + typedef itk::ProcessObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TSkeleton TSkeleton; + + public: + itkNewMacro( Self ); + itkTypeMacro( SkeletonWriter, itk::ProcessObject ); + + itkGetConstMacro( FileName, std::string ); + itkSetMacro( FileName, std::string ); + + public: + const TSkeleton* GetInput( ) const; + void SetInput( const TSkeleton* skeleton ); + virtual void Update( ) override; + + protected: + SkeletonWriter( ); + virtual ~SkeletonWriter( ); + + virtual void GenerateData( ) override; + + private: + // Purposely not implemented + SkeletonWriter( const Self& ); + void operator=( const Self& ); + + protected: + std::string m_FileName; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Base__SkeletonWriter__h__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Base/SkeletonWriter.hxx b/lib/fpa/Base/SkeletonWriter.hxx new file mode 100644 index 0000000..e99d9ce --- /dev/null +++ b/lib/fpa/Base/SkeletonWriter.hxx @@ -0,0 +1,160 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Base__SkeletonWriter__hxx__ +#define __fpa__Base__SkeletonWriter__hxx__ + +#include + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +const _TSkeleton* fpa::Base::SkeletonWriter< _TSkeleton >:: +GetInput( ) const +{ + return( + dynamic_cast< const TSkeleton* >( + this->itk::ProcessObject::GetInput( 0 ) + ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +void fpa::Base::SkeletonWriter< _TSkeleton >:: +SetInput( const _TSkeleton* skeleton ) +{ + this->itk::ProcessObject::SetNthInput( + 0, const_cast< TSkeleton* >( skeleton ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +void fpa::Base::SkeletonWriter< _TSkeleton >:: +Update( ) +{ + TSkeleton* input = const_cast< TSkeleton* >( this->GetInput( ) ); + if( input != NULL ) + { + input->UpdateOutputInformation( ); + input->UpdateOutputData( ); + this->GenerateData( ); + this->ReleaseInputs( ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +fpa::Base::SkeletonWriter< _TSkeleton >:: +SkeletonWriter( ) + : Superclass( ), + m_FileName( "" ) +{ + this->SetNumberOfRequiredInputs( 1 ); +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +fpa::Base::SkeletonWriter< _TSkeleton >:: +~SkeletonWriter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TSkeleton > +void fpa::Base::SkeletonWriter< _TSkeleton >:: +GenerateData( ) +{ + const TSkeleton* sk = this->GetInput( ); + auto mIt = sk->BeginEdgesRows( ); + auto rIt = mIt->second.begin( ); + auto eIt = rIt->second.begin( ); + auto path = *eIt; + + // Write base information + std::stringstream out1, out2; + out1 << TSkeleton::Dimension << std::endl; + auto spa = path->GetSpacing( ); + for( unsigned int d = 0; d < TSkeleton::Dimension; ++d ) + out1 << spa[ d ] << " "; + out1 << std::endl; + auto dir = path->GetDirection( ); + for( unsigned int d = 0; d < TSkeleton::Dimension; ++d ) + for( unsigned int e = 0; e < TSkeleton::Dimension; ++e ) + out1 << dir[ d ][ e ] << " "; + out1 << std::endl; + auto ori = path->GetOrigin( ); + for( unsigned int d = 0; d < TSkeleton::Dimension; ++d ) + out1 << ori[ d ] << " "; + out1 << std::endl; + + // End points + auto end_points = sk->GetEndPoints( ); + out1 << end_points.size( ) << std::endl; + for( auto point : end_points ) + { + for( unsigned int d = 0; d < TSkeleton::Dimension; ++d ) + out1 << point[ d ] << " "; + out1 << std::endl; + + } // rof + + // Bifurcations + auto bifurcations = sk->GetBifurcations( ); + out1 << bifurcations.size( ) << std::endl; + for( auto point : bifurcations ) + { + for( unsigned int d = 0; d < TSkeleton::Dimension; ++d ) + out1 << point[ d ] << " "; + out1 << std::endl; + + } // rof + + // Write paths + unsigned long pathCount = 0; + mIt = sk->BeginEdgesRows( ); + for( ; mIt != sk->EndEdgesRows( ); ++mIt ) + { + auto rIt = mIt->second.begin( ); + for( ; rIt != mIt->second.end( ); ++rIt ) + { + auto eIt = rIt->second.begin( ); + for( ; eIt != rIt->second.end( ); ++eIt ) + { + auto path = *eIt; + pathCount++; + unsigned int size = path->GetSize( ); + out2 << size << std::endl; + for( unsigned int i = 0; i < path->GetSize( ); ++i ) + { + auto v = path->GetVertex( i ); + for( unsigned int d = 0; d < TSkeleton::Dimension; ++d ) + out2 << v[ d ] << " "; + + } // rof + out2 << std::endl; + + } // rof + + } // rof + + } // rof + out1 << pathCount << std::endl << out2.str( ); + + // Real write + std::ofstream file_stream( this->m_FileName.c_str( ), std::ofstream::binary ); + if( !file_stream ) + itkExceptionMacro( + << "Unable to write skeleton to \"" + << this->m_FileName + << "\"" + ); + file_stream.write( out1.str( ).c_str( ), out1.str( ).size( ) ); +} + +#endif // __fpa__Base__SkeletonWriter__hxx__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/SkeletonFilter.h b/lib/fpa/Image/SkeletonFilter.h index 9f0542f..6ed2b32 100644 --- a/lib/fpa/Image/SkeletonFilter.h +++ b/lib/fpa/Image/SkeletonFilter.h @@ -32,7 +32,8 @@ namespace fpa typedef itk::SmartPointer< Self > Pointer; typedef itk::SmartPointer< const Self > ConstPointer; - typedef typename Superclass::TVertex TVertex; + typedef typename Superclass::TMST TMST; + typedef typename Superclass::TVertex TVertex; typedef typename Superclass::TOutputValue TOutputValue; typedef itk::SignedMaurerDistanceMapImageFilter< TImage, TScalarImage > TDistanceMap; diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx index 2b1cb65..9cb3019 100644 --- a/lib/fpa/Image/SkeletonFilter.hxx +++ b/lib/fpa/Image/SkeletonFilter.hxx @@ -6,6 +6,7 @@ #ifndef __fpa__Image__SkeletonFilter__hxx__ #define __fpa__Image__SkeletonFilter__hxx__ +#include #include // ------------------------------------------------------------------------- @@ -138,20 +139,18 @@ 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 - */ + typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue; + + this->Superclass::_SetOutputValue( vertex, value ); + double d = double( this->m_DistanceMap->GetOutput( )->GetPixel( vertex ) ); + if( d >= double( 0 ) ) + { + // Update skeleton candidates + d += double( 1e-5 ); + double v = double( value ) / ( d * d ); + this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, vertex ) ); + + } // fi } // ------------------------------------------------------------------------- @@ -159,51 +158,51 @@ 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 - */ + typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue; + + TMarks* marks = this->GetMarks( ); + TMST* mst = this->GetMinimumSpanningTree( ); + typename TImage::SpacingType 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 + typename _TSkeletonQueue::iterator nIt = this->m_SkeletonQueue.begin( ); + _TSkeletonQueueValue 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 } // ------------------------------------------------------------------------- @@ -211,69 +210,65 @@ 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 - */ + typedef typename TSkeleton::TPath _TPath; + typedef itk::Image< unsigned long, TImage::ImageDimension > _TTagsImage; + + TMST* mst = this->GetMinimumSpanningTree( ); + TSkeleton* 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( TVertex it: end_points ) + { + 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( TVertex it: end_points ) + { + TVertex p = mst->GetParent( it ); + TVertex sIdx = it; + 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 } // ------------------------------------------------------------------------- @@ -281,54 +276,52 @@ 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 - */ + typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt; + + static const double _eps = std::sqrt( double( TImage::ImageDimension + 1 ) ); + const TScalarImage* dmap = this->m_DistanceMap->GetOutput( ); + TMarks* marks = this->GetMarks( ); + typename TImage::SpacingType spac = dmap->GetSpacing( ); + typename TImage::RegionType region = dmap->GetRequestedRegion( ); + + typename TImage::PointType cnt; + dmap->TransformIndexToPhysicalPoint( idx, cnt ); + double r = double( dmap->GetPixel( idx ) ) * _eps; + + TVertex i0, i1; + for( unsigned int d = 0; d < TImage::ImageDimension; ++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 < TImage::ImageDimension; ++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__