1 #ifndef __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
2 #define __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
5 #include <itkImageRegionIteratorWithIndex.h>
7 // -------------------------------------------------------------------------
8 template< class _TImage, class _TMST >
10 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
12 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
16 dynamic_cast< const _TImage* >( this->itk::ProcessObject::GetInput( 1 ) )
20 // -------------------------------------------------------------------------
21 template< class _TImage, class _TMST >
23 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
25 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
29 dynamic_cast< const _TImage* >( this->itk::ProcessObject::GetInput( 2 ) )
33 // -------------------------------------------------------------------------
34 template< class _TImage, class _TMST >
36 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
38 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
42 dynamic_cast< TSkeleton* >(
43 this->itk::ProcessObject::GetOutput(
44 this->GetNumberOfRequiredOutputs( ) - 1
50 // -------------------------------------------------------------------------
51 template< class _TImage, class _TMST >
53 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
54 SetCostsImage( TImage* image )
56 this->itk::ProcessObject::SetNthInput( 1, const_cast< _TImage* >( image ) );
59 // -------------------------------------------------------------------------
60 template< class _TImage, class _TMST >
62 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
63 SetDistanceMap( TImage* image )
65 this->itk::ProcessObject::SetNthInput( 2, const_cast< _TImage* >( image ) );
68 // -------------------------------------------------------------------------
69 template< class _TImage, class _TMST >
70 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
71 ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
73 m_SquaredDistanceMap( false )
75 this->SetNumberOfRequiredInputs( 3 );
76 unsigned int nOuts = this->GetNumberOfRequiredOutputs( );
77 this->SetNumberOfRequiredOutputs( nOuts + 1 );
78 typename TSkeleton::Pointer sk = TSkeleton::New( );
79 this->itk::ProcessObject::SetNthOutput( nOuts, sk.GetPointer( ) );
82 // -------------------------------------------------------------------------
83 template< class _TImage, class _TMST >
84 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
85 ~ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
89 // -------------------------------------------------------------------------
90 template< class _TImage, class _TMST >
92 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
95 // Create auxiliary objects
96 auto image = this->GetCostsImage( );
97 this->m_MarkImage = TMarkImage::New( );
98 this->m_MarkImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
99 this->m_MarkImage->SetRequestedRegion( image->GetRequestedRegion( ) );
100 this->m_MarkImage->SetBufferedRegion( image->GetBufferedRegion( ) );
101 this->m_MarkImage->SetDirection( image->GetDirection( ) );
102 this->m_MarkImage->SetOrigin( image->GetOrigin( ) );
103 this->m_MarkImage->SetSpacing( image->GetSpacing( ) );
104 this->m_MarkImage->Allocate( );
105 this->m_MarkImage->FillBuffer( 0 );
106 this->m_SkeletonImage = TMarkImage::New( );
107 this->m_SkeletonImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
108 this->m_SkeletonImage->SetRequestedRegion( image->GetRequestedRegion( ) );
109 this->m_SkeletonImage->SetBufferedRegion( image->GetBufferedRegion( ) );
110 this->m_SkeletonImage->SetDirection( image->GetDirection( ) );
111 this->m_SkeletonImage->SetOrigin( image->GetOrigin( ) );
112 this->m_SkeletonImage->SetSpacing( image->GetSpacing( ) );
113 this->m_SkeletonImage->Allocate( );
114 this->m_SkeletonImage->FillBuffer( 0 );
117 this->Superclass::GenerateData( );
120 auto sk = this->GetSkeleton( );
121 sk->SetMinimumSpanningTree( this->GetMinimumSpanningTree( ) );
122 auto bIt = this->m_Branches.begin( );
123 for( ; bIt != this->m_Branches.end( ); ++bIt )
124 sk->AddBranch( bIt->first, bIt->second );
127 // -------------------------------------------------------------------------
128 template< class _TImage, class _TMST >
130 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
131 _MarkSkeleton( const TVertex& v, const unsigned long& l )
133 if( this->m_SkeletonImage->GetRequestedRegion( ).IsInside( v ) )
134 this->m_SkeletonImage->SetPixel( v, l );
137 // -------------------------------------------------------------------------
138 template< class _TImage, class _TMST >
140 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
144 const unsigned long& l
147 double rr = r * double( 1.5 );
149 // Get marking region
150 auto rreg = this->m_MarkImage->GetRequestedRegion( );
151 auto spac = this->m_MarkImage->GetSpacing( );
152 TVertex ireg = rreg.GetIndex( );
153 TVertex jreg = ireg + rreg.GetSize( );
155 typename _TImage::SizeType size;
157 for( unsigned int d = 0; d < _TImage::ImageDimension; ++d )
159 unsigned long s = std::ceil( rr / double( spac[ d ] ) );
167 idx[ d ] = v[ d ] - s;
168 jdx[ d ] = v[ d ] + s;
171 if( idx[ d ] < ireg[ d ] ) idx[ d ] = ireg[ d ];
172 if( idx[ d ] > jreg[ d ] ) idx[ d ] = jreg[ d ];
173 if( jdx[ d ] < ireg[ d ] ) jdx[ d ] = ireg[ d ];
174 if( jdx[ d ] > jreg[ d ] ) jdx[ d ] = jreg[ d ];
175 size[ d ] = jdx[ d ] - idx[ d ] + 1;
178 typename _TImage::RegionType region;
179 region.SetIndex( idx );
180 region.SetSize( size );
183 typename _TImage::PointType c;
184 this->m_MarkImage->TransformIndexToPhysicalPoint( v, c );
185 itk::ImageRegionIteratorWithIndex< TMarkImage >
186 spIt( this->m_MarkImage, region );
187 for( spIt.GoToBegin( ); !spIt.IsAtEnd( ); ++spIt )
189 typename _TImage::PointType pnt;
190 this->m_MarkImage->TransformIndexToPhysicalPoint( spIt.GetIndex( ), pnt );
191 if( double( pnt.EuclideanDistanceTo( c ) ) <= rr || all )
192 if( spIt.Get( ) == 0 )
198 // -------------------------------------------------------------------------
199 template< class _TImage, class _TMST >
201 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
202 _Mark( const TVertex& v )
204 return( this->m_MarkImage->GetPixel( v ) );
207 // -------------------------------------------------------------------------
208 template< class _TImage, class _TMST >
209 unsigned long fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
210 _SkeletonMark( const TVertex& v )
212 return( this->m_SkeletonImage->GetPixel( v ) );
215 // -------------------------------------------------------------------------
216 template< class _TImage, class _TMST >
218 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
219 _Radius( const TVertex& v )
221 double d = double( this->GetDistanceMap( )->GetPixel( v ) );
222 double s = std::fabs( d );
223 if( s > double( 0 ) )
225 if( this->m_SquaredDistanceMap )
226 return( s * std::sqrt( std::fabs( d ) ) );
231 #endif // __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__