Движение луча света по поверхности Безье / bezier / Curves

/* By Thomas Mшlhave

  This is file contains the functions for two classes, the CCubicCurve and CCubicPatch classes.
	check http://home10.inet.tele.dk/moelhave for more info
*/


#ifndef CURVES_CPP__
#define CURVES_CPP__
#include <windows.h>
#include <GL\gl.h>
#include <math.h>
#include <stdio.h>
#include "matrix.h"
#include "curves.h"

//little function to convert from Point to Matrix4x1
Matrix4x1 ConvertPoint(Point *p)
{
	Matrix4x1 a;
	a.M[0] = p->x;
	a.M[1] = p->y;
	a.M[2] = p->z;
	a.M[3] = 1.0f;
	return a;
}

/*
Function that subdivides a bezier curve, the two new curves are made at t=Ѕ.
Parameters:
L	 : Left subdivided curve storage var
R	 : Right subdivided curve storage var
*/
void CCubicCurve::SubdivideBezierCurve(CCubicCurve* L, CCubicCurve* R)
{
	if (type != TM_BEZIER)
	{
		//error
		L = NULL;
		R = NULL;
		return;
	}
	
	Matrix4x1 L1, L2, L3, L4, H, R1, R2, R3, R4;
	
	L1.M[0] = Gx.M[0];
	L1.M[1] = Gy.M[0];
	L1.M[2] = Gz.M[0];
	
	L2.M[0] = (Gx.M[0]+Gx.M[1]) / 2.0f;
	L2.M[1] = (Gy.M[0]+Gy.M[1]) / 2.0f;
	L2.M[2] = (Gz.M[0]+Gz.M[1]) / 2.0f;
	L2.M[3] = 1.0f;
	
	H.M[0] = (Gx.M[1] + Gx.M[2])/2.0f;
	H.M[1] = (Gy.M[1] + Gy.M[2])/2.0f;
	H.M[2] = (Gz.M[1] + Gz.M[2])/2.0f;
	
	L3.M[0] = (L2.M[0] + H.M[0])/2.0f;
	L3.M[1] = (L2.M[1] + H.M[1])/2.0f;
	L3.M[2] = (L2.M[2] + H.M[2])/2.0f;
	
	R3.M[0] = (Gx.M[2] + Gx.M[3])/2.0f;
	R3.M[1] = (Gy.M[2] + Gy.M[3])/2.0f;
	R3.M[2] = (Gz.M[2] + Gz.M[3])/2.0f;
	
	R2.M[0] = (H.M[0]+R3.M[0])/2.0f;
	R2.M[1] = (H.M[1]+R3.M[1])/2.0f;
	R2.M[2] = (H.M[2]+R3.M[2])/2.0f;
	
	L4.M[0] = (L3.M[0] + R2.M[0])/2.0f;
	L4.M[1] = (L3.M[1] + R2.M[1])/2.0f;
	L4.M[2] = (L3.M[2] + R2.M[2])/2.0f;
	
	R1 = L4;
	
	R4.M[0] = Gx.M[3];
	R4.M[1] = Gy.M[3];
	R4.M[2] = Gz.M[3];
	
	L->GenerateBezierCoefficients(L1, L2, L3, L4);
	R->GenerateBezierCoefficients(R1, R2, R3, R4);
}

