#define Sgn(num) (num<0?0:1)

inline vec3f norm(vec3f &p1, vec3f &p2, vec3f &p3)
{
	return (p2-p1).cross(p3-p1);
}

inline bool
check_abcd_vf(vec3f &a0, vec3f &b0, vec3f &c0, vec3f &d0,
			  vec3f &a1, vec3f &b1, vec3f &c1, vec3f &d1)
{
	vec3f n0 = norm(a0, b0, c0);
	vec3f n1 = norm(a1, b1, c1);
	vec3f delta = norm(a1-a0, b1-b0, c1-c0);
	vec3f nX = n0+n1-delta;

	vec3f pa0 = d0-a0;
	vec3f pa1 = d1-a1;

	float A = n0.dot(pa0);
	float B = n1.dot(pa1);
	float C = nX.dot(pa0);
	float D = nX.dot(pa1);
	float E = n1.dot(pa0);
	float F = n0.dot(pa1);

	if (A > 0 && B > 0 && (C+F) > 0 && (D+E) > 0)
		return false;

	if (A < 0 && B < 0 && (C+F) < 0 && (D+E) < 0)
		return false;

	vec3f ba0 = b0-a0;
	vec3f ba1 = b1-a1;
	vec3f ca0 = c0-a0;
	vec3f ca1 = c1-a1;

	vec3f A1,B1,C1,D1,E1,F1;
	float M1, M2, M3, M4, M5, N1, N2, N3, N4, N5;
	A1 = ba0.cross(n0);
	B1 = ba1.cross(n1);
	C1 = ba0.cross(nX);
	D1 = ba1.cross(nX);
	E1 = ba0.cross(n1);
	F1 = ba1.cross(n0);
	M1 = A1.dot(ca0);
	M2 = B1.dot(ca1);
	M3 = (C1+F1).dot(ca0)+A1.dot(ca1);
	M4 = (D1+E1).dot(ca0)+(C1+F1).dot(ca1);
	M5 = B1.dot(ca0)+(D1+E1).dot(ca1);

	if(((Sgn(M1)^Sgn(M2)) | (Sgn(M1)^Sgn(M3)) | (Sgn(M1)^Sgn(M4)) | (Sgn(M1)^Sgn(M5))) == 0)
	{
		N1 = A1.dot(pa0);
		N2 = B1.dot(pa1);
		N3 = (C1+F1).dot(pa0)+A1.dot(pa1);
		N4 = (D1+E1).dot(pa0)+(C1+F1).dot(pa1);
		N5 = B1.dot(pa0)+(D1+E1).dot(pa1);
		if((Sgn(M1)^Sgn(N1)) & (Sgn(M1)^Sgn(N2)) & (Sgn(M1)^Sgn(N3)) & (Sgn(M1)^Sgn(N4)) & (Sgn(M1)^Sgn(N5)))
			return false;

		A1 = n0.cross(ca0);
		B1 = n1.cross(ca1);
		C1 = nX.cross(ca0);
		D1 = nX.cross(ca1);
		E1 = n1.cross(ca0);
		F1 = n0.cross(ca1);
		N1 = A1.dot(pa0);
		N2 = B1.dot(pa1);
		N3 = (C1+F1).dot(pa0)+A1.dot(pa1);
		N4 = (D1+E1).dot(pa0)+(C1+F1).dot(pa1);
		N5 = B1.dot(pa0)+(D1+E1).dot(pa1);
		if((Sgn(M1)^Sgn(N1)) & (Sgn(M1)^Sgn(N2)) & (Sgn(M1)^Sgn(N3)) & (Sgn(M1)^Sgn(N4)) & (Sgn(M1)^Sgn(N5)))
			return false;

		vec3f cb0 = c0-b0;
		vec3f cb1 = c1-b1;
		vec3f db0 = d0-b0;
		vec3f db1 = d1-b1;
		A1 = cb0.cross(n0);
		B1 = cb1.cross(n1);
		C1 = cb0.cross(nX);
		D1 = cb1.cross(nX);
		E1 = cb0.cross(n1);
		F1 = cb1.cross(n0);
		N1 = A1.dot(db0);
		N2 = B1.dot(db1);
		N3 = (C1+F1).dot(db0)+A1.dot(db1);
		N4 = (D1+E1).dot(db0)+(C1+F1).dot(db1);
		N5 = B1.dot(db0)+(D1+E1).dot(db1);
		if((Sgn(M1)^Sgn(N1)) & (Sgn(M1)^Sgn(N2)) & (Sgn(M1)^Sgn(N3)) & (Sgn(M1)^Sgn(N4)) & (Sgn(M1)^Sgn(N5)))
			return false;
	}
	return true;
}

