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 fpa::Image::DijkstraWithSphereBacktracking< I, C >::
11 DijkstraWithSphereBacktracking( )
16 // -------------------------------------------------------------------------
17 template< class I, class C >
18 fpa::Image::DijkstraWithSphereBacktracking< I, C >::
19 ~DijkstraWithSphereBacktracking( )
23 // -------------------------------------------------------------------------
24 template< class I, class C >
25 void fpa::Image::DijkstraWithSphereBacktracking< I, C >::
28 const I* img = this->GetInput( );
29 typename I::SpacingType spac = img->GetSpacing( );
30 double max_spac = spac[ 0 ];
31 for( unsigned int d = 1; d < I::ImageDimension; ++d )
33 ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
34 max_spac *= double( 3 );
37 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
39 _TNode seed = this->m_Seeds[ s ];
40 _TRegion region = this->_Region( seed.Vertex, max_spac );
41 itk::ImageRegionConstIteratorWithIndex< I > iIt( img, region );
44 _TPixel max_value = iIt.Get( );
45 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
47 if( iIt.Get( ) > max_value )
49 seed.Vertex = iIt.GetIndex( );
50 seed.Parent = seed.Vertex;
51 max_value = iIt.Get( );
56 this->m_Seeds[ s ] = seed;
61 this->Superclass::_BeforeMainLoop( );
62 this->m_Candidates.clear( );
65 // -------------------------------------------------------------------------
66 template< class I, class C >
67 void fpa::Image::DijkstraWithSphereBacktracking< I, C >::
70 typedef itk::Image< bool, I::ImageDimension > _TMarkImage;
72 // Finish base algorithm
73 this->Superclass::_AfterMainLoop( );
74 this->m_FinalTree.clear( );
75 this->m_EndPoints.clear( );
76 if( this->m_Candidates.size( ) == 0 )
78 this->InvokeEvent( TEndEvent( ) );
81 const I* input = this->GetInput( );
82 typename I::SpacingType spac = input->GetSpacing( );
83 double max_spac = spac[ 0 ];
84 for( unsigned int d = 1; d < I::ImageDimension; ++d )
86 ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac;
87 max_spac *= double( 3 );
89 // Prepare an object to hold marks
90 typename _TMarkImage::Pointer marks = _TMarkImage::New( );
91 marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
92 marks->SetRequestedRegion( input->GetRequestedRegion( ) );
93 marks->SetBufferedRegion( input->GetBufferedRegion( ) );
94 marks->SetOrigin( input->GetOrigin( ) );
95 marks->SetSpacing( input->GetSpacing( ) );
96 marks->SetDirection( input->GetDirection( ) );
98 marks->FillBuffer( false );
100 // Iterate over the candidates, starting from the candidates that
101 // are near thin branches
102 typename _TCandidates::const_reverse_iterator cIt =
103 this->m_Candidates.rbegin( );
104 for( int backId = 0; cIt != this->m_Candidates.rend( ); ++cIt )
106 // If pixel hasn't been visited, continue
107 TVertex v = cIt->second;
108 if( marks->GetPixel( v ) )
111 // Compute nearest start candidate
112 _TRegion region = this->_Region( v, max_spac );
113 itk::ImageRegionConstIteratorWithIndex< I > iIt( input, region );
115 TVertex max_vertex = iIt.GetIndex( );
116 _TPixel max_value = iIt.Get( );
117 for( ++iIt; !iIt.IsAtEnd( ); ++iIt )
119 _TPixel value = iIt.Get( );
120 if( value > max_value )
123 max_vertex = iIt.GetIndex( );
130 if( marks->GetPixel( max_vertex ) )
133 this->m_EndPoints.push_back( max_vertex );
135 // Construct path (at least the part that hasn't been iterated)
138 if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.end( ) )
140 // Mark a sphere around current point as visited
141 double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) );
142 region = this->_Region( max_vertex, dist * double( 1.25 ) );
143 itk::ImageRegionIteratorWithIndex< _TMarkImage >
144 mIt( marks, region );
145 for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
148 // Next vertex in current path
149 this->InvokeEvent( TBacktrackingEvent( max_vertex, backId ) );
150 this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex );
153 max_vertex = this->_Parent( max_vertex );
155 } while( max_vertex != this->_Parent( max_vertex ) );
156 this->m_FinalTree[ max_vertex ] = max_vertex;
157 this->InvokeEvent( TEndBacktrackingEvent( backId ) );
162 // -------------------------------------------------------------------------
163 template< class I, class C >
164 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
165 _UpdateNeigh( _TNode& nn, const _TNode& n )
167 C nc = this->_Cost( nn.Vertex, n.Vertex );
168 if( TCost( 0 ) < nc )
170 nn.Cost = n.Cost + ( TCost( 1 ) / nc );
178 // -------------------------------------------------------------------------
179 template< class I, class C >
180 bool fpa::Image::DijkstraWithSphereBacktracking< I, C >::
181 _UpdateResult( _TNode& n )
183 bool ret = this->Superclass::_UpdateResult( n );
186 TCost nc = this->_Cost( n.Vertex, n.Parent );
187 if( TCost( 0 ) < nc )
189 TCost cc = n.Cost / nc;
190 this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) );
192 this->GetOutput( )->SetPixel( n.Vertex, cc );
200 // -------------------------------------------------------------------------
201 template< class I, class C >
202 typename fpa::Image::DijkstraWithSphereBacktracking< I, C >::
203 _TRegion fpa::Image::DijkstraWithSphereBacktracking< I, C >::
204 _Region( const TVertex& c, const double& r )
206 typename I::ConstPointer input = this->GetInput( );
207 typename I::SpacingType spac = input->GetSpacing( );
208 _TRegion region = input->GetLargestPossibleRegion( );
209 typename I::IndexType idx0 = region.GetIndex( );
210 typename I::IndexType idx1 = idx0 + region.GetSize( );
212 // Compute region size and index
213 typename I::IndexType i0, i1;
215 for( unsigned int d = 0; d < I::ImageDimension; ++d )
217 long s = long( std::ceil( r / double( spac[ d ] ) ) );
218 i0[ d ] = c[ d ] - s;
219 i1[ d ] = c[ d ] + s;
221 if( i0[ d ] < idx0[ d ] ) i0[ d ] = idx0[ d ];
222 if( i1[ d ] < idx0[ d ] ) i1[ d ] = idx0[ d ];
223 if( i0[ d ] > idx1[ d ] ) i0[ d ] = idx1[ d ];
224 if( i1[ d ] > idx1[ d ] ) i1[ d ] = idx1[ d ];
225 size[ d ] = i1[ d ] - i0[ d ];
229 // Prepare region and return it
230 region.SetIndex( i0 );
231 region.SetSize( size );
235 #endif // __FPA__IMAGE__DIJKSTRAWITHSPHEREBACKTRACKING__HXX__