void CCubicCurve::GenerateHermiteCoefficients(Matrix4x1 _p1, Matrix4x1 _p2, Matrix4x1 r1, Matrix4x1 r2)
{
	type = TM_HERMITE;
	
	//build hermite base matrix
	Matrix4x4 mh;
	mh.M[0][0] = 2.0f;
	mh.M[0][1] = -2.0f;
	mh.M[0][2] = 1.0f;
	mh.M[0][3] = 1.0f;
	mh.M[1][0] = -3.0f;
	mh.M[1][1] = 3.0f;
	mh.M[1][2] = -2.0f;
	mh.M[1][3] = -1.0f;
	mh.M[2][0] = 0.0f;
	mh.M[2][1] = 0.0f;
	mh.M[2][2] = 1.0f;
	mh.M[2][3] = 0.0f;
	mh.M[3][0] = 1.0f;
	mh.M[3][1] = 0.0f;
	mh.M[3][2] = 0.0f;
	mh.M[3][3] = 0.0f;
	
	//build geometry matrices
	Gx.M[0] = _p1.M[0];
	Gx.M[1] = _p2.M[0];
	Gx.M[2] = r1.M[0];
	Gx.M[3] = r2.M[0];
	
	Gy.M[0] = _p1.M[1];
	Gy.M[1] = _p2.M[1];
	Gy.M[2] = r1.M[1];
	Gy.M[3] = r2.M[1];
	
	Gz.M[0] = _p1.M[2];
	Gz.M[1] = _p2.M[2];
	Gz.M[2] = r1.M[2];
	Gz.M[3] = r2.M[2];
	
	//computer coefficients
	Matrix4x1 Qx = Gx.MultiplyWith4x4(&mh);
	Matrix4x1 Qy = Gy.MultiplyWith4x4(&mh);
	Matrix4x1 Qz = Gz.MultiplyWith4x4(&mh);
	//store coefficients
	coef[0] = Qx.M[0];
	coef[1] = Qx.M[1];
	coef[2] = Qx.M[2];
	coef[3] = Qx.M[3];
	coef[4] = Qy.M[0];
	coef[5] = Qy.M[1];
	coef[6] = Qy.M[2];
	coef[7] = Qy.M[3];
	coef[8] = Qz.M[0];
	coef[9] = Qz.M[1];
	coef[10] = Qz.M[2];
	coef[11] =	 Qz.M[3];
}

void CCubicCurve::GenerateBezierCoefficients(Matrix4x1 p1, Matrix4x1 p2, Matrix4x1 p3, Matrix4x1 p4)
{
	type = TM_BEZIER;
	
	//build hermite base matrix
	Matrix4x4 mh;
	mh.M[0][0] = -1.0f;
	mh.M[0][1] = 3.0f;
	mh.M[0][2] = -3.0f;
	mh.M[0][3] = 1.0f;
	
	mh.M[1][0] = 3.0f;
	mh.M[1][1] = -6.0f;
	mh.M[1][2] = 3.0f;
	mh.M[1][3] = 0.0f;
	
	mh.M[2][0] = -3.0f;
	mh.M[2][1] = 3.0f;
	mh.M[2][2] = 0.0f;
	mh.M[2][3] = 0.0f;
	
	mh.M[3][0] = 1.0f;
	mh.M[3][1] = 0.0f;
	mh.M[3][2] = 0.0f;
	mh.M[3][3] = 0.0f;
	
	//build geometry matrices
	Gx.M[0] = p1.M[0];
	Gx.M[1] = p2.M[0];
	Gx.M[2] = p3.M[0];
	Gx.M[3] = p4.M[0];
	
	Gy.M[0] = p1.M[1];
	Gy.M[1] = p2.M[1];
	Gy.M[2] = p3.M[1];
	Gy.M[3] = p4.M[1];
	
	Gz.M[0] = p1.M[2];
	Gz.M[1] = p2.M[2];
	Gz.M[2] = p3.M[2];
	Gz.M[3] = p4.M[2];
	
	//computer coefficients
	Matrix4x1 Qx = Gx.MultiplyWith4x4(&mh);
	Matrix4x1 Qy = Gy.MultiplyWith4x4(&mh);
	Matrix4x1 Qz = Gz.MultiplyWith4x4(&mh);
	
	//store coefficients
	coef[0] = Qx.M[0];
	coef[1] = Qx.M[1];
	coef[2] = Qx.M[2];
	coef[3] = Qx.M[3];
	coef[4] = Qy.M[0];
	coef[5] = Qy.M[1];
	coef[6] = Qy.M[2];
	coef[7] = Qy.M[3];
	coef[8] = Qz.M[0];
	coef[9] = Qz.M[1];
	coef[10] = Qz.M[2];
	coef[11] =	 Qz.M[3];
}

