Check out the Pastebin Gadgets Shop. We have thousands of fun, geeky & affordable gadgets on sale :-)Want more features on Pastebin? Sign Up, it's FREE!
tweet

# Untitled

By: a guest on Jul 31st, 2011  |  syntax: C++  |  size: 5.13 KB  |  views: 100  |  expires: Never
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;