#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" );
}