1 #ifndef __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
2 #define __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__
5 #include <itkImageRegionConstIteratorWithIndex.h>
6 #include <itkImageRegionIteratorWithIndex.h>
8 // -------------------------------------------------------------------------
9 template< class I, class C >
10 typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
11 TMarkImage* fpa::Image::DijkstraWithSphereBacktracking< I, C >::
15 dynamic_cast< TMarkImage* >(
16 this->itk::ProcessObject::GetOutput( 1 )
21 // -------------------------------------------------------------------------
22 template< class I, class C >
23 const typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
24 TMarkImage* fpa::Image::DijkstraWithSphereBacktracking< I, C >::
25 GetOutputMarkImage( ) const
28 dynamic_cast< const TMarkImage* >(
29 this->itk::ProcessObject::GetOutput( 1 )
34 // -------------------------------------------------------------------------
35 template< class I, class C >
36 fpa::Image::DijkstraWithSphereBacktracking< I, C >::
37 DijkstraWithSphereBacktracking( )
39 m_NumberOfBranches( 0 )
41 this->SetNumberOfRequiredOutputs( 2 );
42 this->SetNthOutput( 0, I::New( ) );
43 this->SetNthOutput( 1, TMarkImage::New( ) );
46 // -------------------------------------------------------------------------
47 template< class I, class C >
48 fpa::Image::DijkstraWithSphereBacktracking< I, C >::
49 ~DijkstraWithSphereBacktracking( )
53 // -------------------------------------------------------------------------
54 template< class I, class C >
55 void fpa::Image::DijkstraWithSphereBacktracking< I, C >::
58 const I* img = this->GetInput( );
59 typename I::SpacingType spac = img->GetSpacing( );
60 double max_spac = spac[ 0 ];
61 for( unsigned int d = 1; d < I::ImageDimension; ++d )
63 ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
64 max_spac *= double( 30 );
67 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
69 _TNode seed = this->m_Seeds[ s ];
70 _TRegion region = this->_Region( seed.Vertex, max_spac );
71 itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region );
74 _TPixel max_value = iIt.Get( );
75 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
77 if( iIt.Get( ) > max_value )
79 seed.Vertex = iIt.GetIndex( );
80 seed.Parent = seed.Vertex;
81 max_value = iIt.Get( );
86 this->m_Seeds[ s ] = seed;
91 this->Superclass::_BeforeMainLoop( );
92 this->m_Candidates.clear( );
95 // -------------------------------------------------------------------------
96 template< class I, class C >
97 void fpa::Image::DijkstraWithSphereBacktracking< I, C >::
100 // Finish base algorithm
101 this->Superclass::_AfterMainLoop( );
102 this->m_FullTree.clear( );
103 this->m_ReducedTree.clear( );
104 this->m_EndPoints.clear( );
105 this->m_BifurcationPoints.clear( );
106 if( this->m_Candidates.size( ) == 0 )
108 this->InvokeEvent( TEndEvent( ) );
110 // Get some input values
111 const I* input = this->GetInput( );
112 typename I::SpacingType spac = input->GetSpacing( );
113 double max_spac = spac[ 0 ];
114 for( unsigned int d = 1; d < I::ImageDimension; ++d )
116 ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
117 max_spac *= double( 3 );
119 // Prepare an object to hold marks
120 typename TMarkImage::Pointer marks = this->GetOutputMarkImage( );
121 marks->FillBuffer( 0 );
123 // Iterate over the candidates, starting from the candidates that
124 // are near thin branches
125 typename _TCandidates::const_reverse_iterator cIt =
126 this->m_Candidates.rbegin( );
127 this->m_NumberOfBranches = 1;
128 for( ; cIt != this->m_Candidates.rend( ); ++cIt )
130 // If pixel hasn't been visited, continue
131 TVertex v = cIt->second;
132 if( marks->GetPixel( v ) != 0 )
135 // Compute nearest start candidate
136 _TRegion region = this->_Region( v, max_spac );
137 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
139 TVertex max_vertex = iIt.GetIndex( );
140 _TPixel max_value = iIt.Get( );
141 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
143 _TPixel value = iIt.Get( );
144 if( value > max_value )
147 max_vertex = iIt.GetIndex( );
154 if( marks->GetPixel( max_vertex ) != 0 )
157 this->m_EndPoints.push_back( max_vertex );
159 // Construct path (at least the part that hasn't been iterated)
166 if( this->m_FullTree.find( max_vertex ) == this->m_FullTree.end( ) )
168 // Mark a sphere around current point as visited
169 double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
170 region = this->_Region( max_vertex, dist * double( 1.5 ) );
171 itk::ImageRegionIteratorWithIndex< TMarkImage >
172 mIt( marks, region );
173 for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
174 mIt.Set( this->m_NumberOfBranches );
176 // Next vertex in current path
177 this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
178 this->m_FullTree[ max_vertex ] =
179 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
180 // std::cout << "New: " << this->m_NumberOfBranches << std::endl;
184 // A bifurcation point has been reached!
185 this->m_BifurcationPoints.push_back( max_vertex );
186 this->m_NumberOfBranches++;
188 this->m_FullTree[ max_vertex ] =
189 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
190 // std::cout << "Bifurcation: " << this->m_NumberOfBranches << std::endl;
200 this->m_BifurcationPoints.begin( ),
201 this->m_BifurcationPoints.end( ),
203 ) == this->m_BifurcationPoints.end( )
206 // Mark a sphere around current point as visited
207 double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
208 region = this->_Region( max_vertex, dist * double( 1.5 ) );
209 itk::ImageRegionIteratorWithIndex< TMarkImage >
210 mIt( marks, region );
211 for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
212 mIt.Set( this->m_NumberOfBranches );
214 // Next vertex in current path
215 this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) );
216 this->m_FullTree[ max_vertex ] =
217 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
219 // std::cout << "Change: " << this->m_NumberOfBranches << std::endl;
223 // A bifurcation point has been reached!
224 // TODO: this->m_BifurcationPoints.push_back( max_vertex );
225 this->m_NumberOfBranches++;
226 this->m_FullTree[ max_vertex ] =
227 TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches );
228 // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl;
233 max_vertex = this->_Parent( max_vertex );
235 } while( max_vertex != this->_Parent( max_vertex ) );
236 if( start || change )
237 this->m_NumberOfBranches++;
238 this->m_FullTree[ max_vertex ] = TTreeNode( max_vertex, this->m_NumberOfBranches );
240 this->InvokeEvent( TEndBacktrackingEvent( ) );
244 std::map< TMark, unsigned long > histo;
246 typename TTree::iterator treeIt = this->m_FullTree.begin( );
247 treeIt != this->m_FullTree.end( );
250 histo[ treeIt->second.second ]++;
252 std::map< TMark, TMark > changes;
253 TMark last_change = 1;
254 for( TMark i = 1; i <= this->m_NumberOfBranches; ++i )
256 if( histo[ i ] != 0 )
257 changes[ i ] = last_change++;
260 this->m_NumberOfBranches = changes.size( );
263 typename TTree::iterator treeIt = this->m_FullTree.begin( );
264 treeIt != this->m_FullTree.end( );
268 TMark old = treeIt->second.second;
269 treeIt->second.second = changes[ old ];
273 // Construct reduced tree
275 typename TVertices::const_iterator eIt = this->m_EndPoints.begin( );
276 eIt != this->m_EndPoints.end( );
280 typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt );
281 if( tIt != this->m_FullTree.end( ) )
283 TMark id = tIt->second.second;
286 tIt = this->m_FullTree.find( tIt->second.first );
287 if( tIt == this->m_FullTree.end( ) )
290 } while( tIt->second.second == id );
291 this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id );
298 typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( );
299 bIt != this->m_BifurcationPoints.end( );
303 typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt );
304 if( tIt != this->m_FullTree.end( ) )
306 TMark id = tIt->second.second;
309 tIt = this->m_FullTree.find( tIt->second.first );
310 if( tIt == this->m_FullTree.end( ) )
313 } while( tIt->second.second == id );
314 this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id );
321 // -------------------------------------------------------------------------
322 template< class I, class C >
323 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
324 _UpdateNeigh( _TNode& nn, const _TNode& n )
326 C nc = this->_Cost( nn.Vertex, n.Vertex );
327 if( TCost( 0 ) < nc )
329 typename I::PointType pnn, pn;
330 this->GetInput( )->TransformIndexToPhysicalPoint( nn.Vertex, pnn );
331 this->GetInput( )->TransformIndexToPhysicalPoint( n.Vertex, pn );
335 nn.Cost = n.Cost + ( TCost( pnn.EuclideanDistanceTo( pn ) ) / std::pow( nc, 4 ) );
343 // -------------------------------------------------------------------------
344 template< class I, class C >
345 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
346 _UpdateResult( _TNode& n )
348 bool ret = this->Superclass::_UpdateResult( n );
351 TCost nc = this->_Cost( n.Vertex, n.Parent );
352 if( TCost( 0 ) < nc )
354 TCost cc = n.Cost / nc;
355 this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
357 this->GetOutput( )->SetPixel( n.Vertex, cc );
365 // -------------------------------------------------------------------------
366 template< class I, class C >
367 typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
368 _TRegion fpa::Image::DijkstraWithSphereBacktracking< I, C >::
369 _Region( const TVertex& c, const double& r )
371 typename I::ConstPointer input = this->GetInput( );
372 typename I::SpacingType spac = input->GetSpacing( );
373 _TRegion region = input->GetLargestPossibleRegion( );
374 typename I::IndexType idx0 = region.GetIndex( );
375 typename I::IndexType idx1 = idx0 + region.GetSize( );
377 // Compute region size and index
378 typename I::IndexType i0, i1;
380 for( unsigned int d = 0; d < I::ImageDimension; ++d )
382 long s = long( std::ceil( r / double( spac[ d ] ) ) );
383 i0[ d ] = c[ d ] - s;
384 i1[ d ] = c[ d ] + s;
386 if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
387 if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
388 if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
389 if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
390 size[ d ] = i1[ d ] - i0[ d ];
394 // Prepare region and return it
395 region.SetIndex( i0 );
396 region.SetSize( size );
400 #endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__