void CCubicCurve::GenerateBezierCoefficients(Point _p1, Point _p2, Point _p3, Point _p4)
{
	
	type = TM_BEZIER;
	
	Matrix4x1 p1 = ConvertPoint(&_p1);
	Matrix4x1 p2 = ConvertPoint(&_p2);
	Matrix4x1 p3 = ConvertPoint(&_p3);
	Matrix4x1 p4 = ConvertPoint(&_p4);
	
	//build hermite base matrix
	Matrix4x4 mh;
	mh.M[0][0] = -1.0f;
	mh.M[0][1] = 3.0f;
	mh.M[0][2] = -3.0f;
	mh.M[0][3] = 1.0f;
	
	mh.M[1][0] = 3.0f;
	mh.M[1][1] = -6.0f;
	mh.M[1][2] = 3.0f;
	mh.M[1][3] = 0.0f;
	
	mh.M[2][0] = -3.0f;
	mh.M[2][1] = 3.0f;
	mh.M[2][2] = 0.0f;
	mh.M[2][3] = 0.0f;
	
	mh.M[3][0] = 1.0f;
	mh.M[3][1] = 0.0f;
	mh.M[3][2] = 0.0f;
	mh.M[3][3] = 0.0f;
	
	//build geometry matrices
	Gx.M[0] = p1.M[0];
	Gx.M[1] = p2.M[0];
	Gx.M[2] = p3.M[0];
	Gx.M[3] = p4.M[0];
	
	Gy.M[0] = p1.M[1];
	Gy.M[1] = p2.M[1];
	Gy.M[2] = p3.M[1];
	Gy.M[3] = p4.M[1];
	
	Gz.M[0] = p1.M[2];
	Gz.M[1] = p2.M[2];
	Gz.M[2] = p3.M[2];
	Gz.M[3] = p4.M[2];
	
	//computer coefficients
	Matrix4x1 Qx = Gx.MultiplyWith4x4(&mh);
	Matrix4x1 Qy = Gy.MultiplyWith4x4(&mh);
	Matrix4x1 Qz = Gz.MultiplyWith4x4(&mh);
	
	//store coefficients
	coef[0] = Qx.M[0];
	coef[1] = Qx.M[1];
	coef[2] = Qx.M[2];
	coef[3] = Qx.M[3];
	coef[4] = Qy.M[0];
	coef[5] = Qy.M[1];
	coef[6] = Qy.M[2];
	coef[7] = Qy.M[3];
	coef[8] = Qz.M[0];
	coef[9] = Qz.M[1];
	coef[10] = Qz.M[2];
	coef[11] =	 Qz.M[3];
}

/*
Draw the cubic curve
*/
void CCubicCurve::Draw()
{
	//create pointers to coefficients
	float *a = coef;
	float *b = coef+4;
	float *c = coef+8;
	
	//initialise loop vars
	float t = 0.0f;
	float t2, t3;
	float t_step = (1.0f)/(float)steps;
	
	//save old color
	glPushAttrib(GL_CURRENT_BIT);
	glColor4fv(color);
	//main loop, curves drawn as triangle strips
	glBegin(GL_LINE_STRIP);
	for (int j = 0; j <= steps; j++)
	{
		t2 = t*t;
		t3 = t2*t;
		
		float x = a[0]*t3 + a[1]*t2 + a[2]*t  + a[3];
		float y = b[0]*t3 + b[1]*t2 + b[2]*t  + b[3];
		float z = c[0]*t3 + c[1]*t2 + c[2]*t + c[3];
		glVertex3f(x,y,z);
		t+=t_step;
	}
	glEnd();
	
	//restore old color
	glPopAttrib();
}
/*
Draw cotrol points
bezier curve: the points are drawn
hermite curve: the tangent vectors are drawn at the starting points
*/
void CCubicCurve::DrawControlPoints()
{
	glPushAttrib(GL_CURRENT_BIT);
	glPointSize(4);
	glColor4fv(color);
	if (type == TM_BEZIER)
	{
		glBegin(GL_POINTS);
		glVertex3f(Gx.M[0], Gy.M[0], Gz.M[0]);
		glVertex3f(Gx.M[1], Gy.M[1], Gz.M[1]);
		glVertex3f(Gx.M[2], Gy.M[2], Gz.M[2]);
		glVertex3f(Gx.M[3], Gy.M[3], Gz.M[3]);
		glEnd();
		
	}
	else if (type == TM_HERMITE)
	{
		glBegin(GL_LINES);
		glVertex3f(Gx.M[0], Gy.M[0], Gz.M[0]);
		glVertex3f(Gx.M[1], Gy.M[1], Gz.M[1]);
		glVertex3f(Gx.M[2], Gy.M[2], Gz.M[2]);
		glVertex3f(Gx.M[3], Gy.M[3], Gz.M[3]);
		glEnd();
	}
	glPopAttrib();
}

Point CCubicCurve::GetPoint(float t)
{
	Point result;
	float *a = coef;
	float *b = coef+4;
	float *c = coef+8;
	
	float t2, t3;
	t2 = t*t;
	t3 = t2*t;
	
	result.x = a[0]*t3 + a[1]*t2 + a[2]*t  + a[3];
	result.y = b[0]*t3 + b[1]*t2 + b[2]*t  + b[3];
	result.z = c[0]*t3 + c[1]*t2 + c[2]*t + c[3];
	return result;
}



