lib/fpa/Image/ExtractEndPointsAndBifurcationsFromMinimumSpanningTree.hxx

   1 #ifndef __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
   2 #define __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
   3
   4 #include <cmath>
   5 #include <itkImageRegionIteratorWithIndex.h>
   6
   7 // -------------------------------------------------------------------------
   8 template< class _TImage, class _TMST >
   9 const typename
  10 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  11 TImage*
  12 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  13 GetCostsImage( )
  14 {
  15   return(
  16     dynamic_cast< const _TImage* >( this->itk::ProcessObject::GetInput( 1 ) )
  17     );
  18 }
  19
  20 // -------------------------------------------------------------------------
  21 template< class _TImage, class _TMST >
  22 const typename
  23 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  24 TImage*
  25 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  26 GetDistanceMap( )
  27 {
  28   return(
  29     dynamic_cast< const _TImage* >( this->itk::ProcessObject::GetInput( 2 ) )
  30     );
  31 }
  32
  33 // -------------------------------------------------------------------------
  34 template< class _TImage, class _TMST >
  35 void
  36 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  37 SetCostsImage( TImage* image )
  38 {
  39   this->itk::ProcessObject::SetNthInput( 1, const_cast< _TImage* >( image ) );
  40 }
  41
  42 // -------------------------------------------------------------------------
  43 template< class _TImage, class _TMST >
  44 void
  45 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  46 SetDistanceMap( TImage* image )
  47 {
  48   this->itk::ProcessObject::SetNthInput( 2, const_cast< _TImage* >( image ) );
  49 }
  50
  51 // -------------------------------------------------------------------------
  52 template< class _TImage, class _TMST >
  53 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  54 ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
  55   : Superclass( ),
  56     m_SquaredDistanceMap( false )
  57 {
  58   this->SetNumberOfRequiredInputs( 3 );
  59 }
  60
  61 // -------------------------------------------------------------------------
  62 template< class _TImage, class _TMST >
  63 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  64 ~ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
  65 {
  66 }
  67
  68 // -------------------------------------------------------------------------
  69 template< class _TImage, class _TMST >
  70 void
  71 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
  72 GenerateData( )
  73 {
  74   auto image = this->GetCostsImage( );
  75
  76   this->m_MarkImage = TMarkImage::New( );
  77   this->m_MarkImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
  78   this->m_MarkImage->SetRequestedRegion( image->GetRequestedRegion( ) );
  79   this->m_MarkImage->SetBufferedRegion( image->GetBufferedRegion( ) );
  80   this->m_MarkImage->SetDirection( image->GetDirection( ) );
  81   this->m_MarkImage->SetOrigin( image->GetOrigin( ) );
  82   this->m_MarkImage->SetSpacing( image->GetSpacing( ) );
  83   this->m_MarkImage->Allocate( );
  84   this->m_MarkImage->FillBuffer( 0 );
  85
  86   this->m_SkeletonImage = TMarkImage::New( );
  87   this->m_SkeletonImage->SetLargestPossibleRegion( image->GetLargestPossibleRegion( ) );
  88   this->m_SkeletonImage->SetRequestedRegion( image->GetRequestedRegion( ) );
  89   this->m_SkeletonImage->SetBufferedRegion( image->GetBufferedRegion( ) );
  90   this->m_SkeletonImage->SetDirection( image->GetDirection( ) );
  91   this->m_SkeletonImage->SetOrigin( image->GetOrigin( ) );
  92   this->m_SkeletonImage->SetSpacing( image->GetSpacing( ) );
  93   this->m_SkeletonImage->Allocate( );
  94   this->m_SkeletonImage->FillBuffer( 0 );
  95
  96   this->Superclass::GenerateData( );
  97 }
  98
  99 // -------------------------------------------------------------------------
 100 template< class _TImage, class _TMST >
 101 void
 102 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 103 _MarkSkeleton( const TVertex& v, const unsigned long& l )
 104 {
 105   if( this->m_SkeletonImage->GetRequestedRegion( ).IsInside( v ) )
 106     this->m_SkeletonImage->SetPixel( v, l );
 107 }
 108
 109 // -------------------------------------------------------------------------
 110 template< class _TImage, class _TMST >
 111 void
 112 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 113 _MarkSphere(
 114   const TVertex& v,
 115   const double& r,
 116   const unsigned long& l
 117   )
 118 {
 119   double rr = r * double( 1.5 );
 120
 121   // Get marking region
 122   auto rreg = this->m_MarkImage->GetRequestedRegion( );
 123   auto spac = this->m_MarkImage->GetSpacing( );
 124   TVertex ireg = rreg.GetIndex( );
 125   TVertex jreg = ireg + rreg.GetSize( );
 126   TVertex idx, jdx;
 127   typename _TImage::SizeType size;
 128   bool all = false;
 129   for( unsigned int d = 0; d < _TImage::ImageDimension; ++d )
 130   {
 131     unsigned long s = std::ceil( rr / double( spac[ d ] ) );
 132     if( s < 3 )
 133     {
 134       s = 3;
 135       all = true;
 136
 137     } // fi
 138     s += 3;
 139     idx[ d ] = v[ d ] - s;
 140     jdx[ d ] = v[ d ] + s;
 141
 142     jreg[ d ]--;
 143     if( idx[ d ] < ireg[ d ] ) idx[ d ] = ireg[ d ];
 144     if( idx[ d ] > jreg[ d ] ) idx[ d ] = jreg[ d ];
 145     if( jdx[ d ] < ireg[ d ] ) jdx[ d ] = ireg[ d ];
 146     if( jdx[ d ] > jreg[ d ] ) jdx[ d ] = jreg[ d ];
 147     size[ d ] = jdx[ d ] - idx[ d ] + 1;
 148
 149   } // rof
 150   typename _TImage::RegionType region;
 151   region.SetIndex( idx );
 152   region.SetSize( size );
 153
 154   // Mark region
 155   typename _TImage::PointType c;
 156   this->m_MarkImage->TransformIndexToPhysicalPoint( v, c );
 157   itk::ImageRegionIteratorWithIndex< TMarkImage >
 158     spIt( this->m_MarkImage, region );
 159   for( spIt.GoToBegin( ); !spIt.IsAtEnd( ); ++spIt )
 160   {
 161     typename _TImage::PointType pnt;
 162     this->m_MarkImage->TransformIndexToPhysicalPoint( spIt.GetIndex( ), pnt );
 163     if( double( pnt.EuclideanDistanceTo( c ) ) <= rr || all )
 164       if( spIt.Get( ) == 0 )
 165         spIt.Set( l );
 166
 167   } // rof
 168 }
 169
 170 // -------------------------------------------------------------------------
 171 template< class _TImage, class _TMST >
 172 unsigned long
 173 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 174 _Mark( const TVertex& v )
 175 {
 176   return( this->m_MarkImage->GetPixel( v ) );
 177 }
 178
 179 // -------------------------------------------------------------------------
 180 template< class _TImage, class _TMST >
 181 unsigned long fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 182 _SkeletonMark( const TVertex& v )
 183 {
 184   return( this->m_SkeletonImage->GetPixel( v ) );
 185 }
 186
 187 // -------------------------------------------------------------------------
 188 template< class _TImage, class _TMST >
 189 double
 190 fpa::Image::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TImage, _TMST >::
 191 _Radius( const TVertex& v )
 192 {
 193   double d = double( this->GetDistanceMap( )->GetPixel( v ) );
 194   double s = std::fabs( d );
 195   if( s > double( 0 ) )
 196     s /= d;
 197   if( this->m_SquaredDistanceMap )
 198     return( s * std::sqrt( std::fabs( d ) ) );
 199   else
 200     return( d );
 201 }
 202
 203 #endif // __FPA__IMAGE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
 204
 205 // eof - $RCSfile$