lib/cpm/Algorithms/QuadEdge/MeshPlaneCutterFilter.hxx

   1 #ifndef __CPM__ALGORITHMS__QUADEDGE__MESHPLANECUTTERFILTER__HXX__
   2 #define __CPM__ALGORITHMS__QUADEDGE__MESHPLANECUTTERFILTER__HXX__
   3
   4 #include <cmath>
   5 #include <queue>
   6 #include <set>
   7
   8 // -------------------------------------------------------------------------
   9 template< class M >
  10 void cpm::Algorithms::QuadEdge::MeshPlaneCutterFilter< M >::
  11 SetPlaneNormal( const TVector& n )
  12 {
  13   TScalar nn = n.GetNorm( );
  14   if( TScalar( 0 ) < nn )
  15   {
  16     this->m_PlaneNormal = n / nn;
  17     this->Modified( );
  18
  19   } // fi
  20 }
  21
  22 // -------------------------------------------------------------------------
  23 template< class M >
  24 cpm::Algorithms::QuadEdge::MeshPlaneCutterFilter< M >::
  25 MeshPlaneCutterFilter( )
  26   : Superclass( )
  27 {
  28   this->m_PlanePoint.Fill( TScalar( 0 ) );
  29   this->m_PlaneNormal.Fill( TScalar( 0 ) );
  30   this->m_PlaneNormal[ 0 ] = TScalar( 1 );
  31 }
  32
  33 // -------------------------------------------------------------------------
  34 template< class M >
  35 cpm::Algorithms::QuadEdge::MeshPlaneCutterFilter< M >::
  36 ~MeshPlaneCutterFilter( )
  37 {
  38 }
  39
  40 // -------------------------------------------------------------------------
  41 template< class M >
  42 void cpm::Algorithms::QuadEdge::MeshPlaneCutterFilter< M >::
  43 GenerateData( )
  44 {
  45   typedef typename M::TPrimalEdge   _TEdge;
  46   typedef typename M::TQuadEdgeCell _TCell;
  47   typedef typename _TEdge::Iterator _TEdgeIt;
  48
  49   const M* input = this->GetInput( );
  50   M* output = this->GetOutput( );
  51
  52   // Clear output data
  53   output->Initialize( );
  54
  55   // Marked edges
  56   std::set< _TEdge* > marks;
  57
  58   unsigned long nCells = input->GetNumberOfCells( );
  59   unsigned long actualCell = 0;
  60
  61   TVector pp = this->m_PlanePoint;
  62   TVector np = this->m_PlaneNormal;
  63   TVector pl, pm, pml;
  64   TScalar num, den;
  65   while( actualCell < nCells )
  66   {
  67     std::queue< _TEdge* > q;
  68
  69     // Is the next cell a QuadEdgeCell?
  70     typename M::CellAutoPointer cell_ptr;
  71     input->GetCell( actualCell++, cell_ptr );
  72     _TCell* cell = dynamic_cast< _TCell* >( cell_ptr.GetPointer( ) );
  73     if( cell == NULL )
  74       continue;
  75
  76     // Ok, it is... now, does its face needs to be analyzed?
  77     _TEdge* edge = cell->GetEntryPrimalEdge( );
  78     _TEdgeIt fIt = edge->BeginLnext( );
  79     for( ; fIt != edge->EndLnext( ); ++fIt )
  80     {
  81       // Has the edge already been visited?
  82       if( marks.find( *fIt ) != marks.end( ) )
  83         continue;
  84
  85       // No, ok -> intersection candidate
  86       q.push( *fIt );
  87
  88     } // rof
  89
  90     // Main loop
  91     while( !( q.empty( ) ) )
  92     {
  93       // Pick next candidate
  94       edge = q.front( );
  95       q.pop( );
  96
  97       // Check it
  98       if( marks.find( edge ) != marks.end( ) )
  99         continue;
 100
 101       // Mark it (with its symmetric)
 102       marks.insert( edge );
 103       marks.insert( edge->GetSym( ) );
 104
 105       // Now, edge is a true candidate: check intersection
 106       pl = input->GetPoint( edge->GetOrigin( ) ).GetVectorFromOrigin( );
 107       pm = input->GetPoint( edge->GetDestination( ) ).GetVectorFromOrigin( );
 108       pml = pm - pl;
 109       num = ( pp - pl ) * pml;
 110       den = pml * np;
 111       // std::cout << num << " " << den << " " << pp << " " << np << " " << pl << " " << pm << " {" << edge->GetOrigin( ) << "," << edge->GetDestination( ) << "}" << std::endl;
 112       if( !( double( 0 ) < std::fabs( double( den ) ) ) )
 113       {
 114         if( !( double( 0 ) < std::fabs( double( num ) ) ) )
 115         {
 116           // Line lies on plane: both points intersect
 117           // std::cout << pl << " " << pm << " ";
 118
 119         } // fi
 120         // else: Line parallel to plane -> no intersection
 121       }
 122       else
 123       {
 124         num /= den;
 125         if( TScalar( 0 ) <= num && num <= TScalar( 1 ) )
 126         {
 127           // One point intersection
 128           // TPoint pnt( TScalar( 0 )
 129           output->SetPoint(
 130             output->GetNumberOfPoints( ),
 131             ( ( pm - pl ) * num ) + pl
 132             );
 133
 134           // Check over face
 135           fIt = edge->BeginLnext( );
 136           for( ; fIt != edge->EndLnext( ); ++fIt )
 137           {
 138             // Has the edge already been visited?
 139             if( marks.find( *fIt ) != marks.end( ) )
 140               continue;
 141
 142             // No, ok -> intersection candidate
 143             q.push( *fIt );
 144
 145           } // rof
 146         }
 147         // else: no intersection between pm and pl
 148
 149       } // fi
 150
 151     } // elihw
 152
 153   } // elihw
 154 }
 155
 156 #endif // __CPM__ALGORITHMS__QUADEDGE__MESHPLANECUTTERFILTER__HXX__
 157
 158 // eof - $RCSfile$