void CCubicCurve::SetGeometryMatrixGx(Matrix4x1 a)
{
	Gx = a;
}

void CCubicCurve::SetGeometryMatrixGy(Matrix4x1 a)
{
	Gy = a;
}

void CCubicCurve::SetGeometryMatrixGz(Matrix4x1 a)
{
	Gz = a;
}

void CCubicCurve::RecalculateBezier()
{
	Matrix4x4 mh;
	mh.M[0][0] = -1.0f;
	mh.M[0][1] = 3.0f;
	mh.M[0][2] = -3.0f;
	mh.M[0][3] = 1.0f;
	
	mh.M[1][0] = 3.0f;
	mh.M[1][1] = -6.0f;
	mh.M[1][2] = 3.0f;
	mh.M[1][3] = 0.0f;
	
	mh.M[2][0] = -3.0f;
	mh.M[2][1] = 3.0f;
	mh.M[2][2] = 0.0f;
	mh.M[2][3] = 0.0f;
	
	mh.M[3][0] = 1.0f;
	mh.M[3][1] = 0.0f;
	mh.M[3][2] = 0.0f;
	mh.M[3][3] = 0.0f;
	
	//computer coefficients
	Matrix4x1 Qx = Gx.MultiplyWith4x4(&mh);
	Matrix4x1 Qy = Gy.MultiplyWith4x4(&mh);
	Matrix4x1 Qz = Gz.MultiplyWith4x4(&mh);
	
	//store coefficients
	coef[0] = Qx.M[0];
	coef[1] = Qx.M[1];
	coef[2] = Qx.M[2];
	coef[3] = Qx.M[3];
	coef[4] = Qy.M[0];
	coef[5] = Qy.M[1];
	coef[6] = Qy.M[2];
	coef[7] = Qy.M[3];
	coef[8] = Qz.M[0];
	coef[9] = Qz.M[1];
	coef[10] = Qz.M[2];
	coef[11] =	 Qz.M[3];
}

void CCubicCurve::SetColor(float r, float g, float b, float a)
{
	color[0] = r;
	color[1] = g;
	color[2] = b;
	color[3] = a;
}

CCubicCurve::CCubicCurve() : steps(10), type(TM_NOTHING)
{
	color[0] = color[1] = color[2] = 0.5f;
	color[3] = 0.5f;
}

