1 #ifndef __fpa__Image__SkeletonFilter__hxx__
2 #define __fpa__Image__SkeletonFilter__hxx__
6 #include <itkImageRegionConstIteratorWithIndex.h>
8 #include <fpa/Base/Functors/Inverse.h>
9 #include <fpa/Image/Functors/SimpleNeighborhood.h>
10 #include <itkImageRegionIteratorWithIndex.h>
12 // -------------------------------------------------------------------------
13 #define fpa_Image_SkeletonFilter_OutputMacro( o_n, o_t ) \
14 template< class _TImage > \
15 typename fpa::Image::SkeletonFilter< _TImage >:: \
16 o_t* fpa::Image::SkeletonFilter< _TImage >:: \
20 dynamic_cast< o_t* >( \
21 this->itk::ProcessObject::GetOutput( this->m_##o_n##Idx ) \
26 fpa_Image_SkeletonFilter_OutputMacro( Skeleton, TSkeleton );
27 fpa_Image_SkeletonFilter_OutputMacro( EndPoints, TIndices );
28 fpa_Image_SkeletonFilter_OutputMacro( Bifurcations, TIndices );
29 fpa_Image_SkeletonFilter_OutputMacro( Marks, TMarks );
31 // -------------------------------------------------------------------------
32 template< class _TImage >
33 fpa::Image::SkeletonFilter< _TImage >::
37 typedef typename _TImage::PixelType _TPixel;
38 typedef typename _TImage::Superclass _TImageBase;
39 typedef fpa::Image::Functors::SimpleNeighborhood< _TImageBase > _TNeighFunc;
40 typedef fpa::Base::Functors::Inverse< _TPixel, _TPixel > _TInvFunc;
42 unsigned int nOutputs = this->GetNumberOfRequiredOutputs( );
43 this->SetNumberOfRequiredOutputs( nOutputs + 4 );
44 this->m_EndPointsIdx = nOutputs;
45 this->m_BifurcationsIdx = nOutputs + 1;
46 this->m_SkeletonIdx = nOutputs + 2;
47 this->m_MarksIdx = nOutputs + 3;
49 typename TIndices::Pointer end_points = TIndices::New( );
50 typename TIndices::Pointer bifurcations = TIndices::New( );
51 typename TSkeleton::Pointer skeleton = TSkeleton::New( );
52 typename TMarks::Pointer marks = TMarks::New( );
53 this->SetNthOutput( this->m_EndPointsIdx, end_points.GetPointer( ) );
54 this->SetNthOutput( this->m_BifurcationsIdx, bifurcations.GetPointer( ) );
55 this->SetNthOutput( this->m_SkeletonIdx, skeleton.GetPointer( ) );
56 this->SetNthOutput( this->m_MarksIdx, marks.GetPointer( ) );
58 typename _TNeighFunc::Pointer nfunc = _TNeighFunc::New( );
60 this->SetNeighborhoodFunction( nfunc );
62 typename _TInvFunc::Pointer ifunc = _TInvFunc::New( );
63 this->SetCostConversionFunction( ifunc );
66 // -------------------------------------------------------------------------
67 template< class _TImage >
68 fpa::Image::SkeletonFilter< _TImage >::
73 // -------------------------------------------------------------------------
74 template< class _TImage >
75 void fpa::Image::SkeletonFilter< _TImage >::
76 _BeforeGenerateData( )
78 this->Superclass::_BeforeGenerateData( );
79 this->m_SkeletonQueue.clear( );
82 // -------------------------------------------------------------------------
83 template< class _TImage >
84 void fpa::Image::SkeletonFilter< _TImage >::
85 _UpdateResult( const _TQueueNode& n )
87 typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
89 this->Superclass::_UpdateResult( n );
91 // Update skeleton candidates
92 double d = double( this->GetInput( )->GetPixel( n.Vertex ) );
93 if( d >= double( 0 ) )
96 double v = double( n.Result ) / ( d * d );
97 this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) );
102 // -------------------------------------------------------------------------
103 template< class _TImage >
104 void fpa::Image::SkeletonFilter< _TImage >::
105 _AfterGenerateData( )
107 this->Superclass::_AfterGenerateData( );
112 // -------------------------------------------------------------------------
113 template< class _TImage >
114 void fpa::Image::SkeletonFilter< _TImage >::
117 typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
120 typedef itk::ImageRegionConstIteratorWithIndex< _TImage > _TDistMapIt;
121 typedef itk::ImageRegionConstIteratorWithIndex< _TOutput > _TOutputIt;
122 typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
123 typedef typename _TQueue::value_type _TQueueValue;
126 static const double _eps = std::sqrt( double( Self::Dimension + 1 ) );
127 auto input = this->GetInput( );
128 auto mst = this->GetMinimumSpanningTree( );
129 auto marks = this->GetMarks( );
130 auto& end_points = this->GetEndPoints( )->Get( );
133 marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
134 marks->SetRequestedRegion( input->GetRequestedRegion( ) );
135 marks->SetBufferedRegion( input->GetBufferedRegion( ) );
136 marks->SetSpacing( input->GetSpacing( ) );
137 marks->SetOrigin( input->GetOrigin( ) );
138 marks->SetDirection( input->GetDirection( ) );
140 marks->FillBuffer( 0 );
142 // BFS from maximum queue
143 auto region = input->GetRequestedRegion( );
144 while( this->m_SkeletonQueue.size( ) > 0 )
147 auto nIt = this->m_SkeletonQueue.begin( );
149 this->m_SkeletonQueue.erase( nIt );
151 // Mark it and update end-points
152 unsigned char m = marks->GetPixel( n.second );
155 marks->SetPixel( n.second, 1 );
156 end_points.insert( n.second );
159 auto spac = marks->GetSpacing( );
160 typename TMST::TPath::Pointer path;
161 mst->GetPath( path, n.second );
162 for( unsigned long i = 0; i < path->GetSize( ); ++i )
164 TIndex idx = path->GetVertex( i );
165 typename _TImage::PointType cnt;
166 input->TransformIndexToPhysicalPoint( idx, cnt );
167 double r = double( input->GetPixel( idx ) ) * _eps;
170 for( unsigned int d = 0; d < Self::Dimension; ++d )
172 long off = long( std::ceil( r / double( spac[ d ] ) ) );
175 i0[ d ] = idx[ d ] - off;
176 i1[ d ] = idx[ d ] + off;
178 if( i0[ d ] < region.GetIndex( )[ d ] )
179 i0[ d ] = region.GetIndex( )[ d ];
181 if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
182 i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
186 typename _TImage::SizeType size;
187 for( unsigned int d = 0; d < Self::Dimension; ++d )
188 size[ d ] = i1[ d ] - i0[ d ] + 1;
190 typename _TImage::RegionType neighRegion;
191 neighRegion.SetIndex( i0 );
192 neighRegion.SetSize( size );
194 _TMarksIt mIt( marks, neighRegion );
195 for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
197 TIndex w = mIt.GetIndex( );
198 typename _TImage::PointType p;
199 marks->TransformIndexToPhysicalPoint( w, p );
209 // -------------------------------------------------------------------------
210 template< class _TImage >
211 void fpa::Image::SkeletonFilter< _TImage >::
214 auto mst = this->GetMinimumSpanningTree( );
215 auto skeleton = this->GetSkeleton( );
216 auto& end_points = this->GetEndPoints( )->Get( );
218 typedef typename TMST::TPath _TPath;
219 for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
221 typename _TPath::Pointer path;
222 mst->GetPath( path, *eIt );
223 skeleton->AddBranch( path );
228 // -------------------------------------------------------------------------
230 template< class _TImage >
231 void fpa::Image::SkeletonFilter< _TImage >::
233 const TDistanceMap* dmap,
234 const TCostMap* cmap,
236 const TIndicesData& end_points,
237 TIndicesData& bifurcations,
244 #endif // __fpa__Image__SkeletonFilter__hxx__