Want more features on Pastebin? Sign Up, it's FREE!
Guest

Untitled

By: a guest on Jul 31st, 2011  |  syntax: C++  |  size: 5.13 KB  |  views: 95  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include <stdexcept>
  2. #include <vector>
  3. #include <limits>
  4. #include <algorithm>
  5. #include <iostream>
  6.  
  7. typedef float Field;
  8.  
  9. struct Vector {
  10.         Field x,y,z;
  11.         Vector() {}
  12.         Vector( Field x, Field y, Field z ) : x(x),y(y),z(z) {}
  13.         Field scalar( Vector const& v ) const {
  14.                 // Scalar production
  15.                 return x*v.x + y*v.y + z*v.z;
  16.         }
  17. };
  18.  
  19. // just not to be confused by name
  20. typedef Vector Point;
  21.  
  22. struct Plane : public Vector {
  23.         Field d;
  24.         Plane();
  25.         Plane( Point const& a, Point b, Point c ) {
  26.                 // Calculating the plane equation
  27.                 b.x -= a.x; b.y -= a.y; b.z -= a.z;
  28.                 c.x -= a.x; c.y -= a.y; c.z -= a.z;
  29.                 x = b.y*c.z-b.z*c.y;
  30.                 y = b.z*c.x-b.x*c.z;
  31.                 z = b.x*c.y-b.y*c.x;
  32.                 d = - scalar(a);
  33.         }
  34.         bool separates( Point const& a, Point const &b ) {
  35.                 // if substutution of two point into plane equation
  36.                 // gives different signs the points are separated
  37.                 // from each other by the plane
  38.                 return (scalar(a)+d)*(scalar(b)+d) < 0;
  39.         }
  40. };
  41.  
  42. struct Face {
  43.         Point p;
  44.         Plane f;
  45.         Face() {}
  46.         Face( Point const& a, Point const& b ) :
  47.                 p(a), f(a,b,Point( a.x, a.y, (a.z>10)?-a.z:a.z+10 )) {}
  48. };
  49.  
  50. inline Point planes_interception( Plane const &a, Plane const& b, Plane const& c ) {
  51.         // Calculates point where three planes are crossing each other using Cramer's rule
  52.         Field det = -(a.x*(b.y*c.z-b.z*c.y) - a.y*(b.x*c.z-b.z*c.x) + a.z*(b.x*c.y-b.y*c.x));
  53.         // If det is zero - there is no solution
  54.         if ( abs(det) < std::numeric_limits<Field>::epsilon() )
  55.                 throw std::logic_error("solution does not exist");
  56.         Field x = a.d*(b.y*c.z-b.z*c.y) - a.y*(b.d*c.z-b.z*c.d) + a.z*(b.d*c.y-b.y*c.d);
  57.         Field y = a.x*(b.d*c.z-b.z*c.d) - a.d*(b.x*c.z-b.z*c.x) + a.z*(b.x*c.d-b.d*c.x);
  58.         Field z = a.x*(b.y*c.d-b.d*c.y) - a.y*(b.x*c.d-b.d*c.x) + a.d*(b.x*c.y-b.y*c.x);
  59.         return Point( x/det, y/det, z/det );
  60. }
  61.  
  62. struct Polygon {
  63.         // Points are stored by a circle of faces
  64.         // Each face consists of polygon vertex and vertical plane
  65.         // crossing both stored vertex and the next polygon vertex
  66.         typedef std::vector<Face> FaceSet;
  67.         Polygon( Point const& a, Point const& b, Point const& c ) : eq(a,b,c) {
  68.                 faces.push_back( Face(a,b) );
  69.                 faces.push_back( Face(b,c) );
  70.                 faces.push_back( Face(c,a) );
  71.         }
  72.         void overlap( Polygon p ) {
  73.                 if ( faces.empty() || p.faces.empty() ) return;
  74.                 Point q_low, q_hi;
  75.                 size_t i_low = -1, i_hi = 0;
  76.                 for( ; i_hi < faces.size(); ++i_hi ) {
  77.                         try { q_hi = planes_interception( faces[i_hi].f, eq, p.eq ); }
  78.                         catch( std::exception ) { continue; } // there is no solution
  79.                         // checking if the interception point is on polygon boundary
  80.                         // using faces order
  81.                         Face& prev = p.faces[(i_low-1)%faces.size()];
  82.                         Face& next = p.faces[(i_low+1)%faces.size()];
  83.                         if ( prev.f.separates(q_hi,next.p) || next.f.separates(q_hi,prev.p) )
  84.                                 continue;
  85.                         // intersection is inside of polygon area
  86.                         if ( i_low < 0 ) {
  87.                                 i_low = i_hi;
  88.                                 q_low = q_hi;
  89.                         } else break;
  90.                 }
  91.                 Point q = faces[(i_low+i_hi)/2].p;
  92.                 if ( i_low < 0 ) {
  93.                         // current polygon and new polygon are fully overlapped
  94.                         if ( -(q.x*p.eq.x + q.y*p.eq.y + p.eq.d)/p.eq.z > q.z ) {
  95.                                 faces.clear();
  96.                         } else p.faces.clear();
  97.                         return;
  98.                 }
  99.                 // polygons are overlapped only partially
  100.                 // checking if current polygon overlapes the new one
  101.                 if ( -(q.x*p.eq.x + q.y*p.eq.y + p.eq.d)/p.eq.z > q.z ) {
  102.                         // new polygon is above the current one
  103.                         // faces in range (i_low+1,i_hi) should be deleted
  104.                         FaceSet tail( faces.begin(), faces.begin() + (i_low+1) );
  105.                         std::copy( faces.begin() + (i_hi+1), faces.end(), faces.begin() );
  106.                         std::copy( tail.begin(), tail.end(), faces.end() - (i_hi+1) );
  107.                         faces.erase( faces.end()-(i_hi-i_low), faces.end() );
  108.                 } else {
  109.                         // new poligon is under the current one
  110.                         // faces in range (begin,i_low)U(i_hi+1,end) should be deleted
  111.                         std::copy( faces.begin() + (i_low+1), faces.begin() + (i_hi+1), faces.begin() );
  112.                         faces.erase( faces.begin()+(i_hi-i_low), faces.end() );
  113.                         std::swap( q_low, q_hi );
  114.                 }
  115.                 // add new faces not braking the circular order
  116.                 faces.push_back( Face(q_low,q_hi) );
  117.                 faces.push_back( Face(q_hi,faces[0].p) );
  118.         }
  119.         FaceSet faces;
  120.         Plane eq;
  121. };
  122.  
  123. bool is_overlapped( Polygon const& p ) {
  124.         return p.faces.empty();
  125. }
  126.  
  127. struct TopSurface {
  128.         std::vector<Polygon> top;
  129.         void add( Point const& a, Point const& b, Point const& c ) {
  130.                 Polygon t(a,b,c);
  131.                 for( int i = 0; i < top.size(); i++ ) {
  132.                         top[i].overlap( t );
  133.                         t.overlap( top[i] );
  134.                 }
  135.                 top.push_back( t );
  136.                 std::remove_if( top.begin(), top.end(), is_overlapped );
  137.         }
  138.         void print() {
  139.                 for( int i = 0; i < top.size(); i++ )
  140.                         for( int j = 0; j < top[i].faces.size(); j++ ) {
  141.                                 std::cout << "polygon " << i << std::endl;
  142.                                 std::cout << "x: " << top[i].faces[j].p.x << std::endl;
  143.                                 std::cout << "y: " << top[i].faces[j].p.y << std::endl;
  144.                                 std::cout << "z: " << top[i].faces[j].p.z << std::endl;
  145.                         }
  146.         }
  147. };
  148.  
  149. int main() {
  150.         Point a1(0,0,2);
  151.         Point a2(0,1,0);
  152.         Point a3(1,0,0);
  153.  
  154.         Point b1(0,0,0);
  155.         Point b2(0,1,2);
  156.         Point b3(1,0,2);
  157.  
  158.         Point c1(0,0,1.5);
  159.         Point c2(0,1,1.5);
  160.         Point c3(1,0,1.5);
  161.  
  162.         TopSurface s;
  163.         s.add( a1,a2,a3 );
  164.         s.add( b1,b2,b3 );
  165.         s.add( c1,c2,c3 );
  166.         s.print();
  167.         system( "pause" );
  168. }
clone this paste RAW Paste Data