CCubicCurve::~CCubicCurve()
{
	
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
/////////////////// PATCH ////////////////////////////////////
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

CCubicPatch::CCubicPatch() : s_evals(10), curve_evals(10), recompute(true), calc_normals(false)
{  			
}

CCubicPatch::~CCubicPatch()
{
	printf("\nDeallocating memory for patch\n");
	delete[] tris;
	delete[] normals;
}

/* Creates patch, assumes that p is an array of 16 points */
void CCubicPatch::GeneratePatch(Point *p)
{
	printf("\nGenerates patch");
	//one way
	Curve[0].GenerateBezierCoefficients(p[0], p[1], p[2], p[3]);
	Curve[1].GenerateBezierCoefficients(p[4], p[5], p[6], p[7]);
	Curve[2].GenerateBezierCoefficients(p[8], p[9], p[10], p[11]);
	Curve[3].GenerateBezierCoefficients(p[12], p[13], p[14], p[15]);
	
	//the other way round
	Curve[4].GenerateBezierCoefficients(p[0], p[4], p[8], p[12]);
	Curve[5].GenerateBezierCoefficients(p[1], p[5], p[9], p[13]);
	Curve[6].GenerateBezierCoefficients(p[2], p[6], p[10], p[14]);
	Curve[7].GenerateBezierCoefficients(p[3], p[7], p[11], p[15]);
	
	//save points
	for (int i = 0; i < 16; i++)
	{
		ControlPoints[i] = p[i];
	}
}

void CCubicPatch::ChangeEval(int eval, int factor)
{
	if (eval == TM_S_EVALS)
	{
		s_evals += factor;
		if (s_evals < 2)
			s_evals = 2;
		recompute = true;
	} else if (eval == TM_CURVE_EVALS)
	{
		curve_evals += factor;
		if (curve_evals < 2)
			curve_evals = 2;
		recompute = true;		
	}
}

inline void CrossProduct(Point *v, Point *v1, Point *result)
{
	result->x = v->y * v1->z - v1->y * v->z;
	result->y = v->z * v1->x - v1->z * v->x;
	result->z = v->x * v1->y - v1->x * v->y;
	float length = (float)sqrt(result->x*result->x+result->y*result->y+result->z*result->z);
	result->x /= length;
	result->y /= length;
	result->z /= length;
}

void CCubicPatch::ToggleNormals()
{
	calc_normals = !calc_normals;
	recompute = true;
}

void CCubicPatch::Draw()
{
	glPointSize(1);
	
	if (recompute)
	{
		Point* s_points;
		s_points = new Point[s_evals*curve_evals]; ///points from different s values
		Point *dir1 = s_points;
		float step_s = (1.0f)/(float(s_evals-1));
		float step = (1.0f)/((float)(curve_evals-1));
		
		Point p[4];	//storage for control points
		CCubicCurve c;	//curve created in loop
		float s = 0.0f;
		float s_c = 0.0f; //eval for each curve
		for (unsigned int index_s = 0; index_s < s_evals; index_s++, s += step_s)
		{
			c.SetColor(0,1.0f,0,0.5f);
			//get the four points
			p[0] = Curve[4].GetPoint(s);
			p[1] = Curve[5].GetPoint(s);
			p[2] = Curve[6].GetPoint(s);
			p[3] = Curve[7].GetPoint(s);
			c.GenerateBezierCoefficients(p[0], p[1], p[2], p[3]); //generate new curve
			s_c = 0.0f;
			for (unsigned int c1 = 0; c1 < curve_evals; c1++, s_c+=step)
			{
				(*dir1) = c.GetPoint(s_c);
				dir1++;
			}
		}
		
		printf("\nComputes/Recomputes patch triangles");
		unsigned int number_of_triangles = 2*(s_evals-1)*(curve_evals-1);
		
		dir1 = s_points;
		delete[] tris;
		tris = new  Triangle[number_of_triangles];//number of triangles
		
		//create triangles
		int tris_count = 0;
		for (unsigned int j = 0; j < s_evals-1; j++)
		{
			for (unsigned int i = 0; i < curve_evals-1; i++)
			{
				tris[tris_count].verts[0] = dir1[j*curve_evals + i];
				tris[tris_count].verts[1] = dir1[(j+1)*curve_evals + i];
				tris[tris_count].verts[2] = dir1[(j+1)*curve_evals + i+1];
				tris_count++;
				tris[tris_count].verts[2] = dir1[j*curve_evals + i];
				tris[tris_count].verts[1] = dir1[j*curve_evals + i+1];
				tris[tris_count].verts[0] = dir1[(j+1)*curve_evals + i+1];
				tris_count++;
				
			}
		}
					
		if (calc_normals)
		{
			delete[] normals;
			normals = new Point[number_of_triangles];	//allocate space for normals
			
			//calculate normals (one normal for each triangle.
			for (unsigned int index = 0; index < number_of_triangles; index++)
			{
				//create two vectors
				Point v,v1;
				v.x = tris[index].verts[1].x - tris[index].verts[0].x;
				v.y = tris[index].verts[1].y - tris[index].verts[0].y;
				v.z = tris[index].verts[1].z - tris[index].verts[0].z;
				
				v1.x = tris[index].verts[2].x - tris[index].verts[0].x;
				v1.y = tris[index].verts[2].y - tris[index].verts[0].y;
				v1.z = tris[index].verts[2].z - tris[index].verts[0].z;
				
				//calculate and store normal
				CrossProduct(&v, &v1, &normals[index]);
			}
			
		}
		printf("\nNumber of triangles: %d", number_of_triangles);
		delete[] s_points;
		recompute = false;
	}
	
	
	int num_tris = 2*(s_evals-1)*(curve_evals-1);
	//	glEnable(GL_CULL_FACE);
	glBegin(GL_TRIANGLES);
	bool red = true;
	for (int tris_count = 0; tris_count < num_tris; tris_count++)
	{
		if (calc_normals)	//if normals should be specified
			glNormal3f(normals[tris_count].x,normals[tris_count].y,normals[tris_count].z);
		glVertex3f(tris[tris_count].verts[0].x, tris[tris_count].verts[0].y, tris[tris_count].verts[0].z);
		glVertex3f(tris[tris_count].verts[1].x, tris[tris_count].verts[1].y, tris[tris_count].verts[1].z);
		glVertex3f(tris[tris_count].verts[2].x, tris[tris_count].verts[2].y, tris[tris_count].verts[2].z);
	}
	glEnd();
	
	
}
#endif CURVES_CPP__