Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <stdexcept>
- #include <vector>
- #include <limits>
- #include <algorithm>
- #include <iostream>
- typedef float Field;
- struct Vector {
- Field x,y,z;
- Vector() {}
- Vector( Field x, Field y, Field z ) : x(x),y(y),z(z) {}
- Field scalar( Vector const& v ) const {
- // Scalar production
- return x*v.x + y*v.y + z*v.z;
- }
- };
- // just not to be confused by name
- typedef Vector Point;
- struct Plane : public Vector {
- Field d;
- Plane();
- Plane( Point const& a, Point b, Point c ) {
- // Calculating the plane equation
- b.x -= a.x; b.y -= a.y; b.z -= a.z;
- c.x -= a.x; c.y -= a.y; c.z -= a.z;
- x = b.y*c.z-b.z*c.y;
- y = b.z*c.x-b.x*c.z;
- z = b.x*c.y-b.y*c.x;
- d = - scalar(a);
- }
- bool separates( Point const& a, Point const &b ) {
- // if substutution of two point into plane equation
- // gives different signs the points are separated
- // from each other by the plane
- return (scalar(a)+d)*(scalar(b)+d) < 0;
- }
- };
- struct Face {
- Point p;
- Plane f;
- Face() {}
- Face( Point const& a, Point const& b ) :
- p(a), f(a,b,Point( a.x, a.y, (a.z>10)?-a.z:a.z+10 )) {}
- };
- inline Point planes_interception( Plane const &a, Plane const& b, Plane const& c ) {
- // Calculates point where three planes are crossing each other using Cramer's rule
- 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));
- // If det is zero - there is no solution
- if ( abs(det) < std::numeric_limits<Field>::epsilon() )
- throw std::logic_error("solution does not exist");
- 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);
- 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);
- 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);
- return Point( x/det, y/det, z/det );
- }
- struct Polygon {
- // Points are stored by a circle of faces
- // Each face consists of polygon vertex and vertical plane
- // crossing both stored vertex and the next polygon vertex
- typedef std::vector<Face> FaceSet;
- Polygon( Point const& a, Point const& b, Point const& c ) : eq(a,b,c) {
- faces.push_back( Face(a,b) );
- faces.push_back( Face(b,c) );
- faces.push_back( Face(c,a) );
- }
- void overlap( Polygon p ) {
- if ( faces.empty() || p.faces.empty() ) return;
- Point q_low, q_hi;
- size_t i_low = -1, i_hi = 0;
- for( ; i_hi < faces.size(); ++i_hi ) {
- try { q_hi = planes_interception( faces[i_hi].f, eq, p.eq ); }
- catch( std::exception ) { continue; } // there is no solution
- // checking if the interception point is on polygon boundary
- // using faces order
- Face& prev = p.faces[(i_low-1)%faces.size()];
- Face& next = p.faces[(i_low+1)%faces.size()];
- if ( prev.f.separates(q_hi,next.p) || next.f.separates(q_hi,prev.p) )
- continue;
- // intersection is inside of polygon area
- if ( i_low < 0 ) {
- i_low = i_hi;
- q_low = q_hi;
- } else break;
- }
- Point q = faces[(i_low+i_hi)/2].p;
- if ( i_low < 0 ) {
- // current polygon and new polygon are fully overlapped
- if ( -(q.x*p.eq.x + q.y*p.eq.y + p.eq.d)/p.eq.z > q.z ) {
- faces.clear();
- } else p.faces.clear();
- return;
- }
- // polygons are overlapped only partially
- // checking if current polygon overlapes the new one
- if ( -(q.x*p.eq.x + q.y*p.eq.y + p.eq.d)/p.eq.z > q.z ) {
- // new polygon is above the current one
- // faces in range (i_low+1,i_hi) should be deleted
- FaceSet tail( faces.begin(), faces.begin() + (i_low+1) );
- std::copy( faces.begin() + (i_hi+1), faces.end(), faces.begin() );
- std::copy( tail.begin(), tail.end(), faces.end() - (i_hi+1) );
- faces.erase( faces.end()-(i_hi-i_low), faces.end() );
- } else {
- // new poligon is under the current one
- // faces in range (begin,i_low)U(i_hi+1,end) should be deleted
- std::copy( faces.begin() + (i_low+1), faces.begin() + (i_hi+1), faces.begin() );
- faces.erase( faces.begin()+(i_hi-i_low), faces.end() );
- std::swap( q_low, q_hi );
- }
- // add new faces not braking the circular order
- faces.push_back( Face(q_low,q_hi) );
- faces.push_back( Face(q_hi,faces[0].p) );
- }
- FaceSet faces;
- Plane eq;
- };
- bool is_overlapped( Polygon const& p ) {
- return p.faces.empty();
- }
- struct TopSurface {
- std::vector<Polygon> top;
- void add( Point const& a, Point const& b, Point const& c ) {
- Polygon t(a,b,c);
- for( int i = 0; i < top.size(); i++ ) {
- top[i].overlap( t );
- t.overlap( top[i] );
- }
- top.push_back( t );
- std::remove_if( top.begin(), top.end(), is_overlapped );
- }
- void print() {
- for( int i = 0; i < top.size(); i++ )
- for( int j = 0; j < top[i].faces.size(); j++ ) {
- std::cout << "polygon " << i << std::endl;
- std::cout << "x: " << top[i].faces[j].p.x << std::endl;
- std::cout << "y: " << top[i].faces[j].p.y << std::endl;
- std::cout << "z: " << top[i].faces[j].p.z << std::endl;
- }
- }
- };
- int main() {
- Point a1(0,0,2);
- Point a2(0,1,0);
- Point a3(1,0,0);
- Point b1(0,0,0);
- Point b2(0,1,2);
- Point b3(1,0,2);
- Point c1(0,0,1.5);
- Point c2(0,1,1.5);
- Point c3(1,0,1.5);
- TopSurface s;
- s.add( a1,a2,a3 );
- s.add( b1,b2,b3 );
- s.add( c1,c2,c3 );
- s.print();
- system( "pause" );
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement