#include <itkImage.h>
#include <itkImageFileReader.h>
+#include <fpa/Base/SkeletonWriter.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;
-*/
+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[] )
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 );
}
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__SkeletonReader__h__
+#define __fpa__Base__SkeletonReader__h__
+
+#include <itkProcessObject.h>
+
+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 <fpa/Base/SkeletonReader.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __fpa__Base__SkeletonReader__h__
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @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$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__SkeletonWriter__h__
+#define __fpa__Base__SkeletonWriter__h__
+
+#include <itkProcessObject.h>
+
+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 <fpa/Base/SkeletonWriter.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __fpa__Base__SkeletonWriter__h__
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Base__SkeletonWriter__hxx__
+#define __fpa__Base__SkeletonWriter__hxx__
+
+#include <fstream>
+
+// -------------------------------------------------------------------------
+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$
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;
#ifndef __fpa__Image__SkeletonFilter__hxx__
#define __fpa__Image__SkeletonFilter__hxx__
+#include <itkImageRegionIteratorWithIndex.h>
#include <itkMinimumMaximumImageCalculator.h>
// -------------------------------------------------------------------------
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
}
// -------------------------------------------------------------------------
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
}
// -------------------------------------------------------------------------
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
}
// -------------------------------------------------------------------------
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__
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff