/*to render 3D bi-cubic Bezier patches*/
#include <GL/glut.h>
#include "bezier.h"

#define BEZTOL 0.000001

static void midPoint(Point3D p, Point3D q, Point3D *m)
{
	m->x = (p.x + q.x)*0.5;
	m->y = (p.y + q.y)*0.5;
	m->z = (p.z + q.z)*0.5;
}

static float squaredDistance(Point3D p, float a, float b, float c, float d)
{
	float dist;
	
	dist = (a*p.x + b*p.y + c*p.z - d);
	
	return (dist*dist)/(a*a + b*b + c*c);
}

static short coplanar(Point3D p[4][4])
/*this tests if all points are approximately on a plane defined by
p[0][0], p[3][0], p[3][3]. Should also really test if edges are straight*/
{
	Point3D pt[3];
	Vector3D v[2], n;
	float a,b,c,d;
	int i,j;
	
	pt[0] = p[0][0];
	pt[1] = p[3][0];
	pt[2] = p[3][3];

	/*find the two difference vectors*/
	for(i=0;i<2;++i) differencePoint3D(&pt[i+1],&pt[0],&v[i]);

	/*find the cross product and store in n*/
	crossProductVector3D(&v[0],&v[1],&n);

	/*store required values in the plane equation*/
	a = n.x; b = n.y; c = n.z;
	d = dotProductVector3D(&n, (Vector3D *)&pt[0]);
	
	for(i=0;i<4;++i)
	for(j=0;j<4;++j)
		if(squaredDistance(p[i][j],a,b,c,d) > BEZTOL) return 0;
	
	return 1;
}
		

static void splitRow(Point3D p0, Point3D p1, Point3D p2, Point3D p3, Point3D r[7])
{
	Point3D s;

	r[0] = p0;
	
	midPoint(p0,p1,&r[1]);
	
	midPoint(p1,p2,&s);
	midPoint(r[1],s,&r[2]);
	
	midPoint(p2,p3,&r[5]);
	midPoint(s,r[5],&r[4]);
	
	midPoint(r[2],r[4],&r[3]);
	
	r[6] = p3;
}
	
	

static void split(Point3D p[4][4], 
		   Point3D q[4][4], Point3D r[4][4], 
		   Point3D s[4][4], Point3D t[4][4])
{
	Point3D row[4][7], a[7][7];
	int i,j;
	
	/*split each row*/
	for(i=0;i<4;++i) splitRow(p[i][0],p[i][1],p[i][2],p[i][3],row[i]);
	
	/*split each columns*/
	for(j=0;j<7;++j) splitRow(row[0][j],row[1][j],row[2][j],row[3][j],a[j]);
	
	/*now create the 4 output arrays*/
	for(i=0;i<7;++i)
	for(j=0;j<7;++j){
		if(i <= 3 && j <= 3) q[i][j] 	= a[j][i];
		if(i <= 3 && j >= 3) r[i][j-3] 	= a[j][i];
		if(i >= 3 && j <= 3) s[i-3][j] 	= a[j][i];
		if(i >= 3 && j >= 3) t[i-3][j-3] = a[j][i];
	}
}

static void polygon(Point3D p[4][4])
{
	glBegin(GL_LINE_LOOP);
	glVertex3f(p[0][0].x,p[0][0].y,p[0][0].z);
	glVertex3f(p[3][0].x,p[3][0].y,p[3][0].z);
	glVertex3f(p[3][3].x,p[3][3].y,p[3][3].z);
	glVertex3f(p[0][3].x,p[0][3].y,p[0][3].z);
	glEnd();
}

void displayBezier(Point3D p[4][4])
{
	Point3D q[4][4], r[4][4], s[4][4], t[4][4];
	
	if(coplanar(p)) polygon(p);
	else{
		split(p,q,r,s,t);
		displayBezier(q);
		displayBezier(r);
		displayBezier(s);
		displayBezier(t);
	}
}