lib/fpa/Image/SkeletonFilter.hxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #ifndef __fpa__Image__SkeletonFilter__hxx__
   7 #define __fpa__Image__SkeletonFilter__hxx__
   8
   9 #include <itkMinimumMaximumImageCalculator.h>
  10
  11 // -------------------------------------------------------------------------
  12 template< class _TImage, class _TScalar >
  13 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  14 SetInput( const TScalarImage* image )
  15 {
  16   // Do nothing
  17 }
  18
  19 // -------------------------------------------------------------------------
  20 template< class _TImage, class _TScalar >
  21 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  22 SetInput( const TImage* image )
  23 {
  24   this->m_DistanceMap->SetInput( image );
  25   this->Superclass::SetInput( this->m_DistanceMap->GetOutput( ) );
  26 }
  27
  28 // -------------------------------------------------------------------------
  29 template< class _TImage, class _TScalar >
  30 _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  31 GetInput( )
  32 {
  33   return( this->m_DistanceMap->GetInput( ) );
  34 }
  35
  36 // -------------------------------------------------------------------------
  37 template< class _TImage, class _TScalar >
  38 const _TImage* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  39 GetInput( ) const
  40 {
  41   return( this->m_DistanceMap->GetInput( ) );
  42 }
  43
  44 // -------------------------------------------------------------------------
  45 template< class _TImage, class _TScalar >
  46 typename fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  47 TSkeleton* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  48 GetSkeleton( )
  49 {
  50   return(
  51     dynamic_cast< TSkeleton* >(
  52       this->itk::ProcessObject::GetOutput( this->m_SkeletonIdx )
  53       )
  54     );
  55 }
  56
  57 // -------------------------------------------------------------------------
  58 template< class _TImage, class _TScalar >
  59 typename fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  60 TMarks* fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  61 GetMarks( )
  62 {
  63   return(
  64     dynamic_cast< TMarks* >(
  65       this->itk::ProcessObject::GetOutput( this->m_MarksIdx )
  66       )
  67     );
  68 }
  69
  70 // -------------------------------------------------------------------------
  71 template< class _TImage, class _TScalar >
  72 fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  73 SkeletonFilter( )
  74   : Superclass( ),
  75     m_SeedFromMaximumDistance( false )
  76 {
  77   this->m_DistanceMap = TDistanceMap::New( );
  78   this->m_DistanceMap->InsideIsPositiveOn( );
  79   this->m_DistanceMap->SquaredDistanceOff( );
  80   this->m_DistanceMap->UseImageSpacingOn( );
  81
  82   unsigned int nOutputs = this->GetNumberOfRequiredOutputs( );
  83   this->SetNumberOfRequiredOutputs( nOutputs + 2 );
  84   this->m_SkeletonIdx = nOutputs;
  85   this->m_MarksIdx = nOutputs + 1;
  86
  87   typename TSkeleton::Pointer skeleton = TSkeleton::New( );
  88   typename TMarks::Pointer marks = TMarks::New( );
  89   this->SetNthOutput( this->m_SkeletonIdx, skeleton.GetPointer( ) );
  90   this->SetNthOutput( this->m_MarksIdx, marks.GetPointer( ) );
  91 }
  92
  93 // -------------------------------------------------------------------------
  94 template< class _TImage, class _TScalar >
  95 fpa::Image::SkeletonFilter< _TImage, _TScalar >::
  96 ~SkeletonFilter( )
  97 {
  98 }
  99
 100 // -------------------------------------------------------------------------
 101 template< class _TImage, class _TScalar >
 102 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 103 GenerateData( )
 104 {
 105   // Prepare distance map
 106   this->m_DistanceMap->Update( );
 107
 108   // Update start seed
 109   TVertex seed;
 110   if( this->m_SeedFromMaximumDistance )
 111   {
 112     typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax;
 113     typename _TMinMax::Pointer minmax = _TMinMax::New( );
 114     minmax->SetImage( this->m_DistanceMap->GetOutput( ) );
 115     minmax->Compute( );
 116     seed = minmax->GetIndexOfMaximum( );
 117   }
 118   else
 119     seed = *( this->BeginSeeds( ) );
 120   this->m_Seeds.clear( );
 121   this->m_Seeds.insert( seed );
 122
 123   // Prepare skeleton candidates queue
 124   this->m_SkeletonQueue.clear( );
 125
 126   // Go!
 127   this->Superclass::GenerateData( );
 128
 129   // Backtracking
 130   _TAdjacencies A;
 131   std::vector< TVertex > end_points;
 132   this->_EndPoints( end_points, A );
 133   this->_Skeleton( end_points, A );
 134 }
 135
 136 // -------------------------------------------------------------------------
 137 template< class _TImage, class _TScalar >
 138 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 139 _SetOutputValue( const TVertex& vertex, const TOutputValue& value )
 140 {
 141   /* TODO
 142      typedef typename _TSkeletonQueue::value_type _TSkeletonQueueValue;
 143
 144      this->Superclass::_UpdateResult( n );
 145      double d = double( this->GetInput( )->GetPixel( n.Vertex ) );
 146      if( d >= double( 0 ) )
 147      {
 148      // Update skeleton candidates
 149      d += double( 1e-5 );
 150      double v = double( n.Result ) / ( d * d );
 151      this->m_SkeletonQueue.insert( _TSkeletonQueueValue( v, n.Vertex ) );
 152
 153      } // fi
 154   */
 155 }
 156
 157 // -------------------------------------------------------------------------
 158 template< class _TImage, class _TScalar >
 159 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 160 _EndPoints( std::vector< TVertex >& end_points, _TAdjacencies& A )
 161 {
 162   /* TODO
 163      auto marks = this->GetMarks( );
 164      auto mst = this->GetMinimumSpanningTree( );
 165      auto spac = marks->GetSpacing( );
 166
 167      // Some values
 168      marks->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
 169      marks->SetRequestedRegion( mst->GetRequestedRegion( ) );
 170      marks->SetBufferedRegion( mst->GetBufferedRegion( ) );
 171      marks->SetSpacing( mst->GetSpacing( ) );
 172      marks->SetOrigin( mst->GetOrigin( ) );
 173      marks->SetDirection( mst->GetDirection( ) );
 174      marks->Allocate( );
 175      marks->FillBuffer( 0 );
 176
 177      // BFS from maximum queue
 178      while( this->m_SkeletonQueue.size( ) > 0 )
 179      {
 180      // Get node
 181      auto nIt = this->m_SkeletonQueue.begin( );
 182      auto n = *nIt;
 183      this->m_SkeletonQueue.erase( nIt );
 184
 185      // Mark it and update end-points
 186      unsigned char m = marks->GetPixel( n.second );
 187      if( m != 0 )
 188      continue;
 189      marks->SetPixel( n.second, 1 );
 190      end_points.push_back( n.second );
 191
 192      // Mark path
 193      TVertex it = n.second;
 194      TVertex p = mst->GetParent( it );
 195      while( it != p )
 196      {
 197      this->_MarkSphere( it );
 198      it = p;
 199      p = mst->GetParent( it );
 200
 201      } // elihw
 202      this->_MarkSphere( it );
 203      A[ n.second ] = it;
 204
 205      } // elihw
 206   */
 207 }
 208
 209 // -------------------------------------------------------------------------
 210 template< class _TImage, class _TScalar >
 211 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 212 _Skeleton( const std::vector< TVertex >& end_points, _TAdjacencies& A )
 213 {
 214   /* TODO
 215      typedef itk::Image< unsigned long, Self::Dimension > _TTagsImage;
 216      typedef typename TMST::TPath _TPath;
 217
 218      auto mst = this->GetMinimumSpanningTree( );
 219      auto skeleton = this->GetSkeleton( );
 220
 221      // Tag branches
 222      typename _TTagsImage::Pointer tags = _TTagsImage::New( );
 223      tags->SetLargestPossibleRegion( mst->GetLargestPossibleRegion( ) );
 224      tags->SetRequestedRegion( mst->GetRequestedRegion( ) );
 225      tags->SetBufferedRegion( mst->GetBufferedRegion( ) );
 226      tags->Allocate( );
 227      tags->FillBuffer( 0 );
 228      for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
 229      {
 230      TVertex it = *eIt;
 231      TVertex p = mst->GetParent( it );
 232      while( it != p )
 233      {
 234      tags->SetPixel( it, tags->GetPixel( it ) + 1 );
 235      it = p;
 236      p = mst->GetParent( it );
 237
 238      } // elihw
 239      tags->SetPixel( it, tags->GetPixel( it ) + 1 );
 240
 241      } // rof
 242
 243      // Build paths (branches)
 244      for( auto eIt = end_points.begin( ); eIt != end_points.end( ); ++eIt )
 245      {
 246      TVertex it = *eIt;
 247      TVertex p = mst->GetParent( it );
 248      TVertex sIdx = *eIt;
 249      typename _TPath::Pointer path = _TPath::New( );
 250      path->SetReferenceImage( mst );
 251      while( it != p )
 252      {
 253      if( tags->GetPixel( sIdx ) != tags->GetPixel( it ) )
 254      {
 255      // Ok, a new branch has been added
 256      path->AddVertex( it );
 257      skeleton->AddBranch( path );
 258
 259      // Update a new starting index
 260      path = _TPath::New( );
 261      path->SetReferenceImage( mst );
 262      sIdx = it;
 263      }
 264      else
 265      path->AddVertex( it );
 266      it = p;
 267      p = mst->GetParent( it );
 268
 269      } // elihw
 270
 271      // Finally, add last branch
 272      path->AddVertex( it );
 273      skeleton->AddBranch( path );
 274
 275      } // rof
 276   */
 277 }
 278
 279 // -------------------------------------------------------------------------
 280 template< class _TImage, class _TScalar >
 281 void fpa::Image::SkeletonFilter< _TImage, _TScalar >::
 282 _MarkSphere( const TVertex& idx )
 283 {
 284   /* TODO
 285      typedef itk::ImageRegionIteratorWithIndex< TMarks > _TMarksIt;
 286
 287      static const double _eps = std::sqrt( double( Self::Dimension + 1 ) );
 288      auto input = this->GetInput( );
 289      auto marks = this->GetMarks( );
 290      auto spac = input->GetSpacing( );
 291      auto region = input->GetRequestedRegion( );
 292
 293      typename _TImage::PointType cnt;
 294      input->TransformIndexToPhysicalPoint( idx, cnt );
 295      double r = double( input->GetPixel( idx ) ) * _eps;
 296
 297      TVertex i0, i1;
 298      for( unsigned int d = 0; d < Self::Dimension; ++d )
 299      {
 300      long off = long( std::ceil( r / double( spac[ d ] ) ) );
 301      if( off < 3 )
 302      off = 3;
 303      i0[ d ] = idx[ d ] - off;
 304      i1[ d ] = idx[ d ] + off;
 305
 306      if( i0[ d ] < region.GetIndex( )[ d ] )
 307      i0[ d ] = region.GetIndex( )[ d ];
 308
 309      if( i1[ d ] >= region.GetIndex( )[ d ] + region.GetSize( )[ d ] )
 310      i1[ d ] = region.GetIndex( )[ d ] + region.GetSize( )[ d ] - 1;
 311
 312      } // rof
 313
 314      typename _TImage::SizeType size;
 315      for( unsigned int d = 0; d < Self::Dimension; ++d )
 316      size[ d ] = i1[ d ] - i0[ d ] + 1;
 317
 318      typename _TImage::RegionType neighRegion;
 319      neighRegion.SetIndex( i0 );
 320      neighRegion.SetSize( size );
 321
 322      _TMarksIt mIt( marks, neighRegion );
 323      for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt )
 324      {
 325      TVertex w = mIt.GetIndex( );
 326      typename _TImage::PointType p;
 327      marks->TransformIndexToPhysicalPoint( w, p );
 328      mIt.Set( 1 );
 329
 330      } // rof
 331   */
 332 }
 333
 334 #endif // __fpa__Image__SkeletonFilter__hxx__
 335
 336 // eof - $RCSfile$