void CalcEyeFrustumInObjSpace(GEVector3D* p, GECamera* cam){        float t = cam->GetZNear() * tan(geDegToRad(cam->GetFOV()) * 0.5f);	float b = -t;	float r = cam->GetAspect() * t;	float l = -r;	float n = cam->GetZNear();	float f = cam->GetZFar();	float fn = f / n;	float lFar = fn * l;	float rFar = fn * r;	float bFar = fn * b;	float tFar = fn * t;	p[0].Set(l, b, -n);	p[1].Set(r, b, -n);	p[2].Set(r, t, -n);	p[3].Set(l, t, -n);	p[4].Set(lFar, bFar, -f);	p[5].Set(rFar, bFar, -f);	p[6].Set(rFar, tFar, -f);	p[7].Set(lFar, tFar, -f);}void CalcEyeFrustumInL(GEVector3D* eyeFr, GEVector3D* eyeFrInL, GECamera* eye){	GEMatrix4x4 lookAt, proj, eyeProj, camLookAt_inv;	camLookAt_inv.LookAt(eye->Origin, eye->Dir, eye->UpVector); 	camLookAt_inv.Invert(camLookAt_inv);	Light->LightView.LookAt(Light->GetCamera()->Origin, Light->GetCamera()->Dir, Light->GetCamera()->UpVector);	Light->LightProj.Perspective(Light->GetCamera()->GetFOV(), 1.0f, Light->GetCamera()->GetZNear(), Light->GetCamera()->GetZFar());	eyeProj.Multiply(Light->LightProj, Light->LightView);	eyeProj.Multiply(eyeProj, camLookAt_inv);	//////////////////////////////////////////////////////////////////////////	// eyeFrustumInLightSpace = P_L * C_L * C_E^-1 * eyeFrustumInObjSpace	//	// Chapter 3, last par. : The eight corner vertices of E are obtained from 	// those corner vertices of the eye's frustum in the object space multiplied 	// by PL .CL .CE^-1 where CE^-1 is the inverse camera matrix of the eye. We treat 	// E as a flattened 2D object on the front face of the light's unit cube.	//	IsEyeFrustumInLightFrustum = true;	for(int i=0;i<8;i++)	{		eyeProj.TransformHomogenPoint(eyeFr, eyeFrInL);		eyeFrInL.z = 1.0f;		if(fabsf(eyeFrInL.x) > 1.0f || fabsf(eyeFrInL.y) > 1.0f)			IsEyeFrustumInLightFrustum = false;	}	RenderFrustum(eyeFrInL); 	glBegin(GL_POINTS);	glColor3f(1.0f, 0.0f, 0.0f);	for(i=0;i<4;i++) 		glVertex3fv(&eyeFrInL.x);	glColor3f(0.0f, 1.0f, 0.0f);	for(i=4;i<8;i++) 		glVertex3fv(&eyeFrInL.x); 	glEnd();}void CalcCenterLine(GEVector3D* eyeFrInL, GEVector3D* l){	GEVector3D v1 = eyeFrInL[1] - eyeFrInL[0];	GEVector3D v2 = eyeFrInL[2] - eyeFrInL[3];	v1 = eyeFrInL[0] + v1 * 0.5f;	v2 = eyeFrInL[3] + v2 * 0.5f;		l[0] = v1 + (v2 - v1) * 0.5f;	v1 = eyeFrInL[5] - eyeFrInL[4];	v2 = eyeFrInL[6] - eyeFrInL[7];	v1 = eyeFrInL[4] + v1 * 0.5f;	v2 = eyeFrInL[7] + v2 * 0.5f;	l[1] = v1 + (v2 - v1) * 0.5f;	glBegin(GL_LINES);	glVertex3fv(&l[0].x);	glVertex3fv(&l[1].x);	glEnd();}void CalcTrapezoid(GEVector3D* hull, int numHullPoints, GEVector3D* line, GEVector3D* trapezoid){	//////////////////////////////////////////////////////////////////////////	// We need an axis orthogonal to center line.	// Because all the points of the hull are on the z = 1.0f plane, we can take	// the cross product of the center line direction with the z axis vector	// as the perpendicular axis.	//	// FIX : What happens if line[1] = line[0]??? Center line has zero length, 	// and we can't find a valid perpendicular axis! This may happen when we	// are near the dualling frusta case.	//	GEVector3D perpAxis = line[1] - line[0];	perpAxis.Normalize();	if(perpAxis.Length() < 1e-5f)		return;	//////////////////////////////////////////////////////////////////////////	// This is because we always want the perpAxis to point to +x world axis.	// The normal way to do it is to change the vectors in the cross product, 	// but this makes things "symmetrical"!!!	//	if(perpAxis.x < 0.0f)		perpAxis = perpAxis | GEVector3D(0.0f, 0.0f, 1.0f);	else		perpAxis = perpAxis | GEVector3D(0.0f, 0.0f, -1.0f);	//////////////////////////////////////////////////////////////////////////	// If we transform all the points of the hull, so that the above calculated	// perpendicular axis is the x - axis, and the center line is the y - axis,	// we can calculate an axis-aligned bounding box of hull's points. With this	// aabb we can calculate how much we must travel along the center line in case to	// draw the top and the base lines that encloses the whole convex hull.	// We also translate the the hull so the center line starts at (0,0,0). 	//	GEMatrix4x4 rotz, trans;	rotz.MakeRotationMatrix('z', geRadToDeg(-acos(perpAxis.x)));	trans.MakeTranslationMatrix(-line[0].x, -line[0].y, -line[0].z);	rotz.Multiply(rotz, trans);		GEVector3D minv(2.0f, 2.0f, 1.0f), maxv(-2.0f, -2.0f, 1.0f);	GEVector3D transHull[6];	for(int i=0;i<numHullPoints;i++)	{		rotz.TransformPoint(hull, transHull);		//////////////////////////////////////////////////////////////////////////		// We check only y component, because z is always 1.0f, and we are interested,		// in the y-direction only. This is, we only want the top and the bottom of the		// aabbox.		//		if(transHull.y > maxv.y)		maxv.y = transHull.y;		if(transHull.y < minv.y)		minv.y = transHull.y;	}		//////////////////////////////////////////////////////////////////////////	// We do this (the min/max thing for top_y and base_y) because, we always want to	// have the top line touch the near plane, and the base line touch the far plane.	//	float topy = min(fabs(minv.y), fabs(maxv.y));	float basey = max(fabs(minv.y), fabs(maxv.y));	GEVector3D clineDir = line[1] - line[0];	float cLineLen = clineDir.Length();	clineDir.Normalize();	GEVector3D topLine[2];	GEVector3D baseLine[2];	GEVector3D temp;	//////////////////////////////////////////////////////////////////////////	// We move along the center line, starting from the point on the near plane,	// by a distance of -miny to reach the outer point of the hull. Then with the	// perpendicular axis, we form the top line.	//	temp = line[0] - clineDir * topy;	topLine[0] = temp + perpAxis;	topLine[1] = temp - perpAxis;	//////////////////////////////////////////////////////////////////////////	// The same as above, but now we must move +maxy to reach the cross point of 	// the center line and the base line.	//	temp = line[0] + clineDir * basey;	baseLine[0] = temp + perpAxis;	baseLine[1] = temp - perpAxis;	//////////////////////////////////////////////////////////////////////////	// Calculate distance of point q from the top line ('&#951;' in the paper)	// Chapter 6.2.	//	float lineFactor = 0.8f;				// 80% line	float ksi = 1.0f - 2.0f * lineFactor;	// &#958; in the paper	float lamda = cLineLen + fabs(basey) + fabs(topy);	// distance between the top and the base line	float delta = 1.0f - ksi;// + topy;	float n = (lamda * delta + lamda * delta * ksi) / (-lamda - 2.0f * delta - lamda * ksi);		//////////////////////////////////////////////////////////////////////////	// Calculate q	//	temp = line[0] - clineDir * topy;	// temp is the cross between center line and top line.	GEVector3D q = temp + clineDir * n;	//////////////////////////////////////////////////////////////////////////	// Find the two points of the hull which combined with q give us the side 	// edges of the trapezoid. These points, must be points of the hull.	//	//////////////////////////////////////////////////////////////////////////	// For every point on the hull, form a line from q to the point. We know that	// all points lie on z = +1 plane. So, from this line, we form a plane, and we 	// check every other point on the hull against it. In case to have a winner, 	// all other points must lie on the same side of the plane.	//	int index[2] = {-1, -1};	int curIndex = 0;	for(i=0;i<numHullPoints;i++)	{		GEPlane plane(hull, (q - hull) | GEVector3D(0.0f, 0.0f, 1.0f));				int numPos = 0, numNeg = 0;		for(int j=0;j<numHullPoints;j++)		{			if(i == j)				continue;			float dist = plane.ClassifyPointf(hull[j]);			if(dist >= 0.0f)				numPos++;			else				numNeg++;		}		if((numPos > 0 && numNeg == 0) || (numPos == 0 && numNeg > 0))		{			index[curIndex++] = i;			if(curIndex > 2)				break;		}	}	if(curIndex != 2)		MessageBox(NULL, "Severe Error : Couldn't find to extreme point for calculating the trapezoid.", "ERROR", MB_OK);	GEVector3D pMin, pMax;	if(transHull[index[0]].x > transHull[index[1]].x)	{		pMax = hull[index[0]];		pMin = hull[index[1]];	}	else	{		pMax = hull[index[1]];		pMin = hull[index[0]];	}		geIntersectLines2D(trapezoid[0], baseLine[0], baseLine[1], q, pMin);	geIntersectLines2D(trapezoid[1], baseLine[0], baseLine[1], q, pMax);	geIntersectLines2D(trapezoid[2], topLine[0], topLine[1], q, pMax);	geIntersectLines2D(trapezoid[3], topLine[0], topLine[1], q, pMin);	//////////////////////////////////////////////////////////////////////////	// Check trapezoid points order.	// 	GEVector3D dir1 = trapezoid[1] - trapezoid[0];	GEVector3D dir2 = trapezoid[2] - trapezoid[0];	dir1.Normalize();	dir2.Normalize();	dir1 = dir1 | dir2;	dir1.Normalize();	if(!(dir1 == GEVector3D(0.0f, 0.0f, 1.0f)))	{		dir1 = trapezoid[0];		trapezoid[0] = trapezoid[1];		trapezoid[1] = dir1;		dir1 = trapezoid[3];		trapezoid[3] = trapezoid[2];		trapezoid[2] = dir1;	}	//////////////////////////////////////////////////////////////////////////	// Debug visualization	//	glLineWidth(3.0f);	glColor3f(0.0f, 0.0f, 0.0f);	glBegin(GL_LINE_LOOP);	glColor3f(1.0f, 0.0f, 0.0f);	glVertex3fv(&trapezoid[0].x);	glColor3f(0.0f, 1.0f, 0.0f);	glVertex3fv(&trapezoid[1].x);	glColor3f(0.0f, 0.0f, 1.0f);	glVertex3fv(&trapezoid[2].x);	glColor3f(1.0f, 1.0f, 1.0f);	glVertex3fv(&trapezoid[3].x);	glEnd();	glColor3f(1.0f, 1.0f, 1.0f);	glLineWidth(1.0f);	glBegin(GL_LINES);	glColor3f(1.0f, 0.0f, 0.0f);	glVertex3fv(&topLine[0].x);	glVertex3fv(&topLine[1].x);	glColor3f(0.0f, 1.0f, 0.0f);	glVertex3fv(&baseLine[0].x);	glVertex3fv(&baseLine[1].x);	glEnd();	glPointSize(3.0f);	glColor3f(0.5f, 1.0f, 0.5f);	glBegin(GL_POINTS);	glVertex3fv(&q.x);	glEnd();	glColor3f(1.0f, 1.0f, 1.0f);	glPointSize(1.0f);}

GEVector3D eyeFrustum [8];GEVector3D eyeFrustumInL [8];GEVector3D line [2];GEVector3D hull [6];GEVector3D trapezoid [4];GEMatrix4x4* GELight::N_T;CalcEyeFrustumInObjSpace(&eyeFrustum[0], Camera);CalcEyeFrustumInL(&eyeFrustum[0], &eyeFrustumInL[0], Camera);CalcCenterLine(&eyeFrustumInL[0], &line[0]);numHullPoints = CalcConvexHull2D(&eyeFrustumInL[0], &hull[0]);RenderConvexHull(&hull[0], numHullPoints);CalcTrapezoid(&hull[0], numHullPoints, &line[0], &trapezoid[0]);CalcN_T(&trapezoid[0], &Light->N_T, &eyeFrustum[0]);