lib/fpa/Image/SkeletonFilter.hxx

   1 #ifndef __fpa__Image__SkeletonFilter__hxx__
   2 #define __fpa__Image__SkeletonFilter__hxx__
   3
   4 // -------------------------------------------------------------------------
   5 #define __fpa__Image__SkeletonFilter__InputMacro( input_name, input_type, input_id ) \
   6   template< class _TDistanceMap, class _TCostMap >                      \
   7   typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
   8   input_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::  \
   9   Get##input_name( )                                                    \
  10   {                                                                     \
  11     return(                                                             \
  12       dynamic_cast< input_type* >(                                      \
  13         this->Superclass::GetInput( input_id )                          \
  14         )                                                               \
  15       );                                                                \
  16   }                                                                     \
  17   template< class _TDistanceMap, class _TCostMap >                      \
  18   const                                                                 \
  19   typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
  20   input_type*                                                           \
  21   fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::              \
  22   Get##input_name( ) const                                              \
  23   {                                                                     \
  24     return(                                                             \
  25       dynamic_cast< const input_type* >(                                \
  26         this->Superclass::GetInput( input_id )                          \
  27         )                                                               \
  28       );                                                                \
  29   }                                                                     \
  30   template< class _TDistanceMap, class _TCostMap >                      \
  31   void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::         \
  32   Set##input_name( input_type* input )                                  \
  33   {                                                                     \
  34     this->Superclass::SetInput( input_id, input );                      \
  35   }
  36
  37 __fpa__Image__SkeletonFilter__InputMacro( DistanceMap, TDistanceMap, 0 );
  38 __fpa__Image__SkeletonFilter__InputMacro( CostMap, TCostMap, 1 );
  39 __fpa__Image__SkeletonFilter__InputMacro( MinimumSpanningTree, TMST, 2 );
  40
  41 // -------------------------------------------------------------------------
  42 #define __fpa__Image__SkeletonFilter__OutputMacro( output_name, output_type, output_id ) \
  43   template< class _TDistanceMap, class _TCostMap >                      \
  44   typename fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::     \
  45   output_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \
  46   Get##output_name( )                                                   \
  47   {                                                                     \
  48     return(                                                             \
  49       dynamic_cast< output_type* >(                                     \
  50         this->Superclass::GetOutput( output_id )                        \
  51         )                                                               \
  52       );                                                                \
  53   }                                                                     \
  54   template< class _TDistanceMap, class _TCostMap >                      \
  55   const typename                                                        \
  56   fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::              \
  57   output_type* fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >:: \
  58   Get##output_name( ) const                                             \
  59   {                                                                     \
  60     return(                                                             \
  61       dynamic_cast< const output_type* >(                               \
  62         this->Superclass::GetOutput( output_id )                        \
  63         )                                                               \
  64       );                                                                \
  65   }
  66
  67 __fpa__Image__SkeletonFilter__OutputMacro( Skeleton, TSkeleton, 0 );
  68 __fpa__Image__SkeletonFilter__OutputMacro( EndPoints, TIndices, 1 );
  69 __fpa__Image__SkeletonFilter__OutputMacro( Bifurcations, TIndices, 2 );
  70
  71 // -------------------------------------------------------------------------
  72 template< class _TDistanceMap, class _TCostMap >
  73 fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
  74 SkeletonFilter( )
  75   : Superclass( )
  76 {
  77   this->SetNumbeOfRequiredInputs( 3 );
  78   this->SetNumbeOfRequiredOutputs( 3 );
  79
  80   typename TIndices::Pointer end_points = TIndices::New( );
  81   typename TIndices::Pointer bifurcations = TIndices::New( );
  82   typename TSkeleton::Pointer skeleton = TSkeleton::New( );
  83   this->SetNthOutput( 0, skeleton );
  84   this->SetNthOutput( 1, end_points );
  85   this->SetNthOutput( 2, bifurcations );
  86 }
  87
  88 // -------------------------------------------------------------------------
  89 template< class _TDistanceMap, class _TCostMap >
  90 fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
  91 ~SkeletonFilter( )
  92 {
  93 }
  94
  95 // -------------------------------------------------------------------------
  96 template< class _TDistanceMap, class _TCostMap >
  97 void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
  98 GenerateData( )
  99 {
 100   // 0. I/O objects
 101   const TDistanceMap* dmap = this->GetDistanceMap( );
 102   const TCostMap* cmap = this->GetCostMap( );
 103   const TMST* mst = this->GetMinimumSpanningTree( );
 104   TIndices* end_points = this->GetEndPoints( );
 105   TIndices* bifurcations = this->GetBifurcations( );
 106   TSkeleton* skeleton = this->GetSkeleton( );
 107
 108   // 1. Check input correspondance
 109   // TODO
 110
 111   // 2. Detect end-points
 112   this->_EndPoints( dmap, cmap, mst, end_points );
 113
 114   // 3. Build skeleton and keep track of bifurcations
 115   this->_Skeleton( dmap, cmap, mst, end_points, bifurcations, skeleton );
 116 }
 117
 118 // -------------------------------------------------------------------------
 119 template< class _TDistanceMap, class _TCostMap >
 120 void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
 121 _EndPoints(
 122   const TDistanceMap* dmap,
 123   const TCostMap* cmap,
 124   const TMST* mst,
 125   TIndices* end_points
 126   )
 127 {
 128   typedef itk::ImageRegionConstIteratorWithIndex< _TDistanceMap > _TDistMapIt;
 129   typedef itk::ImageRegionConstIteratorWithIndex< _TCostMap > _TCostMapIt;
 130   typedef std::multimap< double, TIndex, std::greater< double > > _TQueue;
 131   typedef typename _TQueue::value_type _TQueueValue;
 132   typedef itk::Image< unsigned char, _TCostMap::ImageDimension > _TMarks;
 133
 134   // Some values
 135   typename _TMarks::Pointer marks = _TMarks::New( );
 136   marks->SetLargestPossibleRegion( dmap->GetLargestPossibleRegion( ) );
 137   marks->SetRequestedRegion( dmap->GetRequestedRegion( ) );
 138   marks->SetBufferedRegion( dmap->GetBufferedRegion( ) );
 139   marks->SetSpacing( dmap->GetSpacing( ) );
 140   marks->SetOrigin( dmap->GetOrigin( ) );
 141   marks->SetDirection( dmap->GetDirection( ) );
 142   marks->Allocate( );
 143   marks->FillBuffer( 0 );
 144
 145   // Create queue
 146   _TQueue queue;
 147   _TDistMapIt dIt( dmap, dmap->GetRequestedRegion( ) );
 148   _TCostMapIt cIt( cmap, cmap->GetRequestedRegion( ) );
 149   dIt.GoToBegin( );
 150   cIt.GoToBegin( );
 151   for( ; !dIt.IsAtEnd( ) && !cIt.IsAtEnd( ); ++dIt, ++cIt )
 152   {
 153     double d = double( dIt.Get( ) );
 154     if( d > 0 )
 155     {
 156       double v = double( cIt.Get( ) ) / ( d * d );
 157       queue.insert( _TQueueValue( v, dIt.GetIndex( ) ) );
 158
 159     } // fi
 160
 161   } // rof
 162
 163   // BFS from maximum queue
 164   while( queue.size( ) > 0 )
 165   {
 166     // Get node
 167     auto nIt = queue.begin( );
 168     auto n = *nIt;
 169     queue.erase( nIt );
 170
 171     unsigned char m = marks->GetPixel( n.second );
 172     if( ( m & 0x01 ) == 0x01 )
 173       continue;
 174
 175     // Mark it and update end-points
 176     m |= 0x01;
 177     marks->SetPixel( n.second, m );
 178     end_points->Get( ).insert( n.second );
 179
 180     // Get path
 181     typename TMST::TPath::Pointer path;
 182     mst->GetPath( path, n.second );
 183     for( unsigned long i = 0; i < path->GetSize( ); ++i )
 184     {
 185       TIndex idx = path->GetVertex( path->GetSize( ) - 1 - i );
 186       typename _TCostMap::PointType cnt;
 187       cmap->TransformIndexToPhysicalPoint( idx, cnt );
 188       double d = double( dmap->GetPixel( idx ) );
 189       d *= double( 2 );
 190
 191       // Mark sphere
 192       std::queue< TIndex > q;
 193       q.push( idx );
 194       while( q.size( ) > 0 )
 195       {
 196         TIndex v = q.front( );
 197         q.pop( );
 198         unsigned char m = marks->GetPixel( v );
 199         if( ( m & 0x02 ) == 0x02 )
 200           continue;
 201         m |= 0x03;
 202         marks->SetPixel( v, m );
 203
 204         for( unsigned int x = 0; x < _TCostMap::ImageDimension; ++x )
 205         {
 206           for( int y = -1; y <= 1; y += 2 )
 207           {
 208             TIndex w = v;
 209             w[ x ] += y;
 210             if( region.IsInside( w ) )
 211             {
 212               typename _TCostMap::PointType p;
 213               cmap->TransformIndexToPhysicalPoint( w, p );
 214               if( cnt.EuclideanDistanceTo( p ) <= d )
 215                 q.push( w );
 216
 217             } // fi
 218
 219           } // rof
 220
 221         } // rof
 222
 223       } // elihw
 224
 225     } // rof
 226
 227   } // elihw
 228 }
 229
 230 // -------------------------------------------------------------------------
 231 template< class _TDistanceMap, class _TCostMap >
 232 void fpa::Image::SkeletonFilter< _TDistanceMap, _TCostMap >::
 233 _Skeleton(
 234   const TDistanceMap* dmap,
 235   const TCostMap* cmap,
 236   const TMST* mst,
 237   TIndices* end_points
 238   TIndices* bifurcations,
 239   TSkeleton* skeleton
 240   )
 241 {
 242 }
 243
 244 #endif // __fpa__Image__SkeletonFilter__hxx__
 245
 246 // eof - $RCSfile$