inline bool
check_abcd_ee(vec3f &a0, vec3f &b0, vec3f &c0, vec3f &d0,
			  vec3f &a1, vec3f &b1, vec3f &c1, vec3f &d1)
{
	vec3f n0 = norm(a0, b0, c0);
	vec3f n1 = norm(a1, b1, c1);
	vec3f delta = norm(a1-a0, b1-b0, c1-c0);
	vec3f nX = n0+n1-delta;

	vec3f pa0 = d0-a0;
	vec3f pa1 = d1-a1;

	float A = n0.dot(pa0);
	float B = n1.dot(pa1);
	float C = nX.dot(pa0);
	float D = nX.dot(pa1);
	float E = n1.dot(pa0);
	float F = n0.dot(pa1);

	if (A > 0 && B > 0 && (C+F) > 0 && (D+E) > 0)
		return false;

	if (A < 0 && B < 0 && (C+F) < 0 && (D+E) < 0)
		return false;

	vec3f ba0 = b0-a0;
	vec3f ba1 = b1-a1;
	vec3f ca0 = c0-a0;
	vec3f ca1 = c1-a1;

	vec3f A1,B1,C1,D1,E1,F1;
	float M1, M2, M3, M4, M5, N1, N2, N3, N4, N5;
	A1 = ba0.cross(n0);
	B1 = ba1.cross(n1);
	C1 = ba0.cross(nX);
	D1 = ba1.cross(nX);
	E1 = ba0.cross(n1);
	F1 = ba1.cross(n0);

	M1 = A1.dot(ca0);
	N1 = A1.dot(pa0);
	M2 = B1.dot(ca1);
	N2 = B1.dot(pa1);

	if(((Sgn(M1)^Sgn(N1)) & (Sgn(M1)^Sgn(N2)) & (Sgn(M1)^Sgn(M2)))==0)
	{
		M3 = (C1+F1).dot(ca0)+A1.dot(ca1);
		N3 = (C1+F1).dot(pa0)+A1.dot(pa1);
		M4 = (D1+E1).dot(ca0)+(C1+F1).dot(ca1);
		N4 = (D1+E1).dot(pa0)+(C1+F1).dot(pa1);
		M5 = B1.dot(ca0)+(D1+E1).dot(ca1);
		N5 = B1.dot(pa0)+(D1+E1).dot(pa1);
		if (((Sgn(M1)^Sgn(M3)) & (Sgn(M1)^Sgn(M4)) & (Sgn(M1)^Sgn(M5)) & (Sgn(M1)^Sgn(N3)) & (Sgn(M1)^Sgn(N4)) & (Sgn(M1)^Sgn(N5)))==0)
			return false;
	}

	n0 = norm(b0, c0, d0);
	n1 = norm(b1, c1, d1);
	delta = norm(b1-b0, c1-c0, d1-d0);
	nX = n0+n1-delta;
	vec3f pc0 = d0-c0;
	vec3f pc1 = d1-c1;
	vec3f bc0 = b0-c0;
	vec3f bc1 = b1-c1;
	A1 = pc0.cross(n0);
	B1 = pc1.cross(n1);
	C1 = pc0.cross(nX);
	D1 = pc1.cross(nX);
	E1 = pc0.cross(n1);
	F1 = pc1.cross(n0);

	M1 = A1.dot(a0-c0);
	N1 = A1.dot(bc0);
	M2 = B1.dot(a1-c1);
	N2 = B1.dot(bc1);

	if(((Sgn(M1)^Sgn(N1)) & (Sgn(M1)^Sgn(N2)) & (Sgn(M1)^Sgn(M2)))==0)
	{

		M3 = (C1+F1).dot(ca0)+A1.dot(ca1);
		N3 = (C1+F1).dot(bc0)+A1.dot(bc1);
		M4 = (D1+F1).dot(ca0)+(C1+F1).dot(ca1);
		N4 = (D1+F1).dot(bc0)+(C1+F1).dot(bc1);

		M5 = B1.dot(ca0)+(D1+E1).dot(ca1);
		N5 = B1.dot(bc0)+(D1+E1).dot(bc1);

		if (((Sgn(M1)^Sgn(M3)) & (Sgn(M1)^Sgn(M4)) & (Sgn(M1)^Sgn(M5)) & (Sgn(M1)^Sgn(N3)) & (Sgn(M1)^Sgn(N4)) & (Sgn(M1)^Sgn(N5)))==0)
			return false;
	}
	return true;
}

bool
DeformModel::check_vf(unsigned int fid1, unsigned int vid2)
{
	unsigned v0 = _tris[fid1].id0();
	unsigned v1 = _tris[fid1].id1();
	unsigned v2 = _tris[fid1].id2();

	vec3f &a0 = _prev_vtxs[v0];
	vec3f &b0 = _prev_vtxs[v1];
	vec3f &c0 = _prev_vtxs[v2];
	vec3f &p0 = _prev_vtxs[vid2];

	vec3f &a1 = _cur_vtxs[v0];
	vec3f &b1 = _cur_vtxs[v1];
	vec3f &c1 = _cur_vtxs[v2];
	vec3f &p1 = _cur_vtxs[vid2];

	bool bo = check_abcd_vf(a0, b0, c0, p0, a1, b1, c1, p1);
}

bool
DeformModel::check_ee(unsigned int e1, unsigned int e2)
{
	unsigned v0 = _edges[e1].vid(0);
	unsigned v1 = _edges[e1].vid(1);
	unsigned w0 = _edges[e2].vid(0);
	unsigned w1 = _edges[e2].vid(1);

	vec3f &a0 = _prev_vtxs[v0];
	vec3f &b0 = _prev_vtxs[v1];
	vec3f &c0 = _prev_vtxs[w0];
	vec3f &d0 = _prev_vtxs[w1];

	vec3f &a1 = _cur_vtxs[v0];
	vec3f &b1 = _cur_vtxs[v1];
	vec3f &c1 = _cur_vtxs[w0];
	vec3f &d1 = _cur_vtxs[w1];

	bool bo = check_abcd_ee(a0, b0, c0, d0, a1, b1, c1, d1);
}