1 #ifndef __fpa__Image__EndPointsFilter__hxx__
2 #define __fpa__Image__EndPointsFilter__hxx__
7 #include <itkImageRegionConstIteratorWithIndex.h>
9 // -------------------------------------------------------------------------
10 template< class _TDistanceMap, class _TCostMap >
11 const typename fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
12 TIndices& fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
13 GetBifurcations( ) const
15 return( this->m_Bifurcations );
18 // -------------------------------------------------------------------------
19 template< class _TDistanceMap, class _TCostMap >
20 const typename fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
21 TIndices& fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
24 return( this->m_EndPoints );
27 // -------------------------------------------------------------------------
28 template< class _TDistanceMap, class _TCostMap >
29 void fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
32 typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
33 typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
34 typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
35 typedef typename _TQueue::value_type _TQueueValue;
36 typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks;
39 typename _TDistanceMap::RegionType region =
40 this->m_DistanceMap->GetRequestedRegion( );
41 typename _TMarks::Pointer marks = _TMarks::New( );
42 marks->SetLargestPossibleRegion( this->m_DistanceMap->GetLargestPossibleRegion( ) );
43 marks->SetRequestedRegion( this->m_DistanceMap->GetRequestedRegion( ) );
44 marks->SetBufferedRegion( this->m_DistanceMap->GetBufferedRegion( ) );
45 marks->SetSpacing( this->m_DistanceMap->GetSpacing( ) );
46 marks->SetOrigin( this->m_DistanceMap->GetOrigin( ) );
47 marks->SetDirection( this->m_DistanceMap->GetDirection( ) );
49 marks->FillBuffer( 0 );
53 _TDistMapIt dIt( this->m_DistanceMap, region );
54 _TCostMapIt cIt( this->m_CostMap, region );
57 for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
59 double d = double( dIt.Get( ) );
62 double v = double( cIt.Get( ) ) / ( d * d );
63 queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
69 // BFS from maximum queue
70 while( queue.size( ) > 0 )
73 auto nIt = queue.begin( );
77 unsigned char m = marks->GetPixel( n.second );
78 if( ( m & 0x01 ) == 0x01 )
83 marks->SetPixel( n.second, m );
84 this->m_EndPoints.insert( n.second );
87 typename TMST::TPath::Pointer path;
88 this->m_MST->GetPath( path, n.second );
89 for( unsigned long i = 0; i < path->GetSize( ); ++i )
91 TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i );
92 typename _TCostMap::PointType cnt;
93 this->m_CostMap->TransformIndexToPhysicalPoint( idx, cnt );
94 double d = double( this->m_DistanceMap->GetPixel( idx ) );
98 std::queue< TIndex > q;
100 while( q.size( ) > 0 )
102 TIndex v = q.front( );
104 unsigned char m = marks->GetPixel( v );
105 if( ( m & 0x02 ) == 0x02 )
108 marks->SetPixel( v, m );
110 for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
112 for( int y = -1; y <= 1; y += 2 )
116 if( region.IsInside( w ) )
118 typename _TCostMap::PointType p;
119 this->m_CostMap->TransformIndexToPhysicalPoint( w, p );
120 if( cnt.EuclideanDistanceTo( p ) <= d )
136 // -------------------------------------------------------------------------
137 template< class _TDistanceMap, class _TCostMap >
138 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
144 // -------------------------------------------------------------------------
145 template< class _TDistanceMap, class _TCostMap >
146 fpa::Image::EndPointsFilter< _TDistanceMap, _TCostMap >::
151 #endif // __fpa__Image__EndPointsFilter__hxx__