lib/fpa/Base/ExtractEndPointsAndBifurcationsFromMinimumSpanningTree.hxx

   1 #ifndef __FPA__BASE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
   2 #define __FPA__BASE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
   3
   4 #include <cmath>
   5 #include <map>
   6 #include <vector>
   7
   8 // -------------------------------------------------------------------------
   9 template< class _TMST >
  10 const typename
  11 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  12 TMinimumSpanningTree*
  13 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  14 GetMinimumSpanningTree( )
  15 {
  16   return(
  17     dynamic_cast< const _TMST* >( this->itk::ProcessObject::GetInput( 0 ) )
  18     );
  19 }
  20
  21 // -------------------------------------------------------------------------
  22 template< class _TMST >
  23 void
  24 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  25 SetMinimumSpanningTree( TMinimumSpanningTree* mst )
  26 {
  27   this->itk::ProcessObject::SetNthInput( 0, const_cast< _TMST* >( mst ) );
  28 }
  29
  30 // -------------------------------------------------------------------------
  31 template< class _TMST >
  32 typename
  33 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  34 TVertices*
  35 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  36 GetEndPoints( )
  37 {
  38   return(
  39     dynamic_cast< TVertices* >( this->itk::ProcessObject::GetOutput( 0 ) )
  40     );
  41 }
  42
  43 // -------------------------------------------------------------------------
  44 template< class _TMST >
  45 typename
  46 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  47 TVertices*
  48 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  49 GetBifurcations( )
  50 {
  51   return(
  52     dynamic_cast< TVertices* >( this->itk::ProcessObject::GetOutput( 1 ) )
  53     );
  54 }
  55
  56 // -------------------------------------------------------------------------
  57 template< class _TMST >
  58 typename
  59 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  60 TVertices*
  61 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  62 GetCollisions( )
  63 {
  64   return(
  65     dynamic_cast< TVertices* >( this->itk::ProcessObject::GetOutput( 2 ) )
  66     );
  67 }
  68
  69 // -------------------------------------------------------------------------
  70 template< class _TMST >
  71 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  72 ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
  73   : Superclass( )
  74 {
  75   this->SetNumberOfRequiredInputs( 1 );
  76   this->SetNumberOfRequiredOutputs( 3 );
  77   typename TVertices::Pointer ep = TVertices::New( );
  78   typename TVertices::Pointer bf = TVertices::New( );
  79   typename TVertices::Pointer co = TVertices::New( );
  80   this->itk::ProcessObject::SetNthOutput( 0, ep.GetPointer( ) );
  81   this->itk::ProcessObject::SetNthOutput( 1, bf.GetPointer( ) );
  82   this->itk::ProcessObject::SetNthOutput( 2, co.GetPointer( ) );
  83 }
  84
  85 // -------------------------------------------------------------------------
  86 template< class _TMST >
  87 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  88 ~ExtractEndPointsAndBifurcationsFromMinimumSpanningTree( )
  89 {
  90 }
  91
  92 // -------------------------------------------------------------------------
  93 template< class _TMST >
  94 void
  95 fpa::Base::ExtractEndPointsAndBifurcationsFromMinimumSpanningTree< _TMST >::
  96 GenerateData( )
  97 {
  98   typedef std::multimap< double, TVertex > _TQueue;
  99   typedef std::pair< TVertex, TVertex >    _TBranch;
 100   typedef std::vector< _TBranch >          _TBranches;
 101   typedef typename _TQueue::value_type     _TQueueValue;
 102
 103   // 0. Useful objects and values
 104   auto mst = this->GetMinimumSpanningTree( );
 105   auto endpoints = this->GetEndPoints( );
 106   auto bifurcations = this->GetBifurcations( );
 107   auto collisions = this->GetCollisions( );
 108   _TBranches branches;
 109   endpoints->Get( ).clear( );
 110   bifurcations->Get( ).clear( );
 111   collisions->Get( ).clear( );
 112
 113   // 1. Get priority queue
 114   auto& q = mst->GetNodeQueue( );
 115
 116   // 2. Main loop: iterate in inverse order (max-priority)
 117   unsigned long label = 0;
 118   for( auto qIt = q.rbegin( ); qIt != q.rend( ); ++qIt )
 119   {
 120     auto gCost = qIt->first;
 121     auto vertex = qIt->second;
 122
 123     // 2.1. Check if the vertex has already been visited
 124     if( this->_Mark( vertex ) > 0 )
 125       continue;
 126
 127     // 2.2. Get path from front seed
 128     auto path = mst->GetPath( vertex );
 129
 130     // 2.3. Prepare new branch data and prepare new end-point
 131     label++;
 132     branches.push_back( _TBranch( vertex, vertex ) );
 133     endpoints->Get( ).push_back( vertex );
 134
 135     // 2.4. Backtracking
 136     auto pIt = path.begin( );
 137     TVertex last_collision;
 138     bool inCollision = false;
 139     while( pIt != path.end( ) )
 140     {
 141       auto mark = this->_SkeletonMark( *pIt );
 142       if( mark > 0 && mark != label )
 143       {
 144         // 2.4.1. A bifurcation has been reached
 145         bifurcations->Get( ).push_back( *pIt );
 146
 147         // Reorder labels
 148         auto coll_branch = branches[ mark ];
 149         branches[ mark  ] = _TBranch( coll_branch.first, *pIt );
 150         branches[ label - 1 ] = _TBranch( qIt->second, *pIt );
 151         branches.push_back( _TBranch( *pIt, coll_branch.second ) );
 152
 153         // Mark skeleton (b,coll_branch_second) with the new label
 154         label++;
 155         while( *pIt != coll_branch.second && pIt != path.end( ) )
 156         {
 157           this->_MarkSkeleton( *pIt, label );
 158           pIt++;
 159
 160         } // elihw
 161
 162         // Force inner loop termination
 163         pIt = path.end( );
 164       }
 165       else
 166       {
 167         if( !inCollision )
 168         {
 169           mark = this->_Mark( *pIt );
 170           if( mark > 0 && mark != label )
 171           {
 172             // 2.4.2. A collision has been reached
 173             last_collision = *pIt;
 174             collisions->Get( ).push_back( *pIt );
 175             std::cout << *pIt << std::endl;
 176             inCollision = true;
 177           }
 178           else
 179           {
 180             // 2.4.3. Just mark sphere and skeleton
 181             double r = this->_Radius( *pIt );
 182             this->_MarkSkeleton( *pIt, label );
 183             this->_MarkSphere( *pIt, r, label );
 184             branches[ label - 1 ].second = *pIt;
 185
 186             // 2.4.4. Is this a seed? -> add it to endpoints
 187             auto sIt = pIt;
 188             sIt++;
 189             if( sIt == path.end( ) )
 190               endpoints->Get( ).push_back( *pIt );
 191
 192           } // fi
 193
 194         } // fi
 195         pIt++;
 196
 197       } // fi
 198
 199     } // elihw
 200
 201   } // rof
 202 }
 203
 204 #endif // __FPA__BASE__EXTRACTENDPOINTSANDBIFURCATIONSFROMMINIMUMSPANNINGTREE__HXX__
 205
 206 // eof - $RCSfile$