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