1 #ifndef __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
2 #define __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
6 #include <itkImageRegionConstIteratorWithIndex.h>
8 // -------------------------------------------------------------------------
9 template< class I, class O >
10 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
11 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
15 dynamic_cast< TLabelImage* >(
16 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
21 // -------------------------------------------------------------------------
22 template< class I, class O >
23 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
24 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
25 GetLabelImage( ) const
28 dynamic_cast< const TLabelImage* >(
29 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
34 // -------------------------------------------------------------------------
35 template< class I, class O >
36 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
37 GraftLabelImage( itk::DataObject* obj )
40 dynamic_cast< TLabelImage* >(
41 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
44 this->GraftNthOutput( this->m_LabelImageIndex, lbl );
47 // -------------------------------------------------------------------------
48 template< class I, class O >
49 fpa::Image::DijkstraWithEndPointDetection< I, O >::
50 DijkstraWithEndPointDetection( )
53 this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( );
54 this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 1 );
55 this->itk::ProcessObject::SetNthOutput(
56 this->m_LabelImageIndex, TLabelImage::New( )
60 // -------------------------------------------------------------------------
61 template< class I, class O >
62 fpa::Image::DijkstraWithEndPointDetection< I, O >::
63 ~DijkstraWithEndPointDetection( )
67 // -------------------------------------------------------------------------
68 template< class I, class O >
69 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
70 _BeforeGenerateData( )
74 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
76 _TNode seed = this->m_Seeds[ s ];
77 _TRegion region = this->_Region( seed.Vertex, max_spac );
78 itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region );
81 _TPixel max_value = iIt.Get( );
82 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
84 if( iIt.Get( ) > max_value )
86 seed.Vertex = iIt.GetIndex( );
87 seed.Parent = seed.Vertex;
88 max_value = iIt.Get( );
93 this->m_Seeds[ s ] = seed;
99 this->Superclass::_BeforeGenerateData( );
100 this->m_Candidates.clear( );
103 // -------------------------------------------------------------------------
104 template< class I, class O >
105 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
106 _AfterGenerateData( )
108 this->Superclass::_AfterGenerateData( );
110 // Finish base algorithm
112 this->m_FullTree.clear( );
113 this->m_ReducedTree.clear( );
115 this->m_EndPoints.clear( );
116 this->m_BifurcationPoints.clear( );
117 if( this->m_Candidates.size( ) == 0 )
119 this->InvokeEvent( TEndEvent( ) );
120 this->InvokeEvent( TStartBacktrackingEvent( ) );
122 // Get some input values
123 const I* input = this->GetInput( );
124 typename I::SpacingType spac = input->GetSpacing( );
125 double max_spac = spac[ 0 ];
126 for( unsigned int d = 1; d < I::ImageDimension; ++d )
128 ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
129 max_spac *= double( 3 );
130 TLabelImage* label = this->GetLabelImage( );
131 label->FillBuffer( 0 );
133 // Prepare an object to hold marks
134 std::set< TVertex, TVertexCompare > tree_marks;
136 typename TMarkImage::Pointer marks = this->GetOutputMarkImage( );
137 marks->FillBuffer( 0 );
140 // Iterate over the candidates, starting from the candidates that
141 // are near thin branches
142 typename _TCandidates::const_reverse_iterator cIt =
143 this->m_Candidates.rbegin( );
144 this->m_NumberOfBranches = 1;
145 std::map< TLabel, std::pair< TVertex, TVertex > > branches;
146 for( ; cIt != this->m_Candidates.rend( ); ++cIt )
148 // If pixel has been already labelled, pass
149 TVertex v = cIt->second;
150 if( label->GetPixel( v ) != 0 )
153 // Compute nearest start candidate
154 _TRegion region = this->_Region( v, max_spac );
155 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
157 TVertex max_vertex = iIt.GetIndex( );
158 _TPixel max_value = iIt.Get( );
159 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
161 _TPixel value = iIt.Get( );
162 if( value > max_value )
165 max_vertex = iIt.GetIndex( );
171 // Re-check labelling
172 if( label->GetPixel( max_vertex ) != 0 )
174 this->m_EndPoints.push_back( max_vertex );
176 // Get the path all the way to seed
177 std::vector< TVertex > path;
178 this->GetMinimumSpanningTree( )->
179 GetPath( path, max_vertex, this->GetSeed( 0 ) );
184 TVertex last_start = max_vertex;
185 typename std::vector< TVertex >::const_iterator pIt = path.begin( );
186 for( ; pIt != path.end( ); ++pIt )
189 TBacktrackingEvent( *pIt, this->m_NumberOfBranches )
193 if( tree_marks.find( *pIt ) == tree_marks.end( ) )
195 tree_marks.insert( *pIt );
197 // Mark a sphere around current point as visited
198 double dist = std::sqrt( double( input->GetPixel( *pIt ) ) );
199 region = this->_Region( max_vertex, dist * double( 1.5 ) );
200 itk::ImageRegionIteratorWithIndex< TLabelImage >
201 lIt( label, region );
202 for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
203 lIt.Set( this->m_NumberOfBranches );
205 // Next vertex in current path
206 // TODO: this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
208 this->m_FullTree[ max_vertex ] =
209 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
214 // A bifurcation point has been reached!
215 branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt );
217 this->m_BifurcationPoints.push_back( *pIt );
218 this->m_NumberOfBranches++;
227 this->m_BifurcationPoints.begin( ),
228 this->m_BifurcationPoints.end( ),
230 ) == this->m_BifurcationPoints.end( )
233 // Mark a sphere around current point as visited
234 double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
235 region = this->_Region( max_vertex, dist * double( 1.5 ) );
236 itk::ImageRegionIteratorWithIndex< TLabelImage >
237 lIt( label, region );
238 for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
239 lIt.Set( this->m_NumberOfBranches );
241 // Next vertex in current path
243 this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
244 this->m_FullTree[ max_vertex ] =
245 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
251 // A bifurcation point has been reached!
252 // TODO: this->m_BifurcationPoints.push_back( max_vertex );
253 branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt );
256 this->m_NumberOfBranches++;
258 this->m_FullTree[ max_vertex ] =
259 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
267 if( start || change )
268 this->m_NumberOfBranches++;
269 this->InvokeEvent( TEndBacktrackingEvent( ) );
274 this->InvokeEvent( TEndBacktrackingEvent( ) );
279 // Re-enumerate labels
280 std::map< TLabel, unsigned long > histo;
282 typename std::set< TVertex, TVertexCompare >::iterator treeIt = tree_marks.begin( );
283 treeIt != tree_marks.end( );
286 histo[ label->GetPixel( *treeIt ) ]++;
289 typename std::map< TLabel, unsigned long >::iterator hIt = histo.begin( );
293 std::cout << hIt->first << " " << hIt->second << std::endl;
296 std::map< TMark, TMark > changes;
297 TMark last_change = 1;
298 for( TMark i = 1; i <= this->m_NumberOfBranches; ++i )
300 if( histo[ i ] != 0 )
301 changes[ i ] = last_change++;
304 this->m_NumberOfBranches = changes.size( );
309 typename TTree::iterator treeIt = this->m_FullTree.begin( );
310 treeIt != this->m_FullTree.end( );
314 TMark old = treeIt->second.second;
315 treeIt->second.second = changes[ old ];
319 // Construct reduced tree
321 typename TVertices::const_iterator eIt = this->m_EndPoints.begin( );
322 eIt != this->m_EndPoints.end( );
326 typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt );
327 if( tIt != this->m_FullTree.end( ) )
329 TMark id = tIt->second.second;
332 tIt = this->m_FullTree.find( tIt->second.first );
333 if( tIt == this->m_FullTree.end( ) )
336 } while( tIt->second.second == id );
337 this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id );
344 typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( );
345 bIt != this->m_BifurcationPoints.end( );
349 typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt );
350 if( tIt != this->m_FullTree.end( ) )
352 TMark id = tIt->second.second;
355 tIt = this->m_FullTree.find( tIt->second.first );
356 if( tIt == this->m_FullTree.end( ) )
359 } while( tIt->second.second == id );
360 this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id );
368 // -------------------------------------------------------------------------
369 template< class I, class O >
370 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
371 _SetResult( const TVertex& v, const _TNode& n )
373 this->Superclass::_SetResult( v, n );
375 TResult vv = TResult( this->_VertexValue( v ) );
376 if( TResult( 0 ) < vv )
377 this->m_Candidates.insert( _TCandidate( n.Result / vv, v ) );
380 // -------------------------------------------------------------------------
381 template< class I, class O >
382 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
383 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
384 _Region( const TVertex& c, const double& r )
386 typename I::ConstPointer input = this->GetInput( );
387 typename I::SpacingType spac = input->GetSpacing( );
388 _TRegion region = input->GetLargestPossibleRegion( );
389 typename I::IndexType idx0 = region.GetIndex( );
390 typename I::IndexType idx1 = idx0 + region.GetSize( );
392 // Compute region size and index
393 typename I::IndexType i0, i1;
395 for( unsigned int d = 0; d < I::ImageDimension; ++d )
397 long s = long( std::ceil( r / double( spac[ d ] ) ) );
398 i0[ d ] = c[ d ] - s;
399 i1[ d ] = c[ d ] + s;
401 if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
402 if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
403 if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
404 if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
405 size[ d ] = i1[ d ] - i0[ d ];
409 // Prepare region and return it
410 region.SetIndex( i0 );
411 region.SetSize( size );
415 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__