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 );
51 // -------------------------------------------------------------------------
52 template< class I, class O >
53 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
54 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
58 dynamic_cast< TUniqueVertices* >(
59 this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
64 // -------------------------------------------------------------------------
65 template< class I, class O >
66 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
67 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
71 dynamic_cast< const TUniqueVertices* >(
72 this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
77 // -------------------------------------------------------------------------
78 template< class I, class O >
79 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
80 GraftEndPoints( itk::DataObject* obj )
82 TUniqueVertices* lbl =
83 dynamic_cast< TUniqueVertices* >(
84 this->itk::ProcessObject::GetOutput( this->m_EndPointsIndex )
87 this->GraftNthOutput( this->m_EndPointsIndex, lbl );
90 // -------------------------------------------------------------------------
91 template< class I, class O >
92 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
93 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
97 dynamic_cast< TUniqueVertices* >(
98 this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
103 // -------------------------------------------------------------------------
104 template< class I, class O >
105 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
106 TUniqueVertices* fpa::Image::DijkstraWithEndPointDetection< I, O >::
107 GetBifurcations( ) const
110 dynamic_cast< const TUniqueVertices* >(
111 this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
116 // -------------------------------------------------------------------------
117 template< class I, class O >
118 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
119 GraftBifurcations( itk::DataObject* obj )
121 TUniqueVertices* lbl =
122 dynamic_cast< TUniqueVertices* >(
123 this->itk::ProcessObject::GetOutput( this->m_BifurcationsIndex )
126 this->GraftNthOutput( this->m_BifurcationsIndex, lbl );
130 // -------------------------------------------------------------------------
131 template< class I, class O >
132 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
133 TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >::
137 dynamic_cast< TBranches* >(
138 this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
143 // -------------------------------------------------------------------------
144 template< class I, class O >
145 const typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
146 TBranches* fpa::Image::DijkstraWithEndPointDetection< I, O >::
150 dynamic_cast< const TBranches* >(
151 this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
156 // -------------------------------------------------------------------------
157 template< class I, class O >
158 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
159 GraftBranches( itk::DataObject* obj )
162 dynamic_cast< TBranches* >(
163 this->itk::ProcessObject::GetOutput( this->m_BranchesIndex )
166 this->GraftNthOutput( this->m_BranchesIndex, lbl );
169 // -------------------------------------------------------------------------
170 template< class I, class O >
171 fpa::Image::DijkstraWithEndPointDetection< I, O >::
172 DijkstraWithEndPointDetection( )
174 m_CorrectSeeds( true ),
175 m_CorrectEndPoints( true ),
176 m_SafetyNeighborhoodSize( 3 )
178 this->m_LabelImageIndex = this->GetNumberOfRequiredOutputs( );
179 this->m_EndPointsIndex = this->m_LabelImageIndex + 1;
180 this->m_BifurcationsIndex = this->m_LabelImageIndex + 2;
181 this->m_BranchesIndex = this->m_LabelImageIndex + 3;
182 this->SetNumberOfRequiredOutputs( this->m_LabelImageIndex + 4 );
183 this->itk::ProcessObject::SetNthOutput(
184 this->m_LabelImageIndex, TLabelImage::New( )
186 this->itk::ProcessObject::SetNthOutput(
187 this->m_EndPointsIndex, TUniqueVertices::New( )
189 this->itk::ProcessObject::SetNthOutput(
190 this->m_BifurcationsIndex, TUniqueVertices::New( )
192 this->itk::ProcessObject::SetNthOutput(
193 this->m_BranchesIndex, TBranches::New( )
197 // -------------------------------------------------------------------------
198 template< class I, class O >
199 fpa::Image::DijkstraWithEndPointDetection< I, O >::
200 ~DijkstraWithEndPointDetection( )
204 // -------------------------------------------------------------------------
205 template< class I, class O >
206 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
207 _BeforeGenerateData( )
209 if( this->m_CorrectSeeds )
211 const I* input = this->GetInput( );
212 typename I::SpacingType spac = input->GetSpacing( );
213 double max_spac = double( spac[ 0 ] );
214 for( unsigned int d = 1; d < I::ImageDimension; ++d )
216 ( double( spac[ d ] ) > max_spac )? double( spac[ d ] ): max_spac;
217 max_spac *= double( 3 );
220 for( unsigned int s = 0; s < this->m_SeedVertices.size( ); ++s )
222 TVertex seed = this->m_SeedVertices[ s ];
223 _TNode n = this->m_Seeds[ seed ];
224 _TRegion region = this->_Region( seed, max_spac );
225 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
228 _TPixel max_value = iIt.Get( );
229 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
231 if( iIt.Get( ) > max_value )
233 this->m_SeedVertices[ s ] = iIt.GetIndex( );
234 max_value = iIt.Get( );
239 this->m_Seeds.erase( seed );
240 n.Parent = this->m_SeedVertices[ s ];
241 this->m_Seeds[ this->m_SeedVertices[ s ] ] = n;
247 // End initialization
248 this->Superclass::_BeforeGenerateData( );
249 this->m_Candidates.clear( );
252 // -------------------------------------------------------------------------
253 template< class I, class O >
254 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
255 _AfterGenerateData( )
257 // Finish base algorithm
258 this->Superclass::_AfterGenerateData( );
260 // Check if backtracking has any sense
261 if( this->m_Candidates.size( ) == 0 )
263 this->InvokeEvent( TEndEvent( ) );
264 this->InvokeEvent( TStartBacktrackingEvent( ) );
266 // Detect endpoints and bifurcations
267 this->_EndPointsAndBifurcations( );
270 this->_FindBranches( );
272 // Label pixels and branches
276 // -------------------------------------------------------------------------
277 template< class I, class O >
278 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
279 _SetResult( const TVertex& v, const _TNode& n )
281 this->Superclass::_SetResult( v, n );
282 this->m_Candidates.insert( _TCandidate( n.Result, v ) );
285 // -------------------------------------------------------------------------
286 template< class I, class O >
287 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
288 _EndPointsAndBifurcations( )
290 // Prepare backtracking objects
291 TUniqueVertices* endpoints = this->GetEndPoints( );
292 TUniqueVertices* bifurcations = this->GetBifurcations( );
294 bifurcations->Clear( );
296 // Get some input values
297 TVertex seed = this->GetSeed( 0 );
298 const I* input = this->GetInput( );
299 const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
301 // Prepare a pixel marks structure
302 typedef std::map< TVertex, bool, TVertexCompare > _TPixelMarks;
303 _TPixelMarks pixel_marks;
305 // Object to hold tree-nodes marks
306 typedef std::set< TVertex, TVertexCompare > _TTreeMarks;
307 _TTreeMarks tree_marks;
309 // Iterate over the candidates, starting from the candidates that
310 // are near thin branches (high accumulated cost)
311 typename _TCandidates::const_reverse_iterator cIt =
312 this->m_Candidates.rbegin( );
313 for( ; cIt != this->m_Candidates.rend( ); ++cIt )
315 // If pixel has already been visited, pass
316 TVertex v = cIt->second;
317 if( pixel_marks[ v ] )
320 // Correct it to nearest start candidate (high distance value)
321 double vr = std::sqrt( double( input->GetPixel( v ) ) );
322 if( this->m_CorrectEndPoints )
323 v = this->_MaxInRegion( input, v, vr * double( 1.5 ) );
325 // Now, check for real marking conditions
326 // 1. Has it been visited by dijkstra?
327 if( this->_Node( v ).Label == Self::FarLabel )
329 // 2. Is it already marked?
330 if( pixel_marks[ v ] )
332 // 3. Is it already an endpoint?
333 if( endpoints->Find( v ) )
335 // 4. Ok, it is completely new!
336 endpoints->Insert( v );
338 // Get the path all the way to global seed
340 mst->GetPath( path, v, seed );
342 // Backtracking to find endpoints and bifurcations
343 bool adding_new_points = true;
344 typename TVertices::const_iterator pIt = path.begin( );
345 for( ; pIt != path.end( ) && adding_new_points; ++pIt )
347 this->InvokeEvent( TBacktrackingEvent( *pIt, 1 ) );
348 if( tree_marks.find( *pIt ) == tree_marks.end( ) )
350 // Mark current point as a tree point
351 tree_marks.insert( *pIt );
353 // Mark a region around current point as visited
354 vr = std::sqrt( double( input->GetPixel( *pIt ) ) );
355 _TRegion region = this->_Region( *pIt, vr );
356 if( region.GetNumberOfPixels( ) > 0 )
358 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
359 for( iIt.GoToBegin( ); !iIt.IsAtEnd( ); ++iIt )
360 pixel_marks[ iIt.GetIndex( ) ] = true;
363 pixel_marks[ *pIt ] = true;
367 // A bifurcation point has been reached!
370 bifurcations->Insert( *pIt );
371 adding_new_points = false;
378 this->InvokeEvent( TEndBacktrackingEvent( ) );
383 // -------------------------------------------------------------------------
384 template< class I, class O >
385 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
388 // Configure and obtain inputs
389 TVertex seed = this->GetSeed( 0 );
390 const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
391 const TUniqueVertices* endpoints = this->GetEndPoints( );
392 const TUniqueVertices* bifurcations = this->GetBifurcations( );
393 TBranches* branches = this->GetBranches( );
396 // Reconstruct pixels
397 typename TUniqueVertices::ConstIterator eIt = endpoints->Begin( );
398 for( ; eIt != endpoints->End( ); ++eIt )
400 // Get the path all the way to global seed
402 mst->GetPath( path, *eIt, seed );
404 TVertex start_vertex = *eIt;
405 typename TVertices::const_iterator pIt = path.begin( );
406 for( ; pIt != path.end( ); ++pIt )
408 if( bifurcations->Find( *pIt ) )
410 branches->SetValue( start_vertex, *pIt, 1 );
417 // Finish with branches to global seed
418 branches->SetValue( start_vertex, seed, 1 );
423 // -------------------------------------------------------------------------
424 template< class I, class O >
425 void fpa::Image::DijkstraWithEndPointDetection< I, O >::
428 // Configure and obtain inputs
429 const I* input = this->GetInput( );
430 const TMinimumSpanningTree* mst = this->GetMinimumSpanningTree( );
431 TBranches* branches = this->GetBranches( );
432 TLabelImage* label = this->GetLabelImage( );
433 label->FillBuffer( 0 );
436 typename TBranches::Iterator bIt = branches->Begin( );
437 TLabel actual_label = 0;
438 for( ; bIt != branches->End( ); ++bIt )
440 typename TBranches::RowIterator brIt = branches->Begin( bIt );
441 for( ; brIt != branches->End( bIt ); ++brIt )
444 brIt->second = actual_label;
447 mst->GetPath( path, bIt->first, brIt->first );
448 typename TVertices::const_iterator pIt = path.begin( );
449 for( ; pIt != path.end( ); ++pIt )
451 double d = std::sqrt( double( input->GetPixel( *pIt ) ) );
452 _TRegion region = this->_Region( *pIt, d );
453 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
454 itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
457 for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
459 // Mask real input image
461 !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) &&
462 lIt.Get( ) == TLabel( 0 )
464 lIt.Set( actual_label );
473 this->m_NumberOfBranches = actual_label;
475 // Label bifurcations
476 /* TODO: Is this necessary?
477 typename TUniqueVertices::const_iterator bifIt =
478 this->m_Bifurcations.begin( );
479 for( ; bifIt != this->m_Bifurcations.end( ); ++bifIt )
482 double d = std::sqrt( double( input->GetPixel( *bifIt ) ) );
483 _TRegion region = this->_Region( *bifIt, d );
484 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
485 itk::ImageRegionIteratorWithIndex< TLabelImage > lIt( label, region );
488 for( ; !iIt.IsAtEnd( ); ++iIt, ++lIt )
490 // Mask real input image
491 if( !( iIt.Get( ) < ( typename I::PixelType )( 0 ) ) )
492 lIt.Set( actual_label );
500 // -------------------------------------------------------------------------
501 template< class I, class O >
502 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
503 _TRegion fpa::Image::DijkstraWithEndPointDetection< I, O >::
504 _Region( const TVertex& c, const double& r )
506 typename I::ConstPointer input = this->GetInput( );
507 typename I::SpacingType spac = input->GetSpacing( );
508 _TRegion region = input->GetLargestPossibleRegion( );
509 typename I::IndexType idx0 = region.GetIndex( );
510 typename I::IndexType idx1 = idx0 + region.GetSize( );
512 // Compute region size and index
513 typename I::IndexType i0, i1;
515 for( unsigned int d = 0; d < I::ImageDimension; ++d )
517 // NOTE: 3 is a minimum neighborhood size
519 long( std::ceil( r / double( spac[ d ] ) ) ) +
520 long( this->m_SafetyNeighborhoodSize );
521 i0[ d ] = c[ d ] - s;
522 i1[ d ] = c[ d ] + s;
524 if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
525 if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
526 if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
527 if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
528 size[ d ] = i1[ d ] - i0[ d ];
532 // Prepare region and return it
533 region.SetIndex( i0 );
534 region.SetSize( size );
539 // -------------------------------------------------------------------------
540 template< class I, class O >
542 typename fpa::Image::DijkstraWithEndPointDetection< I, O >::
543 TVertex fpa::Image::DijkstraWithEndPointDetection< I, O >::
544 _MaxInRegion( const _T* image, const TVertex& v, const double& r )
546 typedef itk::ImageRegionConstIteratorWithIndex< _T > _TIt;
548 _TIt iIt( image, this->_Region( v, r ) );
550 TVertex max_vertex = iIt.GetIndex( );
551 typename _T::PixelType max_value = iIt.Get( );
552 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
554 typename _T::PixelType value = iIt.Get( );
555 if( value > max_value )
558 max_vertex = iIt.GetIndex( );
563 return( max_vertex );
566 #endif // __FPA__IMAGE__DIJKSTRAWITHENDPOINTDETECTION__HXX__