From b9128432da0c510bb4ac9a165cc84c065b5e34e8 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Sat, 15 Oct 2016 15:31:01 -0500 Subject: [PATCH] Skeleton filter added. --- lib/fpa/Image/SkeletonFilter.h | 109 ++++++++++++++ lib/fpa/Image/SkeletonFilter.hxx | 246 +++++++++++++++++++++++++++++++ 2 files changed, 355 insertions(+) create mode 100644 lib/fpa/Image/SkeletonFilter.h create mode 100644 lib/fpa/Image/SkeletonFilter.hxx diff --git a/lib/fpa/Image/SkeletonFilter.h b/lib/fpa/Image/SkeletonFilter.h new file mode 100644 index 0000000..4d36fd1 --- /dev/null +++ b/lib/fpa/Image/SkeletonFilter.h @@ -0,0 +1,109 @@ +#ifndef __fpa__Image__SkeletonFilter__h__ +#define __fpa__Image__SkeletonFilter__h__ + +#include +#include +#include +#include + +namespace fpa +{ + namespace Image + { + /** + */ + template< class _TDistanceMap, class _TCostMap > + class SkeletonFilter + : public itk::ProcessObject + { + public: + typedef SkeletonFilter Self; + typedef itk::ProcessObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + +#ifdef ITK_USE_CONCEPT_CHECKING + itkConceptMacro( + SameDimension, + ( itk::Concept::SameDimension< _TDistanceMap::ImageDimension, _TCostMap::ImageDimension > ) + ); +#endif + typedef _TDistanceMap TDistanceMap; + typedef _TCostMap TCostMap; + typedef typename TCostMap::IndexType TIndex; + typedef MinimumSpanningTree< TCostMap::ImageDimension > TMST; + typedef + cpExtensions::DataStructures::Skeleton< TCostMap::ImageDimension > + TSkeleton; + + typedef + itk::Functor::IndexLexicographicCompare< _TCostMap::ImageDimension > + TIndexCompare; + typedef std::set< TIndex, TIndexCompare > TIndicesData; + typedef itk::SimpleDataObjectDecorator< TIndicesData > TIndices; + + public: + itkNewMacro( Self ); + itkTypeMacro( fpa::Image::SkeletonFilter, itk::Object ); + + public: + _TDistanceMap* GetDistanceMap( ); + const _TDistanceMap* GetDistanceMap( ) const; + void SetDistanceMap( _TDistanceMap* dmap ); + + _TCostMap* GetCostMap( ); + const _TCostMap* GetCostMap( ) const; + void SetCostMap( _TCostMap* cmap ); + + TMST* GetMinimumSpanningTree( ); + const TMST* GetMinimumSpanningTree( ) const; + void SetMinimumSpanningTree( TMST* mst ); + + TIndices* GetEndPoints( ); + const TIndices* GetEndPoints( ) const; + + TIndices* GetBifurcations( ); + const TIndices* GetBifurcations( ) const; + + TSkeleton* GetSkeleton( ); + const TSkeleton* GetSkeleton( ) const; + + protected: + SkeletonFilter( ); + virtual ~SkeletonFilter( ); + + virtual void GenerateData( ) fpa_OVERRIDE; + + void _EndPoints( + const TDistanceMap* dmap, + const TCostMap* cmap, + const TMST* mst, + TIndices* end_points + ); + + void _Skeleton( + const TDistanceMap* dmap, + const TCostMap* cmap, + const TMST* mst, + TIndices* end_points + TIndices* bifurcations, + TSkeleton* skeleton + ); + + private: + // Purposely not defined + SkeletonFilter( const Self& other ); + Self& operator=( const Self& other ); + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __fpa__Image__SkeletonFilter__h__ + +// eof - $RCSfile$ diff --git a/lib/fpa/Image/SkeletonFilter.hxx b/lib/fpa/Image/SkeletonFilter.hxx new file mode 100644 index 0000000..65a3fe2 --- /dev/null +++ b/lib/fpa/Image/SkeletonFilter.hxx @@ -0,0 +1,246 @@ +#ifndef __fpa__Image__SkeletonFilter__hxx__ +#define __fpa__Image__SkeletonFilter__hxx__ + +// ------------------------------------------------------------------------- +#define __fpa__Image__SkeletonFilter__InputMacro( input_name, input_type, input_id ) \ + template< class _TDistanceMap, class _TCostMap > \ + typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + input_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Get##input_name( ) \ + { \ + return( \ + dynamic_cast< input_type* >( \ + this->Superclass::GetInput( input_id ) \ + ) \ + ); \ + } \ + template< class _TDistanceMap, class _TCostMap > \ + const \ + typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + input_type* \ + fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Get##input_name( ) const \ + { \ + return( \ + dynamic_cast< const input_type* >( \ + this->Superclass::GetInput( input_id ) \ + ) \ + ); \ + } \ + template< class _TDistanceMap, class _TCostMap > \ + void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Set##input_name( input_type* input ) \ + { \ + this->Superclass::SetInput( input_id, input ); \ + } + +__fpa__Image__SkeletonFilter__InputMacro( DistanceMap, TDistanceMap, 0 ); +__fpa__Image__SkeletonFilter__InputMacro( CostMap, TCostMap, 1 ); +__fpa__Image__SkeletonFilter__InputMacro( MinimumSpanningTree, TMST, 2 ); + +// ------------------------------------------------------------------------- +#define __fpa__Image__SkeletonFilter__OutputMacro( output_name, output_type, output_id ) \ + template< class _TDistanceMap, class _TCostMap > \ + typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + output_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Get##output_name( ) \ + { \ + return( \ + dynamic_cast< output_type* >( \ + this->Superclass::GetOutput( output_id ) \ + ) \ + ); \ + } \ + template< class _TDistanceMap, class _TCostMap > \ + const typename \ + fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + output_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \ + Get##output_name( ) const \ + { \ + return( \ + dynamic_cast< const output_type* >( \ + this->Superclass::GetOutput( output_id ) \ + ) \ + ); \ + } + +__fpa__Image__SkeletonFilter__OutputMacro( Skeleton, TSkeleton, 0 ); +__fpa__Image__SkeletonFilter__OutputMacro( EndPoints, TIndices, 1 ); +__fpa__Image__SkeletonFilter__OutputMacro( Bifurcations, TIndices, 2 ); + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: +SkeletonFilter( ) + : Superclass( ) +{ + this->SetNumbeOfRequiredInputs( 3 ); + this->SetNumbeOfRequiredOutputs( 3 ); + + typename TIndices::Pointer end_points = TIndices::New( ); + typename TIndices::Pointer bifurcations = TIndices::New( ); + typename TSkeleton::Pointer skeleton = TSkeleton::New( ); + this->SetNthOutput( 0, skeleton ); + this->SetNthOutput( 1, end_points ); + this->SetNthOutput( 2, bifurcations ); +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: +~SkeletonFilter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: +GenerateData( ) +{ + // 0. I/O objects + const TDistanceMap* dmap = this->GetDistanceMap( ); + const TCostMap* cmap = this->GetCostMap( ); + const TMST* mst = this->GetMinimumSpanningTree( ); + TIndices* end_points = this->GetEndPoints( ); + TIndices* bifurcations = this->GetBifurcations( ); + TSkeleton* skeleton = this->GetSkeleton( ); + + // 1. Check input correspondance + // TODO + + // 2. Detect end-points + this->_EndPoints( dmap, cmap, mst, end_points ); + + // 3. Build skeleton and keep track of bifurcations + this->_Skeleton( dmap, cmap, mst, end_points, bifurcations, skeleton ); +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: +_EndPoints( + const TDistanceMap* dmap, + const TCostMap* cmap, + const TMST* mst, + TIndices* end_points + ) +{ + typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt; + typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt; + typedef std::multimap< double, TIndex, std::greater< double > > _TQueue; + typedef typename _TQueue::value_type _TQueueValue; + typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks; + + // Some values + typename _TMarks::Pointer marks = _TMarks::New( ); + marks->SetLargestPossibleRegion( dmap->GetLargestPossibleRegion( ) ); + marks->SetRequestedRegion( dmap->GetRequestedRegion( ) ); + marks->SetBufferedRegion( dmap->GetBufferedRegion( ) ); + marks->SetSpacing( dmap->GetSpacing( ) ); + marks->SetOrigin( dmap->GetOrigin( ) ); + marks->SetDirection( dmap->GetDirection( ) ); + marks->Allocate( ); + marks->FillBuffer( 0 ); + + // Create queue + _TQueue queue; + _TDistMapIt dIt( dmap, dmap->GetRequestedRegion( ) ); + _TCostMapIt cIt( cmap, cmap->GetRequestedRegion( ) ); + dIt.GoToBegin( ); + cIt.GoToBegin( ); + for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt ) + { + double d = double( dIt.Get( ) ); + if( d > 0 ) + { + double v = double( cIt.Get( ) ) / ( d * d ); + queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) ); + + } // fi + + } // rof + + // BFS from maximum queue + while( queue.size( ) > 0 ) + { + // Get node + auto nIt = queue.begin( ); + auto n = *nIt; + queue.erase( nIt ); + + unsigned char m = marks->GetPixel( n.second ); + if( ( m & 0x01 ) == 0x01 ) + continue; + + // Mark it and update end-points + m |= 0x01; + marks->SetPixel( n.second, m ); + end_points->Get( ).insert( n.second ); + + // Get path + typename TMST::TPath::Pointer path; + mst->GetPath( path, n.second ); + for( unsigned long i = 0; i < path->GetSize( ); ++i ) + { + TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i ); + typename _TCostMap::PointType cnt; + cmap->TransformIndexToPhysicalPoint( idx, cnt ); + double d = double( dmap->GetPixel( idx ) ); + d *= double( 2 ); + + // Mark sphere + std::queue< TIndex > q; + q.push( idx ); + while( q.size( ) > 0 ) + { + TIndex v = q.front( ); + q.pop( ); + unsigned char m = marks->GetPixel( v ); + if( ( m & 0x02 ) == 0x02 ) + continue; + m |= 0x03; + marks->SetPixel( v, m ); + + for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x ) + { + for( int y = -1; y <= 1; y += 2 ) + { + TIndex w = v; + w[ x ] += y; + if( region.IsInside( w ) ) + { + typename _TCostMap::PointType p; + cmap->TransformIndexToPhysicalPoint( w, p ); + if( cnt.EuclideanDistanceTo( p ) <= d ) + q.push( w ); + + } // fi + + } // rof + + } // rof + + } // elihw + + } // rof + + } // elihw +} + +// ------------------------------------------------------------------------- +template< class _TDistanceMap, class _TCostMap > +void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: +_Skeleton( + const TDistanceMap* dmap, + const TCostMap* cmap, + const TMST* mst, + TIndices* end_points + TIndices* bifurcations, + TSkeleton* skeleton + ) +{ +} + +#endif // __fpa__Image__SkeletonFilter__hxx__ + +// eof - $RCSfile$ -- 2.45.1