1 #ifndef __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
2 #define __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
5 #include <itkImageRegionConstIteratorWithIndex.h>
7 // -------------------------------------------------------------------------
8 template< class I, class O >
9 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
10 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
14 dynamic_cast< TLabelImage* >(
15 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
20 // -------------------------------------------------------------------------
21 template< class I, class O >
22 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
23 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
24 GetLabelImage( ) const
27 dynamic_cast< const TLabelImage* >(
28 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
33 // -------------------------------------------------------------------------
34 template< class I, class O >
35 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
36 GraftLabelImage( itk::DataObject* obj )
39 dynamic_cast< TLabelImage* >(
40 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
43 this->GraftNthOutput( this->m_LabelImageIndex, lbl );
46 // -------------------------------------------------------------------------
47 template< class I, class O >
48 fpa::Image::DijkstraWithEndPointDetection< I, O >::
49 DijkstraWithEndPointDetection( )
52 this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( );
53 this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 1 );
54 this->itk::ProcessObject::SetNthOutput(
55 this->m_LabelImageIndex, TLabelImage::New( )
59 // -------------------------------------------------------------------------
60 template< class I, class O >
61 fpa::Image::DijkstraWithEndPointDetection< I, O >::
62 ~DijkstraWithEndPointDetection( )
66 // -------------------------------------------------------------------------
67 template< class I, class O >
68 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
69 _BeforeGenerateData( )
73 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
75 _TNode seed = this->m_Seeds[ s ];
76 _TRegion region = this->_Region( seed.Vertex, max_spac );
77 itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region );
80 _TPixel max_value = iIt.Get( );
81 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
83 if( iIt.Get( ) > max_value )
85 seed.Vertex = iIt.GetIndex( );
86 seed.Parent = seed.Vertex;
87 max_value = iIt.Get( );
92 this->m_Seeds[ s ] = seed;
98 this->Superclass::_BeforeGenerateData( );
99 this->m_Candidates.clear( );
102 // -------------------------------------------------------------------------
103 template< class I, class O >
104 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
105 _AfterGenerateData( )
107 this->Superclass::_AfterGenerateData( );
109 // Finish base algorithm
111 this->m_FullTree.clear( );
112 this->m_ReducedTree.clear( );
114 this->m_EndPoints.clear( );
115 this->m_BifurcationPoints.clear( );
116 if( this->m_Candidates.size( ) == 0 )
118 this->InvokeEvent( TEndEvent( ) );
119 this->InvokeEvent( TStartBacktrackingEvent( ) );
121 // Get some input values
122 const I* input = this->GetInput( );
123 typename I::SpacingType spac = input->GetSpacing( );
124 double max_spac = spac[ 0 ];
125 for( unsigned int d = 1; d < I::ImageDimension; ++d )
127 ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
128 max_spac *= double( 3 );
129 TLabelImage* label = this->GetLabelImage( );
130 label->FillBuffer( 0 );
132 // Prepare an object to hold marks
134 typename TMarkImage::Pointer marks = this->GetOutputMarkImage( );
135 marks->FillBuffer( 0 );
138 // Iterate over the candidates, starting from the candidates that
139 // are near thin branches
140 typename _TCandidates::const_reverse_iterator cIt =
141 this->m_Candidates.rbegin( );
142 this->m_NumberOfBranches = 1;
143 std::map< TLabel, std::pair< TVertex, TVertex > > branches;
144 for( ; cIt != this->m_Candidates.rend( ); ++cIt )
146 // If pixel hasn't been visited, continue
147 TVertex v = cIt->second;
148 if( label->GetPixel( v ) != 0 )
151 // Compute nearest start candidate
152 _TRegion region = this->_Region( v, max_spac );
153 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
155 TVertex max_vertex = iIt.GetIndex( );
156 _TPixel max_value = iIt.Get( );
157 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
159 _TPixel value = iIt.Get( );
160 if( value > max_value )
163 max_vertex = iIt.GetIndex( );
170 if( label->GetPixel( max_vertex ) != 0 )
172 this->m_EndPoints.push_back( max_vertex );
174 // Get the path all the way to seed
175 std::vector< TVertex > path;
176 this->GetMinimumSpanningTree( )->
177 GetPath( path, max_vertex, this->GetSeed( 0 ) );
182 TVertex last_start = max_vertex;
183 typename std::vector< TVertex >::const_iterator pIt = path.begin( );
184 for( ; pIt != path.end( ); ++pIt )
187 TBacktrackingEvent( *pIt, this->m_NumberOfBranches )
191 TLabel lbl = label->GetPixel( *pIt );
192 if( lbl == 0 || lbl == this->m_NumberOfBranches )
194 // Mark a sphere around current point as visited
195 double dist = std::sqrt( double( input->GetPixel( *pIt ) ) );
196 region = this->_Region( max_vertex, dist * double( 1.1 ) );
197 itk::ImageRegionIteratorWithIndex< TLabelImage >
198 lIt( label, region );
199 for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
201 if( lIt.Get( ) == 0 )
202 lIt.Set( this->m_NumberOfBranches );
206 // Next vertex in current path
207 // TODO: this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
209 this->m_FullTree[ max_vertex ] =
210 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
211 std::cout << "New: " << this->m_NumberOfBranches << std::endl;
216 // A bifurcation point has been reached!
217 branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt );
218 std::cout << "bif " << this->m_NumberOfBranches << std::endl;
220 this->m_BifurcationPoints.push_back( *pIt );
221 this->m_NumberOfBranches++;
230 this->m_BifurcationPoints.begin( ),
231 this->m_BifurcationPoints.end( ),
233 ) == this->m_BifurcationPoints.end( )
236 // Mark a sphere around current point as visited
237 double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
238 region = this->_Region( max_vertex, dist * double( 1.1 ) );
239 itk::ImageRegionIteratorWithIndex< TLabelImage >
240 lIt( label, region );
241 for( lIt.GoToBegin( ); !lIt.IsAtEnd( ); ++lIt )
242 lIt.Set( this->m_NumberOfBranches );
244 // Next vertex in current path
246 this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
247 this->m_FullTree[ max_vertex ] =
248 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
251 // std::cout << "Change: " << this->m_NumberOfBranches << std::endl;
255 // A bifurcation point has been reached!
256 // TODO: this->m_BifurcationPoints.push_back( max_vertex );
257 branches[ this->m_NumberOfBranches ] = std::pair< TVertex, TVertex >( last_start, *pIt );
258 std::cout << "change " << this->m_NumberOfBranches << std::endl;
261 this->m_NumberOfBranches++;
263 this->m_FullTree[ max_vertex ] =
264 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
266 // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl;
273 if( start || change )
274 this->m_NumberOfBranches++;
275 this->InvokeEvent( TEndBacktrackingEvent( ) );
280 this->InvokeEvent( TEndBacktrackingEvent( ) );
285 std::cout << this->m_NumberOfBranches << " " << branches.size( ) << std::endl;
286 std::cout << this->m_EndPoints.size( ) << " " << this->m_BifurcationPoints.size( ) << std::endl;
289 // Re-enumerate labels
291 std::map< TLabel, unsigned long > histo;
293 typename TTree::iterator treeIt = this->m_FullTree.begin( );
294 treeIt != this->m_FullTree.end( );
297 histo[ treeIt->second.second ]++;
299 std::map< TMark, TMark > changes;
300 TMark last_change = 1;
301 for( TMark i = 1; i <= this->m_NumberOfBranches; ++i )
303 if( histo[ i ] != 0 )
304 changes[ i ] = last_change++;
307 this->m_NumberOfBranches = changes.size( );
312 typename TTree::iterator treeIt = this->m_FullTree.begin( );
313 treeIt != this->m_FullTree.end( );
317 TMark old = treeIt->second.second;
318 treeIt->second.second = changes[ old ];
322 // Construct reduced tree
324 typename TVertices::const_iterator eIt = this->m_EndPoints.begin( );
325 eIt != this->m_EndPoints.end( );
329 typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt );
330 if( tIt != this->m_FullTree.end( ) )
332 TMark id = tIt->second.second;
335 tIt = this->m_FullTree.find( tIt->second.first );
336 if( tIt == this->m_FullTree.end( ) )
339 } while( tIt->second.second == id );
340 this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id );
347 typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( );
348 bIt != this->m_BifurcationPoints.end( );
352 typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt );
353 if( tIt != this->m_FullTree.end( ) )
355 TMark id = tIt->second.second;
358 tIt = this->m_FullTree.find( tIt->second.first );
359 if( tIt == this->m_FullTree.end( ) )
362 } while( tIt->second.second == id );
363 this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id );
371 // -------------------------------------------------------------------------
372 template< class I, class O >
373 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
374 _SetResult( const TVertex& v, const _TNode& n )
376 this->Superclass::_SetResult( v, n );
378 TResult vv = TResult( this->_VertexValue( v ) );
379 if( TResult( 0 ) < vv )
380 this->m_Candidates.insert( _TCandidate( n.Result / vv, v ) );
383 // -------------------------------------------------------------------------
384 template< class I, class O >
385 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
386 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
387 _Region( const TVertex& c, const double& r )
389 typename I::ConstPointer input = this->GetInput( );
390 typename I::SpacingType spac = input->GetSpacing( );
391 _TRegion region = input->GetLargestPossibleRegion( );
392 typename I::IndexType idx0 = region.GetIndex( );
393 typename I::IndexType idx1 = idx0 + region.GetSize( );
395 // Compute region size and index
396 typename I::IndexType i0, i1;
398 for( unsigned int d = 0; d < I::ImageDimension; ++d )
400 long s = long( std::ceil( r / double( spac[ d ] ) ) );
401 i0[ d ] = c[ d ] - s;
402 i1[ d ] = c[ d ] + s;
404 if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
405 if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
406 if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
407 if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
408 size[ d ] = i1[ d ] - i0[ d ];
412 // Prepare region and return it
413 region.SetIndex( i0 );
414 region.SetSize( size );
418 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__