1 #ifndef __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
2 #define __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__
6 #include <itkImageRegionConstIteratorWithIndex.h>
7 #include <itkImageRegionIteratorWithIndex.h>
9 // -------------------------------------------------------------------------
10 template< class I, class O >
11 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
12 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
16 dynamic_cast< TLabelImage* >(
17 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
22 // -------------------------------------------------------------------------
23 template< class I, class O >
24 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
25 TLabelImage* fpa::Image::DijkstraWithEndPointDetection< I, O >::
26 GetLabelImage( ) const
29 dynamic_cast< const TLabelImage* >(
30 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
35 // -------------------------------------------------------------------------
36 template< class I, class O >
37 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
38 GraftLabelImage( itk::DataObject* obj )
41 dynamic_cast< TLabelImage* >(
42 this->itk::ProcessObject::GetOutput( this->m_LabelImageIndex )
45 this->GraftNthOutput( this->m_LabelImageIndex, lbl );
48 // -------------------------------------------------------------------------
49 template< class I, class O >
50 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
51 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
55 dynamic_cast< TUniqueVertices* >(
56 this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
61 // -------------------------------------------------------------------------
62 template< class I, class O >
63 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
64 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
68 dynamic_cast< const TUniqueVertices* >(
69 this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
74 // -------------------------------------------------------------------------
75 template< class I, class O >
76 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
77 GraftEndPoints( itk::DataObject* obj )
79 TUniqueVertices* lbl =
80 dynamic_cast< TUniqueVertices* >(
81 this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
84 this->GraftNthOutput( this->m_EndPointsIndex, lbl );
87 // -------------------------------------------------------------------------
88 template< class I, class O >
89 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
90 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
94 dynamic_cast< TUniqueVertices* >(
95 this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
100 // -------------------------------------------------------------------------
101 template< class I, class O >
102 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
103 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
104 GetBifurcations( ) const
107 dynamic_cast< const TUniqueVertices* >(
108 this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
113 // -------------------------------------------------------------------------
114 template< class I, class O >
115 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
116 GraftBifurcations( itk::DataObject* obj )
118 TUniqueVertices* lbl =
119 dynamic_cast< TUniqueVertices* >(
120 this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
123 this->GraftNthOutput( this->m_BifurcationsIndex, lbl );
127 // -------------------------------------------------------------------------
128 template< class I, class O >
129 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
130 TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >::
134 dynamic_cast< TBranches* >(
135 this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
140 // -------------------------------------------------------------------------
141 template< class I, class O >
142 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
143 TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >::
147 dynamic_cast< const TBranches* >(
148 this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
153 // -------------------------------------------------------------------------
154 template< class I, class O >
155 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
156 GraftBranches( itk::DataObject* obj )
159 dynamic_cast< TBranches* >(
160 this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
163 this->GraftNthOutput( this->m_BranchesIndex, lbl );
166 // -------------------------------------------------------------------------
167 template< class I, class O >
168 fpa::Image::DijkstraWithEndPointDetection< I, O >::
169 DijkstraWithEndPointDetection( )
171 m_CorrectSeeds( true ),
172 m_CorrectEndPoints( true ),
173 m_SafetyNeighborhoodSize( 3 )
175 this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( );
176 this->m_EndPointsIndex = this->m_LabelImageIndex + 1;
177 this->m_BifurcationsIndex = this->m_LabelImageIndex + 2;
178 this->m_BranchesIndex = this->m_LabelImageIndex + 3;
179 this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 4 );
180 this->itk::ProcessObject::SetNthOutput(
181 this->m_LabelImageIndex, TLabelImage::New( )
183 this->itk::ProcessObject::SetNthOutput(
184 this->m_EndPointsIndex, TUniqueVertices::New( )
186 this->itk::ProcessObject::SetNthOutput(
187 this->m_BifurcationsIndex, TUniqueVertices::New( )
189 this->itk::ProcessObject::SetNthOutput(
190 this->m_BranchesIndex, TBranches::New( )
194 // -------------------------------------------------------------------------
195 template< class I, class O >
196 fpa::Image::DijkstraWithEndPointDetection< I, O >::
197 ~DijkstraWithEndPointDetection( )
201 // -------------------------------------------------------------------------
202 template< class I, class O >
203 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
204 _BeforeGenerateData( )
206 if( this->m_CorrectSeeds )
208 const I* input = this->GetInput( );
209 typename I::SpacingType spac = input->GetSpacing( );
210 double max_spac = double( spac[ 0 ] );
211 for( unsigned int d = 1; d < I::ImageDimension; ++d )
213 ( double( spac[ d ] ) > max_spac )? double( spac[ d ] ): max_spac;
214 max_spac *= double( 3 );
217 for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s )
219 TVertex seed = this->m_SeedVertices[ s ];
220 _TNode n = this->m_Seeds[ seed ];
221 _TRegion region = this->_Region( seed, max_spac );
222 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
225 _TPixel max_value = iIt.Get( );
226 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
228 if( iIt.Get( ) > max_value )
230 this->m_SeedVertices[ s ] = iIt.GetIndex( );
231 max_value = iIt.Get( );
236 this->m_Seeds.erase( seed );
237 n.Parent = this->m_SeedVertices[ s ];
238 this->m_Seeds[ this->m_SeedVertices[ s ] ] = n;
244 // End initialization
245 this->Superclass::_BeforeGenerateData( );
246 this->m_Candidates.clear( );
249 // -------------------------------------------------------------------------
250 template< class I, class O >
251 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
252 _AfterGenerateData( )
254 // Finish base algorithm
255 this->Superclass::_AfterGenerateData( );
257 // Check if backtracking has any sense
258 if( this->m_Candidates.size( ) == 0 )
260 this->InvokeEvent( TEndEvent( ) );
261 this->InvokeEvent( TStartBacktrackingEvent( ) );
263 // Detect endpoints and bifurcations
264 this->_EndPointsAndBifurcations( );
267 this->_FindBranches( );
269 // Label pixels and branches
273 // -------------------------------------------------------------------------
274 template< class I, class O >
275 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
276 _SetResult( const TVertex& v, const _TNode& n )
278 this->Superclass::_SetResult( v, n );
279 this->m_Candidates.insert( _TCandidate( n.Result, v ) );
282 // -------------------------------------------------------------------------
283 template< class I, class O >
284 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
285 _EndPointsAndBifurcations( )
287 // Prepare backtracking objects
288 TUniqueVertices* endpoints = this->GetEndPoints( );
289 TUniqueVertices* bifurcations = this->GetBifurcations( );
291 bifurcations->Clear( );
293 // Get some input values
294 TVertex seed = this->GetSeed( 0 );
295 const I* input = this->GetInput( );
296 const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
298 // Prepare a pixel marks structure
299 typedef std::map< TVertex, bool, TVertexCompare > _TPixelMarks;
300 _TPixelMarks pixel_marks;
302 // Object to hold tree-nodes marks
303 typedef std::set< TVertex, TVertexCompare > _TTreeMarks;
304 _TTreeMarks tree_marks;
306 // Iterate over the candidates, starting from the candidates that
307 // are near thin branches (high accumulated cost)
308 typename _TCandidates::const_reverse_iterator cIt =
309 this->m_Candidates.rbegin( );
310 for( ; cIt != this->m_Candidates.rend( ); ++cIt )
312 // If pixel has already been visited, pass
313 TVertex v = cIt->second;
314 if( pixel_marks[ v ] )
317 // Correct it to nearest start candidate (high distance value)
318 double vr = std::sqrt( double( input->GetPixel( v ) ) );
319 if( this->m_CorrectEndPoints )
320 v = this->_MaxInRegion( input, v, vr * double( 1.5 ) );
322 // Now, check for real marking conditions
323 // 1. Has it been visited by dijkstra?
324 if( this->_Node( v ).Label == Self::FarLabel )
326 // 2. Is it already marked?
327 if( pixel_marks[ v ] )
329 // 3. Is it already an endpoint?
330 if( endpoints->Find( v ) )
332 // 4. Ok, it is completely new!
333 endpoints->Insert( v );
335 // Get the path all the way to global seed
336 TVertices path = mst->GetPath( v, seed );
338 // Backtracking to find endpoints and bifurcations
339 bool adding_new_points = true;
340 typename TVertices::const_iterator pIt = path.begin( );
341 for( ; pIt != path.end( ) && adding_new_points; ++pIt )
343 this->InvokeEvent( TBacktrackingEvent( *pIt, 1 ) );
344 if( tree_marks.find( *pIt ) == tree_marks.end( ) )
346 // Mark current point as a tree point
347 tree_marks.insert( *pIt );
349 // Mark a region around current point as visited
350 vr = std::sqrt( double( input->GetPixel( *pIt ) ) );
351 _TRegion region = this->_Region( *pIt, vr );
352 if( region.GetNumberOfPixels( ) > 0 )
354 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
355 for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt )
356 pixel_marks[ iIt.GetIndex( ) ] = true;
359 pixel_marks[ *pIt ] = true;
363 // A bifurcation point has been reached!
366 bifurcations->Insert( *pIt );
367 adding_new_points = false;
374 this->InvokeEvent( TEndBacktrackingEvent( ) );
379 // -------------------------------------------------------------------------
380 template< class I, class O >
381 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
384 // Configure and obtain inputs
385 TVertex seed = this->GetSeed( 0 );
386 const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
387 const TUniqueVertices* endpoints = this->GetEndPoints( );
388 const TUniqueVertices* bifurcations = this->GetBifurcations( );
389 TBranches* branches = this->GetBranches( );
392 // Reconstruct pixels
393 typename TUniqueVertices::ConstIterator eIt = endpoints->Begin( );
394 for( ; eIt != endpoints->End( ); ++eIt )
396 // Get the path all the way to global seed
397 TVertices path = mst->GetPath( *eIt, seed );
399 TVertex start_vertex = *eIt;
400 typename TVertices::const_iterator pIt = path.begin( );
401 for( ; pIt != path.end( ); ++pIt )
403 if( bifurcations->Find( *pIt ) )
405 branches->SetValue( start_vertex, *pIt, 1 );
412 // Finish with branches to global seed
413 branches->SetValue( start_vertex, seed, 1 );
418 // -------------------------------------------------------------------------
419 template< class I, class O >
420 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
423 // Configure and obtain inputs
424 const I* input = this->GetInput( );
425 const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
426 const TUniqueVertices* bifurcations = this->GetBifurcations( );
427 TBranches* branches = this->GetBranches( );
428 TLabelImage* label = this->GetLabelImage( );
429 label->FillBuffer( 0 );
432 typename TBranches::Iterator bIt = branches->Begin( );
433 TLabel actual_label = 0;
434 for( ; bIt != branches->End( ); ++bIt )
436 typename TBranches::RowIterator brIt = branches->Begin( bIt );
437 for( ; brIt != branches->End( bIt ); ++brIt )
440 brIt->second = actual_label;
442 TVertices path = mst->GetPath( bIt->first, brIt->first );
443 typename TVertices::const_iterator pIt = path.begin( );
444 for( ; pIt != path.end( ); ++pIt )
446 double d = std::sqrt( double( input->GetPixel( *pIt ) ) );
447 _TRegion region = this->_Region( *pIt, d );
448 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
449 itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
452 for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
454 // Mask real input image
456 !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) &&
457 lIt.Get( ) == TLabel( 0 )
459 lIt.Set( actual_label );
468 this->m_NumberOfBranches = actual_label;
470 // Label bifurcations
471 typename TUniqueVertices::ConstIterator bifIt =
472 bifurcations->Begin( );
473 for( ; bifIt != bifurcations->End( ); ++bifIt )
476 double d = std::sqrt( double( input->GetPixel( *bifIt ) ) );
477 _TRegion region = this->_Region( *bifIt, d * 1.5 );
478 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
479 itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
482 for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
484 // Mask real input image
485 if( !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) )
486 lIt.Set( actual_label );
493 // -------------------------------------------------------------------------
494 template< class I, class O >
495 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
496 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
497 _Region( const TVertex& c, const double& r )
499 typename I::ConstPointer input = this->GetInput( );
500 typename I::SpacingType spac = input->GetSpacing( );
501 _TRegion region = input->GetLargestPossibleRegion( );
502 typename I::IndexType idx0 = region.GetIndex( );
503 typename I::IndexType idx1 = idx0 + region.GetSize( );
505 // Compute region size and index
506 typename I::IndexType i0, i1;
508 for( unsigned int d = 0; d < I::ImageDimension; ++d )
510 // NOTE: 3 is a minimum neighborhood size
512 long( std::ceil( r / double( spac[ d ] ) ) ) +
513 long( this->m_SafetyNeighborhoodSize );
514 i0[ d ] = c[ d ] - s;
515 i1[ d ] = c[ d ] + s;
517 if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
518 if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
519 if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
520 if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
521 size[ d ] = i1[ d ] - i0[ d ];
525 // Prepare region and return it
526 region.SetIndex( i0 );
527 region.SetSize( size );
532 // -------------------------------------------------------------------------
533 template< class I, class O >
535 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
536 TVertex fpa::Image::DijkstraWithEndPointDetection< I, O >::
537 _MaxInRegion( const _T* image, const TVertex& v, const double& r )
539 typedef itk::ImageRegionConstIteratorWithIndex< _T > _TIt;
541 _TIt iIt( image, this->_Region( v, r ) );
543 TVertex max_vertex = iIt.GetIndex( );
544 typename _T::PixelType max_value = iIt.Get( );
545 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
547 typename _T::PixelType value = iIt.Get( );
548 if( value > max_value )
551 max_vertex = iIt.GetIndex( );
556 return( max_vertex );
559 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__