1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
6 #ifndef __fpa__Image__SkeletonFilter__hxx__
7 #define __fpa__Image__SkeletonFilter__hxx__
9 #include <itkImageRegionIteratorWithIndex.h>
10 #include <itkMinimumMaximumImageCalculator.h>
12 #include <fpa/Image/Functors/VertexIdentity.h>
13 #include <fpa/Base/Functors/InvertValue.h>
15 // -------------------------------------------------------------------------
16 template< class _TImage, class _TScalar >
17 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
18 SetInput( const TScalarImage* image )
23 // -------------------------------------------------------------------------
24 template< class _TImage, class _TScalar >
25 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
26 SetInput( const TImage* image )
28 this->m_DistanceMap->SetInput( image );
29 this->Superclass::SetInput( this->m_DistanceMap->GetOutput( ) );
32 // -------------------------------------------------------------------------
33 template< class _TImage, class _TScalar >
34 _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
37 return( this->m_DistanceMap->GetInput( ) );
40 // -------------------------------------------------------------------------
41 template< class _TImage, class _TScalar >
42 const _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
45 return( this->m_DistanceMap->GetInput( ) );
48 // -------------------------------------------------------------------------
49 template< class _TImage, class _TScalar >
50 typename fpa::Image::SkeletonFilter< _TImage, _TScalar >::
51 TSkeleton* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
55 dynamic_cast< TSkeleton* >(
56 this->itk::ProcessObject::GetOutput( this->m_SkeletonIdx )
61 // -------------------------------------------------------------------------
62 template< class _TImage, class _TScalar >
63 typename fpa::Image::SkeletonFilter< _TImage, _TScalar >::
64 TMarks* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
68 dynamic_cast< TMarks* >(
69 this->itk::ProcessObject::GetOutput( this->m_MarksIdx )
74 // -------------------------------------------------------------------------
75 template< class _TImage, class _TScalar >
76 fpa::Image::SkeletonFilter< _TImage, _TScalar >::
79 m_SeedFromMaximumDistance( false )
81 typedef fpa::Image::Functors::VertexIdentity< TScalarImage, _TScalar > _TVertexFunc;
82 typedef fpa::Base::Functors::InvertValue< _TScalar, _TScalar > _TValueFunc;
84 typename _TVertexFunc::Pointer vertex_func = _TVertexFunc::New( );
85 typename _TValueFunc::Pointer value_func = _TValueFunc::New( );
86 value_func->SetAlpha( 1 );
87 value_func->SetBeta( 1 );
89 this->SetFunctor( vertex_func );
90 this->SetFunctor( value_func );
92 this->m_DistanceMap = TDistanceMap::New( );
93 this->m_DistanceMap->InsideIsPositiveOn( );
94 this->m_DistanceMap->SquaredDistanceOff( );
95 this->m_DistanceMap->UseImageSpacingOn( );
97 unsigned int nOutputs = this->GetNumberOfRequiredOutputs( );
98 this->SetNumberOfRequiredOutputs( nOutputs + 2 );
99 this->m_SkeletonIdx = nOutputs;
100 this->m_MarksIdx = nOutputs + 1;
102 typename TSkeleton::Pointer skeleton = TSkeleton::New( );
103 typename TMarks::Pointer marks = TMarks::New( );
104 this->SetNthOutput( this->m_SkeletonIdx, skeleton.GetPointer( ) );
105 this->SetNthOutput( this->m_MarksIdx, marks.GetPointer( ) );
108 // -------------------------------------------------------------------------
109 template< class _TImage, class _TScalar >
110 fpa::Image::SkeletonFilter< _TImage, _TScalar >::
115 // -------------------------------------------------------------------------
116 template< class _TImage, class _TScalar >
117 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
120 // Prepare distance map
121 this->m_DistanceMap->Update( );
125 if( this->m_SeedFromMaximumDistance )
127 typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
128 typename _TMinMax::Pointer minmax = _TMinMax::New( );
129 minmax->SetImage( this->m_DistanceMap->GetOutput( ) );
131 seed = minmax->GetIndexOfMaximum( );
134 seed = *( this->BeginSeeds( ) );
135 this->m_Seeds.clear( );
136 this->m_Seeds.insert( seed );
138 // Prepare skeleton candidates queue
139 this->m_SkeletonQueue.clear( );
142 this->Superclass::GenerateData( );
146 std::vector< TVertex > end_points;
147 this->_EndPoints( end_points, A );
148 this->_Skeleton( end_points, A );
151 // -------------------------------------------------------------------------
152 template< class _TImage, class _TScalar >
153 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
154 _SetOutputValue( const TVertex& vertex, const TOutputValue& value )
156 typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
158 this->Superclass::_SetOutputValue( vertex, value );
159 double d = double( this->m_DistanceMap->GetOutput( )->GetPixel( vertex ) );
160 if( d >= double( 0 ) )
162 // Update skeleton candidates
164 double v = double( value ) / ( d * d );
165 this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, vertex ) );
170 // -------------------------------------------------------------------------
171 template< class _TImage, class _TScalar >
172 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
173 _EndPoints( std::vector< TVertex >& end_points, _TAdjacencies& A )
175 typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
177 TMarks* marks = this->GetMarks( );
178 TMST* mst = this->GetMinimumSpanningTree( );
179 typename TImage::SpacingType spac = marks->GetSpacing( );
182 marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
183 marks->SetRequestedRegion( mst->GetRequestedRegion( ) );
184 marks->SetBufferedRegion( mst->GetBufferedRegion( ) );
185 marks->SetSpacing( mst->GetSpacing( ) );
186 marks->SetOrigin( mst->GetOrigin( ) );
187 marks->SetDirection( mst->GetDirection( ) );
189 marks->FillBuffer( 0 );
191 // BFS from maximum queue
192 while( this->m_SkeletonQueue.size( ) > 0 )
195 typename _TSkeletonQueue::iterator nIt = this->m_SkeletonQueue.begin( );
196 _TSkeletonQueueValue n = *nIt;
197 this->m_SkeletonQueue.erase( nIt );
199 // Mark it and update end-points
200 unsigned char m = marks->GetPixel( n.second );
203 marks->SetPixel( n.second, 1 );
204 end_points.push_back( n.second );
207 TVertex it = n.second;
208 TVertex p = mst->GetParent( it );
211 this->_MarkSphere( it );
213 p = mst->GetParent( it );
216 this->_MarkSphere( it );
222 // -------------------------------------------------------------------------
223 template< class _TImage, class _TScalar >
224 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
225 _Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A )
227 typedef typename TSkeleton::TPath _TPath;
228 typedef itk::Image< unsigned long, TImage::ImageDimension > _TTagsImage;
230 TMST* mst = this->GetMinimumSpanningTree( );
231 TSkeleton* skeleton = this->GetSkeleton( );
234 typename _TTagsImage::Pointer tags = _TTagsImage::New( );
235 tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
236 tags->SetRequestedRegion( mst->GetRequestedRegion( ) );
237 tags->SetBufferedRegion( mst->GetBufferedRegion( ) );
239 tags->FillBuffer( 0 );
240 typename std::vector< TVertex >::const_iterator eIt = end_points.begin( );
241 for( ; eIt != end_points.end( ); ++eIt )
244 TVertex p = mst->GetParent( it );
247 tags->SetPixel( it, tags->GetPixel( it ) + 1 );
249 p = mst->GetParent( it );
252 tags->SetPixel( it, tags->GetPixel( it ) + 1 );
256 // Build paths (branches)
257 eIt = end_points.begin( );
258 for( ; eIt != end_points.end( ); ++eIt )
261 TVertex p = mst->GetParent( it );
263 typename _TPath::Pointer path = _TPath::New( );
264 path->SetReferenceImage( mst );
267 if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) )
269 // Ok, a new branch has been added
270 path->AddVertex( it );
271 skeleton->AddBranch( path );
273 // Update a new starting index
274 path = _TPath::New( );
275 path->SetReferenceImage( mst );
279 path->AddVertex( it );
281 p = mst->GetParent( it );
285 // Finally, add last branch
286 path->AddVertex( it );
287 skeleton->AddBranch( path );
292 // -------------------------------------------------------------------------
293 template< class _TImage, class _TScalar >
294 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
295 _MarkSphere( const TVertex& idx )
297 typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
299 static const double _eps = std::sqrt( double( TImage::ImageDimension + 1 ) );
300 const TScalarImage* dmap = this->m_DistanceMap->GetOutput( );
301 TMarks* marks = this->GetMarks( );
302 typename TImage::SpacingType spac = dmap->GetSpacing( );
303 typename TImage::RegionType region = dmap->GetRequestedRegion( );
305 typename TImage::PointType cnt;
306 dmap->TransformIndexToPhysicalPoint( idx, cnt );
307 double r = double( dmap->GetPixel( idx ) ) * _eps;
310 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
312 long off = long( std::ceil( r / double( spac[ d ] ) ) );
315 i0[ d ] = idx[ d ] - off;
316 i1[ d ] = idx[ d ] + off;
318 if( i0[ d ] < region.GetIndex( )[ d ] )
319 i0[ d ] = region.GetIndex( )[ d ];
321 if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
322 i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
326 typename TImage::SizeType size;
327 for( unsigned int d = 0; d < TImage::ImageDimension; ++d )
328 size[ d ] = i1[ d ] - i0[ d ] + 1;
330 typename TImage::RegionType neighRegion;
331 neighRegion.SetIndex( i0 );
332 neighRegion.SetSize( size );
334 _TMarksIt mIt( marks, neighRegion );
335 for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
337 TVertex w = mIt.GetIndex( );
338 typename TImage::PointType p;
339 marks->TransformIndexToPhysicalPoint( w, p );
345 #endif // __fpa__Image__SkeletonFilter__hxx__