1 #include "SDL/SDL_opengl.h"
2 #include <math.h>
3 #include "raytracing.h"
4 
5 #ifndef pi
6 #define pi 3.14159265358979
7 #endif
8 
9 
10 const double radian = (double) (pi / 180);
11 const double epsilon = (double) (0.05);
12 
13 
14 typedef double SCALAR;
15 
16 class VECTOR
17 {
18 	public:
19 		VECTOR(double sx = 0, double sy = 0, double sz = 0);
20 		VECTOR(const VECTOR& Vector);
21 		~VECTOR();
22 
23 		double GetMagnitude();
24 		GLvoid Normalize();
25 		GLvoid Reset();
Set(double sx,double sy,double sz)26 		GLvoid Set(double sx, double sy, double sz) {x = sx, y = sy, z = sz;}
27 		GLvoid CrossVector(VECTOR vect);
28 		double DotProduct(VECTOR vect);
29 
30         //equal within an error e
nearlyEquals(const VECTOR & v,const SCALAR e) const31 		const bool nearlyEquals( const VECTOR& v, const SCALAR e ) const
32 		{
33 			return fabs(x-v.x)<e && fabs(y-v.y)<e && fabs(z-v.z)<e;
34 		}
35 
36         //cross product
cross(const VECTOR & v) const37 		const VECTOR cross( const VECTOR& v ) const
38 		{
39             //Davis, Snider, "Introduction to Vector Analysis", p. 44
40 			return VECTOR( y*v.z - z*v.y, z*v.x - x*v.z, x*v.y - y*v.x );
41 		}
42 
43         //scalar dot product
dot(const VECTOR & v) const44 		const SCALAR dot( const VECTOR& v ) const
45 		{
46 			return x*v.x + y*v.y + z*v.z;
47 		}
48 
49         //length
length() const50 		const SCALAR length() const
51 		{
52 			return (SCALAR)sqrt( (double)this->dot(*this) );
53 		}
54 
55         //unit vector
unit() const56 		const VECTOR unit() const
57 		{
58 			return (*this) / length();
59 		}
60 
61         //make this a unit vector
normalize()62 		void normalize()
63 		{
64 			(*this) /= length();
65 		}
66 
67         //Members
68 		double x;
69 		double y;
70 		double z;
71 
72         //index a component
73         //NOTE: returning a reference allows
74         //you to assign the indexed element
operator [](const long i)75 		SCALAR& operator [] ( const long i )
76 		{
77 			return *((&x) + i);
78 		}
79 
80         //compare
operator ==(const VECTOR & v) const81 		const bool operator == ( const VECTOR& v ) const
82 		{
83 			return (v.x==x && v.y==y && v.z==z);
84 		}
85 
operator !=(const VECTOR & v) const86 		const bool operator != ( const VECTOR& v ) const
87 		{
88 			return !(v == *this);
89 		}
90 
91         //negate
operator -() const92 		const VECTOR operator - () const
93 		{
94 			return VECTOR( -x, -y, -z );
95 		}
96 
97         //assign
operator =(const VECTOR & v)98 		const VECTOR& operator = ( const VECTOR& v )
99 		{
100 			x = v.x;
101 			y = v.y;
102 			z = v.z;
103 			return *this;
104 		}
105 
106         //increment
operator +=(const VECTOR & v)107 		const VECTOR& operator += ( const VECTOR& v )
108 		{
109 			x+=v.x;
110 			y+=v.y;
111 			z+=v.z;
112 			return *this;
113 		}
114 
115         //decrement
operator -=(const VECTOR & v)116 		const VECTOR& operator -= ( const VECTOR& v )
117 		{
118 			x-=v.x;
119 			y-=v.y;
120 			z-=v.z;
121 			return *this;
122 		}
123 
124         //self-multiply
operator *=(const SCALAR & s)125 		const VECTOR& operator *= ( const SCALAR& s )
126 		{
127 			x*=s;
128 			y*=s;
129 			z*=s;
130 			return *this;
131 		}
132 
133         //self-divide
operator /=(const SCALAR & s)134 		const VECTOR& operator /= ( const SCALAR& s )
135 		{
136 			const SCALAR r = 1 / s;
137 			x *= r;
138 			y *= r;
139 			z *= r;
140 			return *this;
141 		}
142 
143         //add
operator +(const VECTOR & v) const144 		const VECTOR operator + ( const VECTOR& v ) const
145 		{
146 			return VECTOR(x + v.x, y + v.y, z + v.z);
147 		}
148 
149         //subtract
operator -(const VECTOR & v) const150 		const VECTOR operator - ( const VECTOR& v ) const
151 		{
152 			return VECTOR(x - v.x, y - v.y, z - v.z);
153 		}
154 
155         //post-multiply by a scalar
operator *(const SCALAR & s) const156 		const VECTOR operator * ( const SCALAR& s ) const
157 		{
158 			return VECTOR( x*s, y*s, z*s );
159 		}
160 
161         //pre-multiply by a scalar
operator *(const SCALAR & s,const VECTOR & v)162 		friend inline const VECTOR operator * ( const SCALAR& s, const VECTOR& v )
163 		{
164 			return v * s;
165 		}
166 
167         //divide
operator /(SCALAR s) const168 		const VECTOR operator / (SCALAR s) const
169 		{
170 			s = 1/s;
171 			return VECTOR( s*x, s*y, s*z );
172 		}
173 };
174 
175 
176 
177 class VERTEX
178 {
179 	public:
180 		VERTEX(double x = 0, double y = 0, double z = 0, double nx = 0, double ny = 0, double nz = 0);
181 		~VERTEX();
182 
183 		VECTOR coords;
184 		VECTOR normal;
185 		double u;
186 		double v;
187 
188         // operator overloading
operator =(const VERTEX & Vertex)189 		VERTEX& operator = (const VERTEX& Vertex)
190 		{
191 			coords = Vertex.coords;
192 			normal = Vertex.normal;
193 			u = Vertex.u;
194 			v = Vertex.v;
195 			return *this;
196 		}
197 };
198 
199 
200 class QUAD
201 {
202 	public:
203 		QUAD();
204 		~QUAD();
205 
206 		VECTOR GetNormal();
207 		GLvoid SetNormal();
208 
209 		VERTEX Vertex[4];
210 		int numVertices;
211 };
212 
213 
214 class QUAT
215 {
216 	public:
217 		QUAT(double sx = 0, double sy = 0, double sz = 0, double sw = 1);
218 		~QUAT();
219 
220 		GLvoid Reset();
221 		GLvoid CopyQuat(QUAT q);
Set(double sx,double sy,double sz,double sw)222 		GLvoid Set(double sx, double sy, double sz, double sw) {x = sx, y = sy, z = sz, w = sw;}
223 		GLvoid AxisAngleToQuat(VECTOR axis, double theta);
224 		GLvoid EulerToQuat(double pitch, double yaw, double roll);
225 		GLvoid NormaliseQuat();
226 		double MagnitudeQuat();
227 		GLvoid MultQuat(QUAT q);
228 
229 		double x;
230 		double y;
231 		double z;
232 		double w;
233 };
234 
235 
236 
237 
238 class MATRIX
239 {
240 	public:
241 		MATRIX();
242 		~MATRIX();
243 
244 		GLvoid LoadIdentity();
245 		GLvoid CopyMatrix(double m[16]);
246 		GLvoid MultMatrix(double m[16]);
247 		GLvoid MatrixInverse();
248 		GLvoid MatrixFromAxisAngle(VECTOR axis, double theta);
249 		GLvoid QuatToMatrix(QUAT quat);
250 
251 		double Element[16];
252 };
253 
254 
255 
256 
257 
258 class OBJECT
259 {
260 	public:
261 		OBJECT();
262 		~OBJECT();
263 
264 		GLvoid Reset();
265 		GLvoid Rotate();
266 		GLvoid Draw();
267 		GLvoid UpdatePosition();
268 		GLvoid UpdatePosition(double x, double y, double z);
269 		GLvoid MoveX();
270 		GLvoid MoveY();
271 		GLvoid MoveZ();
272 		GLvoid MoveX(double x);
273 		GLvoid MoveY(double y);
274 		GLvoid MoveZ(double z);
275 		VECTOR GetXUnit();
276 		VECTOR GetYUnit();
277 		VECTOR GetZUnit();
278 
279 		QUAT Orientation;
280 		VECTOR Position;
281 		double Delta_x;   //Rotation deltas
282 		double Delta_y;
283 		double Delta_z;
284 		double Movement_x;    //Movement displacements
285 		double Movement_y;
286 		double Movement_z;
287 		double Multiplier;
288 };
289 
290 
291 class CAMERA : public OBJECT
292 {
293 	public:
294 		CAMERA();
295 		~CAMERA();
296 
297 		GLvoid Reset();
298 		GLvoid Update();
299 		GLvoid Apply();
300 
301 		MATRIX Matrix;
302 };
303 
304 
305 
306 
307 
308 
309 
310 
VECTOR(double sx,double sy,double sz)311 VECTOR::VECTOR(double sx, double sy, double sz)
312 :
313 	x(sx),
314 	y(sy),
315 	z(sz)
316 {
317 }
318 
VECTOR(const VECTOR & Vector)319 VECTOR::VECTOR(const VECTOR& Vector)
320 {
321 	x = Vector.x;
322 	y = Vector.y;
323 	z = Vector.z;
324 }
325 
~VECTOR()326 VECTOR::~VECTOR()
327 {
328 }
329 
Reset()330 GLvoid VECTOR::Reset()
331 {
332 	x = 0;
333 	y = 0;
334 	z = 0;
335 }
336 
DotProduct(VECTOR vect)337 double VECTOR::DotProduct(VECTOR vect)
338 {
339   	double dot;
340   	dot = vect.x * x + vect.y * y + vect.z * z;
341   	return dot;
342 }
343 
CrossVector(VECTOR vect)344 GLvoid VECTOR::CrossVector(VECTOR vect)
345 {
346   	VECTOR temp = *this;
347   	x = vect.y * temp.z - vect.z * temp.y;
348   	y = vect.z * temp.x - vect.x * temp.z;
349   	z = vect.x * temp.y - vect.y * temp.x;
350 }
351 
GetMagnitude()352 double VECTOR::GetMagnitude()
353 {
354 	double presque_rien=0.0000000000001f;
355     double magnitude = sqrt(x * x + y * y + z * z);
356     if (magnitude > presque_rien)
357         return magnitude;
358     else
359         return  presque_rien;
360 }
361 
Normalize()362 GLvoid VECTOR::Normalize()
363 {
364     double magnitude = this->GetMagnitude();
365     x /= magnitude;
366     y /= magnitude;
367     z /= magnitude;
368 }
369 
370 
371 
VERTEX(double sx,double sy,double sz,double snx,double sny,double snz)372 VERTEX::VERTEX(double sx, double sy, double sz, double snx, double sny, double snz)
373 {
374 	coords.x = sx;
375 	coords.y = sy;
376 	coords.z = sz;
377 	normal.x = snx;
378 	normal.y = sny;
379 	normal.z = snz;
380 }
381 
~VERTEX()382 VERTEX::~VERTEX()
383 {
384 }
385 
386 
QUAD()387 QUAD::QUAD()
388 {
389 	numVertices=4;
390 }
391 
~QUAD()392 QUAD::~QUAD()
393 {
394 }
395 
GetNormal()396 VECTOR QUAD::GetNormal()
397 {
398         VECTOR temp;
399     	double ux;
400     	double uy;
401     	double uz;
402     	double vx;
403     	double vy;
404     	double vz;
405       	ux = Vertex[1].coords.x - Vertex[0].coords.x;
406       	uy = Vertex[1].coords.y - Vertex[0].coords.y;
407       	uz = Vertex[1].coords.z - Vertex[0].coords.z;
408       	vx = Vertex[2].coords.x - Vertex[0].coords.x;
409       	vy = Vertex[2].coords.y - Vertex[0].coords.y;
410       	vz = Vertex[2].coords.z - Vertex[0].coords.z;
411       	temp.x = (uy*vz)-(vy*uz);
412       	temp.y = (uz*vx)-(vz*ux);
413       	temp.z = (ux*vy)-(vx*uy);
414         return temp;
415 }
416 
SetNormal()417 GLvoid QUAD::SetNormal()
418 {
419     	double ux;
420     	double uy;
421     	double uz;
422     	double vx;
423     	double vy;
424     	double vz;
425       	ux = Vertex[1].coords.x - Vertex[0].coords.x;
426       	uy = Vertex[1].coords.y - Vertex[0].coords.y;
427       	uz = Vertex[1].coords.z - Vertex[0].coords.z;
428       	vx = Vertex[2].coords.x - Vertex[0].coords.x;
429       	vy = Vertex[2].coords.y - Vertex[0].coords.y;
430       	vz = Vertex[2].coords.z - Vertex[0].coords.z;
431       	Vertex[0].normal.x = (uy*vz)-(vy*uz);
432       	Vertex[0].normal.y = (uz*vx)-(vz*ux);
433       	Vertex[0].normal.z = (ux*vy)-(vx*uy);
434         Vertex[1].normal = Vertex[0].normal;
435       	Vertex[2].normal = Vertex[0].normal;
436 		Vertex[3].normal = Vertex[0].normal;
437 }
438 
439 
440 
MATRIX()441 MATRIX::MATRIX()
442 {
443 	LoadIdentity();
444 }
445 
~MATRIX()446 MATRIX::~MATRIX()
447 {
448 }
449 
LoadIdentity()450 GLvoid MATRIX::LoadIdentity()
451 {
452 	Element[0]=1.0f;
453     	Element[1]=0.0f;
454     	Element[2]=0.0f;
455     	Element[3]=0.0f;
456 
457     	Element[4]=0.0f;
458     	Element[5]=1.0f;
459     	Element[6]=0.0f;
460     	Element[7]=0.0f;
461 
462     	Element[8]=0.0f;
463     	Element[9]=0.0f;
464     	Element[10]=1.0f;
465     	Element[11]=0.0f;
466 
467     	Element[12]=0.0f;
468     	Element[13]=0.0f;
469     	Element[14]=0.0f;
470     	Element[15]=1.0f;
471 }
472 
CopyMatrix(double m[16])473 GLvoid MATRIX::CopyMatrix(double m[16])
474 {
475     	Element[0 ] = m[0 ];
476     	Element[1 ] = m[1 ];
477     	Element[2 ] = m[2 ];
478     	Element[3 ] = m[3 ];
479     	Element[4 ] = m[4 ];
480     	Element[5 ] = m[5 ];
481     	Element[6 ] = m[6 ];
482     	Element[7 ] = m[7 ];
483     	Element[8 ] = m[8 ];
484     	Element[9 ] = m[9 ];
485     	Element[10] = m[10];
486     	Element[11] = m[11];
487     	Element[12] = m[12];
488     	Element[13] = m[13];
489     	Element[14] = m[14];
490     	Element[15] = m[15];
491 }
492 
MultMatrix(double m[])493 GLvoid MATRIX::MultMatrix(double m[])
494 {
495     	MATRIX temp;
496 
497     	temp.CopyMatrix(this->Element);
498 
499   	Element[0] = temp.Element[0 ]*m[0 ]
500         	+ temp.Element[4 ]*m[1 ]
501         	+ temp.Element[8 ]*m[2 ]
502         	+ temp.Element[12]*m[3 ];
503 
504   	Element[1] = temp.Element[1 ]*m[0 ]
505         	+ temp.Element[5 ]*m[1 ]
506         	+ temp.Element[9 ]*m[2 ]
507         	+ temp.Element[13]*m[3 ];
508 
509   	Element[2] = temp.Element[2 ]*m[0 ]
510        		+ temp.Element[6 ]*m[1 ]
511        		+ temp.Element[10]*m[2 ]
512        		+ temp.Element[14]*m[3 ];
513 
514 	Element[3] = temp.Element[3 ]*m[0 ]
515        		+ temp.Element[7 ]*m[1 ]
516        		+ temp.Element[11]*m[2 ]
517        		+ temp.Element[15]*m[3 ];
518 
519   	Element[4] = temp.Element[0 ]*m[4 ]
520        		+ temp.Element[4 ]*m[5 ]
521        		+ temp.Element[8 ]*m[6 ]
522        		+ temp.Element[12]*m[7 ];
523 
524   	Element[5] = temp.Element[1 ]*m[4 ]
525        		+ temp.Element[5 ]*m[5 ]
526        		+ temp.Element[9 ]*m[6 ]
527        		+ temp.Element[13]*m[7 ];
528 
529   	Element[6] = temp.Element[2 ]*m[4 ]
530        		+ temp.Element[6 ]*m[5 ]
531        		+ temp.Element[10]*m[6 ]
532        		+ temp.Element[14]*m[7 ];
533 
534   	Element[7] = temp.Element[3 ]*m[4 ]
535        		+ temp.Element[7 ]*m[5 ]
536        		+ temp.Element[11]*m[6 ]
537        		+ temp.Element[15]*m[7 ];
538 
539   	Element[8] = temp.Element[0 ]*m[8 ]
540        		+ temp.Element[4 ]*m[9 ]
541        		+ temp.Element[8 ]*m[10]
542        		+ temp.Element[12]*m[11];
543 
544   	Element[9] = temp.Element[1 ]*m[8 ]
545        		+ temp.Element[5 ]*m[9 ]
546        		+ temp.Element[9 ]*m[10]
547        		+ temp.Element[13]*m[11];
548 
549   	Element[10]= temp.Element[2 ]*m[8 ]
550        		+ temp.Element[6 ]*m[9 ]
551        		+ temp.Element[10]*m[10]
552        		+ temp.Element[14]*m[11];
553 
554   	Element[11]= temp.Element[3 ]*m[8 ]
555        		+ temp.Element[7 ]*m[9 ]
556        		+ temp.Element[11]*m[10]
557        		+ temp.Element[15]*m[11];
558 
559   	Element[12]= temp.Element[0 ]*m[12]
560        		+ temp.Element[4 ]*m[13]
561        		+ temp.Element[8 ]*m[14]
562        		+ temp.Element[12]*m[15];
563 
564   	Element[13]= temp.Element[1 ]*m[12]
565        		+ temp.Element[5 ]*m[13]
566        		+ temp.Element[9 ]*m[14]
567        		+ temp.Element[13]*m[15];
568 
569   	Element[14]= temp.Element[2 ]*m[12]
570        		+ temp.Element[6 ]*m[13]
571        		+ temp.Element[10]*m[14]
572        		+ temp.Element[14]*m[15];
573 
574   	Element[15]= temp.Element[3 ]*m[12]
575        		+ temp.Element[7 ]*m[13]
576        		+ temp.Element[11]*m[14]
577        		+ temp.Element[15]*m[15];
578 }
579 
MatrixInverse()580 GLvoid MATRIX::MatrixInverse()
581  {
582   	MATRIX temp;
583 
584   	temp.CopyMatrix(this->Element);
585 
586   	Element[0 ] = temp.Element[0 ];
587   	Element[1 ] = temp.Element[4 ];
588   	Element[2 ] = temp.Element[8 ];
589 
590   	Element[4 ] = temp.Element[1 ];
591   	Element[5 ] = temp.Element[5 ];
592   	Element[6 ] = temp.Element[9 ];
593 
594   	Element[8 ] = temp.Element[2 ];
595   	Element[9 ] = temp.Element[6 ];
596   	Element[10] = temp.Element[10];
597 
598   	Element[12] *= -1.0f;
599   	Element[13] *= -1.0f;
600   	Element[14] *= -1.0f;
601 }
602 
MatrixFromAxisAngle(VECTOR axis,double theta)603 GLvoid MATRIX::MatrixFromAxisAngle(VECTOR axis, double theta)
604 {
605         QUAT q;
606         double halfTheta =  theta * 0.5;
607         double cosHalfTheta =  (cos(halfTheta));
608         double sinHalfTheta =  sin(halfTheta);
609         double xs, ys, zs, wx, wy, wz, xx, xy, xz, yy, yz, zz;
610         q.x = axis.x * sinHalfTheta;
611         q.y = axis.y * sinHalfTheta;
612         q.z = axis.z * sinHalfTheta;
613         q.w = cosHalfTheta;
614         xs = q.x * 2;  ys = q.y * 2;  zs = q.z * 2;
615         wx = q.w * xs; wy = q.w * ys; wz = q.w * zs;
616         xx = q.x * xs; xy = q.x * ys; xz = q.x * zs;
617         yy = q.y * ys; yz = q.y * zs; zz = q.z * zs;
618         Element[0] = 1 - (yy + zz);
619         Element[1] = xy - wz;
620         Element[2] = xz + wy;
621         Element[4] = xy + wz;
622         Element[5] = 1 - (xx + zz);
623         Element[6] = yz - wx;
624         Element[8] = xz - wy;
625         Element[9] = yz + wx;
626         Element[10] = 1 - (xx + yy);
627         Element[12] = Element[13] = Element[14] = Element[3] = Element[7] = Element[11] = 0;
628         Element[15] = 1;
629 }
630 
QuatToMatrix(QUAT quat)631 GLvoid MATRIX::QuatToMatrix(QUAT quat)
632 {
633   	double wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
634   	// calculate coefficients
635   	x2 = quat.x + quat.x;
636   	y2 = quat.y + quat.y;
637   	z2 = quat.z + quat.z;
638   	xx = quat.x * x2;
639   	xy = quat.x * y2;
640   	xz = quat.x * z2;
641   	yy = quat.y * y2;
642   	yz = quat.y * z2;
643   	zz = quat.z * z2;
644   	wx = quat.w * x2;
645   	wy = quat.w * y2;
646   	wz = quat.w * z2;
647   	Element[0] = (double) (1.0 - (yy + zz));
648   	Element[1] = (double) (xy - wz);
649   	Element[2] = (double) (xz + wy);
650   	Element[3] = 0.0;
651   	Element[4] = xy + wz;
652   	Element[5] = (double) (1.0 - (xx + zz));
653   	Element[6] = yz - wx;
654   	Element[7] = 0.0;
655   	Element[8] = xz - wy;
656   	Element[9] = yz + wx;
657   	Element[10] = (double) (1.0 - (xx + yy));
658   	Element[11] = 0.0;
659   	Element[12] = 0;
660   	Element[13] = 0;
661   	Element[14] = 0;
662   	Element[15] = 1;
663 }
664 
665 
666 
667 
668 
669 
QUAT(double sx,double sy,double sz,double sw)670 QUAT::QUAT(double sx, double sy, double sz, double sw)
671 :
672 	x(sx),
673 	y(sy),
674 	z(sz),
675     w(sw)
676 {
677 }
678 
~QUAT()679 QUAT::~QUAT()
680 {
681 }
682 
Reset()683 GLvoid QUAT::Reset()
684 {
685 	x = 0.0;
686 	y = 0.0;
687 	z = 0.0;
688 	w = 1.0;
689 }
690 
CopyQuat(QUAT q)691 GLvoid QUAT::CopyQuat(QUAT q)
692 {
693 	x = q.x;
694 	y = q.y;
695 	z = q.z;
696 	w = q.w;
697 }
698 
AxisAngleToQuat(VECTOR axis,double theta)699 GLvoid QUAT::AxisAngleToQuat(VECTOR axis, double theta)
700 {
701         double halfTheta = (double) (theta * 0.5);
702         double cosHalfTheta = (double) (cos(halfTheta));
703         double sinHalfTheta = (double) (sin(halfTheta));
704         x = axis.x * sinHalfTheta;
705         y = axis.y * sinHalfTheta;
706         z = axis.z * sinHalfTheta;
707         w = cosHalfTheta;
708 }
709 
EulerToQuat(double roll,double pitch,double yaw)710 GLvoid QUAT::EulerToQuat(double roll, double pitch, double yaw)
711 {
712     	double cr, cp, cy, sr, sp, sy, cpcy, spsy;  // calculate trig identities
713     	cr = (double) (cos(roll/2));
714     	cp = (double) (cos(pitch/2));
715     	cy = (double) (cos(yaw/2));
716     	sr = (double) (sin(roll/2));
717     	sp = (double) (sin(pitch/2));
718     	sy = (double) (sin(yaw/2));
719     	cpcy = cp * cy;
720     	spsy = sp * sy;
721     	w = cr * cpcy + sr * spsy;
722     	x = sr * cpcy - cr * spsy;
723     	y = cr * sp * cy + sr * cp * sy;
724     	z = cr * cp * sy - sr * sp * cy;
725 }
726 
MagnitudeQuat()727 double QUAT::MagnitudeQuat()
728 {
729   	return( (double) (sqrt(w*w+x*x+y*y+z*z)));
730 }
731 
NormaliseQuat()732 GLvoid QUAT::NormaliseQuat()
733 {
734   	double Mag;
735   	Mag = MagnitudeQuat();
736   	w = w/Mag;
737   	x = x/Mag;
738   	y = y/Mag;
739   	z = z/Mag;
740 }
741 
MultQuat(QUAT q)742 GLvoid QUAT::MultQuat(QUAT q)
743 {
744   	QUAT q3;
745   	VECTOR vectorq1;
746   	VECTOR vectorq2;
747   	vectorq1.x = x;
748   	vectorq1.y = y;
749   	vectorq1.z = z;
750   	vectorq2.x = q.x;
751   	vectorq2.y = q.y;
752   	vectorq2.z = q.z;
753 
754   	VECTOR tempvec1;
755   	VECTOR tempvec2;
756   	VECTOR tempvec3;
757   	tempvec1 = vectorq1;
758   	q3.w = (w*q.w) - tempvec1.DotProduct(vectorq2);
759   	tempvec1.CrossVector(vectorq2);
760   	tempvec2.x = w * q.x;
761   	tempvec2.y = w * q.y;
762   	tempvec2.z = w * q.z;
763   	tempvec3.x = q.w * x;
764   	tempvec3.y = q.w * y;
765   	tempvec3.z = q.w * z;
766   	q3.x = tempvec1.x + tempvec2.x + tempvec3.x;
767   	q3.y = tempvec1.y + tempvec2.y + tempvec3.y;
768   	q3.z = tempvec1.z + tempvec2.z + tempvec3.z;
769   	CopyQuat(q3);
770 }
771 
772 
773 
OBJECT()774 OBJECT::OBJECT()
775 {
776 }
777 
~OBJECT()778 OBJECT::~OBJECT()
779 {
780 }
781 
Reset()782 GLvoid OBJECT::Reset()
783 {
784 	Orientation.Reset();
785 	Position.Reset();
786         Delta_x = 0.0;
787         Delta_y = 0.0;
788         Delta_z = 0.0;
789         Movement_x = 0.0;
790         Movement_y = 0.0;
791         Movement_z = 0.0;
792 }
793 
Rotate()794 GLvoid OBJECT::Rotate()
795 {
796 	QUAT temp_quat;
797 	temp_quat.EulerToQuat(Delta_x * radian * Multiplier, Delta_y * radian * Multiplier, Delta_z * radian * Multiplier);
798 	Orientation.MultQuat(temp_quat);
799 }
800 
Draw()801 GLvoid OBJECT::Draw()
802 {
803   // Should probably be a pure virtual
804 }
805 
UpdatePosition()806 GLvoid OBJECT::UpdatePosition()
807 {
808     if (Movement_x != 0)
809         MoveX();
810     if (Movement_y != 0)
811         MoveY();
812     if (Movement_z != 0)
813         MoveZ();
814 
815     Movement_x = 0.0;
816     Movement_y = 0.0;
817     Movement_z = 0.0;
818 }
819 
UpdatePosition(double x,double y,double z)820 GLvoid OBJECT::UpdatePosition(double x, double y, double z)
821 {
822     if (x != 0)
823         MoveX(x);
824     if (y != 0)
825         MoveY(y);
826     if (z != 0)
827         MoveZ(z);
828 }
829 
MoveX()830 GLvoid OBJECT::MoveX()
831 {
832     double DirX;
833     double DirY;
834     double DirZ;
835     double W;
836     double X;
837     double Y;
838     double Z;
839     QUAT TempQuat;
840     QUAT TempQuat2;
841     TempQuat = Orientation;
842     TempQuat2.EulerToQuat(0.0,(double) (-90.0*(pi/180)), 0.0);
843     TempQuat.MultQuat(TempQuat2);
844     W = TempQuat.w;
845     X = TempQuat.x;
846     Y = TempQuat.y;
847     Z = TempQuat.z;
848     DirX = (double)(2.0 * ( X * Z - W * Y ));
849     DirY = (double)(2.0 * ( Y * Z + W * X ));
850     DirZ = (double)(1.0 - 2.0 * ( X * X + Y * Y ));
851     Position.x += DirX * Movement_x * Multiplier;
852     Position.y += DirY * Movement_x * Multiplier;
853     Position.z += DirZ * Movement_x * Multiplier;
854 }
855 
MoveY()856 GLvoid OBJECT::MoveY()
857 {
858     double DirX;
859     double DirY;
860     double DirZ;
861     double W;
862     double X;
863     double Y;
864     double Z;
865     QUAT TempQuat;
866     QUAT TempQuat2;
867     TempQuat = Orientation;
868     TempQuat2.EulerToQuat((double)(90.0*(pi/180)), 0.0, 0.0);
869     TempQuat.MultQuat(TempQuat2);
870     W = TempQuat.w;
871     X = TempQuat.x;
872     Y = TempQuat.y;
873     Z = TempQuat.z;
874     DirX = (double)(2.0 * ( X * Z - W * Y ));
875     DirY = (double)(2.0 * ( Y * Z + W * X ));
876     DirZ = (double)(1.0 - 2.0 * ( X * X + Y * Y ));
877     Position.x += DirX * Movement_y * Multiplier;
878     Position.y += DirY * Movement_y * Multiplier;
879     Position.z += DirZ * Movement_y * Multiplier;
880 }
881 
MoveZ()882 GLvoid OBJECT::MoveZ()
883 {
884     double DirX;
885     double DirY;
886     double DirZ;
887     double W = Orientation.w;
888     double X = Orientation.x;
889     double Y = Orientation.y;
890     double Z = Orientation.z;
891     DirX = 2.0 * ( X * Z - W * Y );
892     DirY = 2.0 * ( Y * Z + W * X );
893     DirZ = 1.0 - 2.0 * ( X * X + Y * Y );
894     Position.x += DirX * Movement_z * Multiplier;
895     Position.y += DirY * Movement_z * Multiplier;
896     Position.z += DirZ * Movement_z * Multiplier;
897 }
898 
MoveX(double x)899 GLvoid OBJECT::MoveX(double x)
900 {
901     double DirX;
902     double DirY;
903     double DirZ;
904     double W;
905     double X;
906     double Y;
907     double Z;
908     QUAT TempQuat;
909     QUAT TempQuat2;
910     TempQuat = Orientation;
911     TempQuat2.EulerToQuat(0.0, -90.0*(pi/180), 0.0);
912     TempQuat.MultQuat(TempQuat2);
913     W = TempQuat.w;
914     X = TempQuat.x;
915     Y = TempQuat.y;
916     Z = TempQuat.z;
917     DirX = 2.0 * ( X * Z - W * Y );
918     DirY = 2.0 * ( Y * Z + W * X );
919     DirZ = 1.0 - 2.0 * ( X * X + Y * Y );
920     Position.x += DirX * x;
921     Position.y += DirY * x;
922     Position.z += DirZ * x;
923 }
924 
MoveY(double y)925 GLvoid OBJECT::MoveY(double y)
926 {
927     double DirX;
928     double DirY;
929     double DirZ;
930     double W;
931     double X;
932     double Y;
933     double Z;
934     QUAT TempQuat;
935     QUAT TempQuat2;
936     TempQuat = Orientation;
937     TempQuat2.EulerToQuat(90.0*(pi/180), 0.0, 0.0);
938     TempQuat.MultQuat(TempQuat2);
939     W = TempQuat.w;
940     X = TempQuat.x;
941     Y = TempQuat.y;
942     Z = TempQuat.z;
943     DirX = 2.0 * ( X * Z - W * Y );
944     DirY = 2.0 * ( Y * Z + W * X );
945     DirZ = 1.0 - 2.0 * ( X * X + Y * Y );
946     Position.x += DirX * y;
947     Position.y += DirY * y;
948     Position.z += DirZ * y;
949 }
950 
MoveZ(double z)951 GLvoid OBJECT::MoveZ(double z)
952 {
953     double DirX;
954     double DirY;
955     double DirZ;
956     double W = Orientation.w;
957     double X = Orientation.x;
958     double Y = Orientation.y;
959     double Z = Orientation.z;
960     DirX = 2.0 * ( X * Z - W * Y );
961     DirY = 2.0 * ( Y * Z + W * X );
962     DirZ = 1.0 - 2.0 * ( X * X + Y * Y );
963     Position.x += DirX * z;
964     Position.y += DirY * z;
965     Position.z += DirZ * z;
966 }
967 
GetXUnit()968 VECTOR OBJECT::GetXUnit()
969 {
970     double DirX;
971     double DirY;
972     double DirZ;
973     double W;
974     double X;
975     double Y;
976     double Z;
977     QUAT TempQuat;
978     QUAT TempQuat2;
979     TempQuat = Orientation;
980     TempQuat2.EulerToQuat(0.0, -90.0*(pi/180), 0.0);
981     TempQuat.MultQuat(TempQuat2);
982     W = TempQuat.w;
983     X = TempQuat.x;
984     Y = TempQuat.y;
985     Z = TempQuat.z;
986     DirX = 2.0 * ( X * Z - W * Y );
987     DirY = 2.0 * ( Y * Z + W * X );
988     DirZ = 1.0 - 2.0 * ( X * X + Y * Y );
989     VECTOR Unit;
990     Unit.x += DirX * 1;
991     Unit.y += DirY * 1;
992     Unit.z += DirZ * 1;
993     return Unit;
994 }
995 
GetYUnit()996 VECTOR OBJECT::GetYUnit()
997 {
998     double DirX;
999     double DirY;
1000     double DirZ;
1001     double W;
1002     double X;
1003     double Y;
1004     double Z;
1005     QUAT TempQuat;
1006     QUAT TempQuat2;
1007     TempQuat = Orientation;
1008     TempQuat2.EulerToQuat(90.0*(pi/180), 0.0, 0.0);
1009     TempQuat.MultQuat(TempQuat2);
1010     W = TempQuat.w;
1011     X = TempQuat.x;
1012     Y = TempQuat.y;
1013     Z = TempQuat.z;
1014     DirX = 2.0 * ( X * Z - W * Y );
1015     DirY = 2.0 * ( Y * Z + W * X );
1016     DirZ = 1.0 - 2.0 * ( X * X + Y * Y );
1017     VECTOR Unit;
1018     Unit.x += DirX * 1;
1019     Unit.y += DirY * 1;
1020     Unit.z += DirZ * 1;
1021     return Unit;
1022 }
1023 
GetZUnit()1024 VECTOR OBJECT::GetZUnit()
1025 {
1026     double DirX;
1027     double DirY;
1028     double DirZ;
1029     double W = Orientation.w;
1030     double X = Orientation.x;
1031     double Y = Orientation.y;
1032     double Z = Orientation.z;
1033     DirX = 2.0 * ( X * Z - W * Y );
1034     DirY = 2.0 * ( Y * Z + W * X );
1035     DirZ = 1.0 - 2.0 * ( X * X + Y * Y );
1036     VECTOR Unit;
1037     Unit.x += DirX * 1;
1038     Unit.y += DirY * 1;
1039     Unit.z += DirZ * 1;
1040     return Unit;
1041 }
1042 
1043 
1044 
CAMERA()1045 CAMERA::CAMERA()
1046 {
1047 }
1048 
~CAMERA()1049 CAMERA::~CAMERA()
1050 {
1051 }
1052 
Reset()1053 GLvoid CAMERA::Reset()
1054 {
1055 	Orientation.Reset();
1056 	Position.Reset();
1057         Delta_x = 0.0;
1058         Delta_y = 0.0;
1059         Delta_z = 0.0;
1060         Matrix.LoadIdentity();
1061 }
1062 
Update()1063 GLvoid CAMERA::Update()
1064 {
1065     Rotate();
1066 
1067     UpdatePosition();
1068 }
1069 
Apply()1070 GLvoid CAMERA::Apply()
1071 {
1072     Matrix.QuatToMatrix(Orientation);
1073     Matrix.MatrixInverse();
1074 
1075     glLoadMatrixd(Matrix.Element);
1076     glTranslated(-Position.x,-Position.y,-Position.z);
1077 }
1078 
1079 
1080 
LoadIdentity(double m[])1081 GLvoid LoadIdentity(double m[])
1082 {
1083     	m[0]=1.0f;
1084     	m[1]=0.0f;
1085     	m[2]=0.0f;
1086     	m[3]=0.0f;
1087 
1088     	m[4]=0.0f;
1089     	m[5]=1.0f;
1090     	m[6]=0.0f;
1091     	m[7]=0.0f;
1092 
1093     	m[8]=0.0f;
1094     	m[9]=0.0f;
1095     	m[10]=1.0f;
1096     	m[11]=0.0f;
1097 
1098     	m[12]=0.0f;
1099     	m[13]=0.0f;
1100     	m[14]=0.0f;
1101     	m[15]=1.0f;
1102 }
1103 
CopyMatrix(double m[],double n[])1104 GLvoid CopyMatrix(double m[], double n[])
1105 {
1106     	m[0 ] = n[0 ];
1107     	m[1 ] = n[1 ];
1108     	m[2 ] = n[2 ];
1109     	m[3 ] = n[3 ];
1110     	m[4 ] = n[4 ];
1111     	m[5 ] = n[5 ];
1112     	m[6 ] = n[6 ];
1113     	m[7 ] = n[7 ];
1114     	m[8 ] = n[8 ];
1115     	m[9 ] = n[9 ];
1116     	m[10] = n[10];
1117     	m[11] = n[11];
1118     	m[12] = n[12];
1119     	m[13] = n[13];
1120     	m[14] = n[14];
1121     	m[15] = n[15];
1122 }
1123 
MultMatrix(double m[],double n[])1124 GLvoid MultMatrix(double m[], double n[])
1125 {
1126     	double temp[16];
1127 
1128     	CopyMatrix(temp, m);
1129   	m[0] = temp[0 ]*n[0 ]
1130        		+ temp[4 ]*n[1 ]
1131        		+ temp[8 ]*n[2 ]
1132        		+ temp[12]*n[3 ];
1133 
1134   	m[1] = temp[1 ]*n[0 ]
1135        		+ temp[5 ]*n[1 ]
1136        		+ temp[9 ]*n[2 ]
1137        		+ temp[13]*n[3 ];
1138 
1139   	m[2] = temp[2 ]*n[0 ]
1140        		+ temp[6 ]*n[1 ]
1141        		+ temp[10]*n[2 ]
1142        		+ temp[14]*n[3 ];
1143 
1144   	m[3] = temp[3 ]*n[0 ]
1145        		+ temp[7 ]*n[1 ]
1146        		+ temp[11]*n[2 ]
1147        		+ temp[15]*n[3 ];
1148 
1149   	m[4] = temp[0 ]*n[4 ]
1150        		+ temp[4 ]*n[5 ]
1151        		+ temp[8 ]*n[6 ]
1152        		+ temp[12]*n[7 ];
1153 
1154   	m[5] = temp[1 ]*n[4 ]
1155        		+ temp[5 ]*n[5 ]
1156        		+ temp[9 ]*n[6 ]
1157        		+ temp[13]*n[7 ];
1158 
1159   	m[6] = temp[2 ]*n[4 ]
1160        		+ temp[6 ]*n[5 ]
1161        		+ temp[10]*n[6 ]
1162        		+ temp[14]*n[7 ];
1163 
1164   	m[7] = temp[3 ]*n[4 ]
1165        		+ temp[7 ]*n[5 ]
1166        		+ temp[11]*n[6 ]
1167        		+ temp[15]*n[7 ];
1168 
1169   	m[8] = temp[0 ]*n[8 ]
1170        		+ temp[4 ]*n[9 ]
1171        		+ temp[8 ]*n[10]
1172        		+ temp[12]*n[11];
1173 
1174   	m[9] = temp[1 ]*n[8 ]
1175        		+ temp[5 ]*n[9 ]
1176        		+ temp[9 ]*n[10]
1177        		+ temp[13]*n[11];
1178 
1179   	m[10]= temp[2 ]*n[8 ]
1180        		+ temp[6 ]*n[9 ]
1181        		+ temp[10]*n[10]
1182        		+ temp[14]*n[11];
1183 
1184   	m[11]= temp[3 ]*n[8 ]
1185        		+ temp[7 ]*n[9 ]
1186        		+ temp[11]*n[10]
1187        		+ temp[15]*n[11];
1188 
1189   	m[12]= temp[0 ]*n[12]
1190        		+ temp[4 ]*n[13]
1191        		+ temp[8 ]*n[14]
1192        		+ temp[12]*n[15];
1193 
1194   	m[13]= temp[1 ]*n[12]
1195        		+ temp[5 ]*n[13]
1196        		+ temp[9 ]*n[14]
1197        		+ temp[13]*n[15];
1198 
1199   	m[14]= temp[2 ]*n[12]
1200        		+ temp[6 ]*n[13]
1201        		+ temp[10]*n[14]
1202        		+ temp[14]*n[15];
1203 
1204   	m[15]= temp[3 ]*n[12]
1205        		+ temp[7 ]*n[13]
1206        		+ temp[11]*n[14]
1207        		+ temp[15]*n[15];
1208 }
1209 
MatrixInverse(double m[])1210 GLvoid MatrixInverse(double m[])
1211 {
1212 	double n[16];
1213 
1214   	CopyMatrix(n, m);
1215   	m[0 ] = n[0 ];
1216   	m[1 ] = n[4 ];
1217   	m[2 ] = n[8 ];
1218 
1219   	m[4 ] = n[1 ];
1220   	m[5 ] = n[5 ];
1221   	m[6 ] = n[9 ];
1222 
1223   	m[8 ] = n[2 ];
1224   	m[9 ] = n[6 ];
1225   	m[10] = n[10];
1226 
1227   	m[12] *= -1.0f;
1228   	m[13] *= -1.0f;
1229   	m[14] *= -1.0f;
1230 }
1231 
1232     /* The following routine converts an angle and a unit axis vector
1233         * to a matrix, returning the corresponding unit quaternion at no
1234         * extra cost. It is written in such a way as to allow both fixed
1235         * point and doubleing point versions to be created by appropriate
1236         * definitions of FPOINT, ANGLE, VECTOR, QUAT, MATRIX, MUL, HALF,
1237         * TWICE, COS, SIN, ONE, and ZERO.
1238         * The following is an example of doubleing point definitions.*/
AxisAngleToMatrix(VECTOR axis,double theta,double m[16])1239 QUAT AxisAngleToMatrix(VECTOR axis, double theta, double m[16])
1240 {
1241         QUAT q;
1242         double halfTheta = theta * 0.5;
1243         double cosHalfTheta = cos(halfTheta);
1244         double sinHalfTheta = sin(halfTheta);
1245         double xs, ys, zs, wx, wy, wz, xx, xy, xz, yy, yz, zz;
1246         q.x = axis.x * sinHalfTheta;
1247         q.y = axis.y * sinHalfTheta;
1248         q.z = axis.z * sinHalfTheta;
1249         q.w = cosHalfTheta;
1250         xs = q.x * 2;  ys = q.y * 2;  zs = q.z * 2;
1251         wx = q.w * xs; wy = q.w * ys; wz = q.w * zs;
1252         xx = q.x * xs; xy = q.x * ys; xz = q.x * zs;
1253         yy = q.y * ys; yz = q.y * zs; zz = q.z * zs;
1254         m[0] = 1 - (yy + zz);
1255         m[1] = xy - wz;
1256         m[2] = xz + wy;
1257         m[4] = xy + wz;
1258         m[5] = 1 - (xx + zz);
1259         m[6] = yz - wx;
1260         m[8] = xz - wy;
1261         m[9] = yz + wx;
1262         m[10] = 1 - (xx + yy);
1263         /* Fill in remainder of 4x4 homogeneous transform matrix. */
1264         m[12] = m[13] = m[14] = m[3] = m[7] = m[11] = 0;
1265         m[15] = 1;
1266         return (q);
1267 }
1268 
DotProduct(VECTOR vec1,VECTOR vec2)1269 double DotProduct(VECTOR vec1, VECTOR vec2)
1270 {
1271     /*
1272     Dot Product of two Vectors.
1273 
1274     U = (Ux, Uy, Uz)
1275     V = (Vx, Vy, Vz)
1276     U*V = UxVx + UyVy + UzVz
1277     U*V = |U||V|cos(t) (where t is the angle theta between the two vectors)
1278     */
1279   	double dot;
1280   	dot = vec1.x * vec2.x + vec1.y * vec2.y + vec1.z * vec2.z;
1281   	return dot;
1282 }
1283 
CrossVector(VECTOR vec1,VECTOR vec2)1284 VECTOR CrossVector(VECTOR vec1, VECTOR vec2)
1285 {
1286     /*
1287     Cross Product of Two Vectors.
1288 
1289     a �b = ( a.y * b.z - a.z * b.y,
1290 
1291     a.z * b.x - a.x * b.z,
1292 
1293     a.x * b.y - a.y * b.x )
1294 
1295     | a �b | = |a| * |b| * sin()
1296     */
1297   	VECTOR vec3;
1298   	vec3.x = vec1.y * vec2.z - vec1.z * vec2.y;
1299   	vec3.y = vec1.z * vec2.x - vec1.x * vec2.z;
1300   	vec3.z = vec1.x * vec2.y - vec1.y * vec2.x;
1301   	return vec3;
1302 }
1303 
EulerToQuat(double roll,double pitch,double yaw,QUAT * quat)1304 void EulerToQuat(double roll, double pitch, double yaw, QUAT * quat)
1305 {
1306     /*
1307     Euler Angle to Quarternion.
1308 
1309     q = qyaw qpitch qroll where:
1310 
1311     qyaw = [cos(f /2), (0, 0, sin(f /2)]
1312     qpitch = [cos (q/2), (0, sin(q/2), 0)]
1313     qroll = [cos (y/2), (sin(y/2), 0, 0)]
1314     */
1315 	double cr, cp, cy, sr, sp, sy, cpcy, spsy;  // calculate trig identities
1316     	cr = cos(roll/2);
1317 	cp = cos(pitch/2);
1318     	cy = cos(yaw/2);
1319     	sr = sin(roll/2);
1320 	sp = sin(pitch/2);
1321     	sy = sin(yaw/2);
1322     	cpcy = cp * cy;
1323     	spsy = sp * sy;
1324 	quat->w = cr * cpcy + sr * spsy;
1325     	quat->x = sr * cpcy - cr * spsy;
1326 	quat->y = cr * sp * cy + sr * cp * sy;
1327     	quat->z = cr * cp * sy - sr * sp * cy;
1328 }
1329 
MagnitudeQuat(QUAT q1)1330 double MagnitudeQuat(QUAT q1)
1331 {
1332   	return( sqrt(q1.w*q1.w+q1.x*q1.x+q1.y*q1.y+q1.z*q1.z));
1333 }
1334 
NormaliseQuat(QUAT q1)1335 QUAT NormaliseQuat(QUAT q1)
1336 {
1337   	QUAT q2;
1338   	double Mag;
1339   	Mag = MagnitudeQuat(q1);
1340   	q2.w = q1.w/Mag;
1341   	q2.x = q1.x/Mag;
1342   	q2.y = q1.y/Mag;
1343   	q2.z = q1.z/Mag;
1344   	return q2;
1345 }
1346 
QuatToMatrix(QUAT quat,double m[16])1347 GLvoid QuatToMatrix(QUAT quat, double m[16])
1348 {
1349   	double wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
1350   	// calculate coefficients
1351   	x2 = quat.x + quat.x;
1352   	y2 = quat.y + quat.y;
1353   	z2 = quat.z + quat.z;
1354   	xx = quat.x * x2;
1355   	xy = quat.x * y2;
1356   	xz = quat.x * z2;
1357   	yy = quat.y * y2;
1358   	yz = quat.y * z2;
1359   	zz = quat.z * z2;
1360   	wx = quat.w * x2;
1361   	wy = quat.w * y2;
1362   	wz = quat.w * z2;
1363   	m[0] = 1.0 - (yy + zz);
1364   	m[1] = xy - wz;
1365   	m[2] = xz + wy;
1366   	m[3] = 0.0;
1367   	m[4] = xy + wz;
1368   	m[5] = 1.0 - (xx + zz);
1369   	m[6] = yz - wx;
1370   	m[7] = 0.0;
1371   	m[8] = xz - wy;
1372   	m[9] = yz + wx;
1373   	m[10] = 1.0 - (xx + yy);
1374   	m[11] = 0.0;
1375   	m[12] = 0;
1376   	m[13] = 0;
1377   	m[14] = 0;
1378   	m[15] = 1;
1379 }
1380 
MultQuat(QUAT q1,QUAT q2)1381 QUAT MultQuat(QUAT q1, QUAT q2)
1382 {
1383     /*
1384     Multiplication of two Quarternions.
1385 
1386     qq = [ww - v  v, v x v + wv +wv]
1387     (  is vector dot product and x is vector cross product )
1388     */
1389   	QUAT q3;
1390   	VECTOR vectorq1;
1391   	VECTOR vectorq2;
1392   	vectorq1.x = q1.x;
1393   	vectorq1.y = q1.y;
1394   	vectorq1.z = q1.z;
1395   	vectorq2.x = q2.x;
1396   	vectorq2.y = q2.y;
1397   	vectorq2.z = q2.z;
1398 
1399   	VECTOR tempvec1;
1400   	VECTOR tempvec2;
1401   	VECTOR tempvec3;
1402   	q3.w = (q1.w*q2.w) - DotProduct(vectorq1, vectorq2);
1403   	tempvec1 = CrossVector(vectorq1, vectorq2);
1404   	tempvec2.x = q1.w * q2.x;
1405   	tempvec2.y = q1.w * q2.y;
1406   	tempvec2.z = q1.w * q2.z;
1407   	tempvec3.x = q2.w * q1.x;
1408   	tempvec3.y = q2.w * q1.y;
1409   	tempvec3.z = q2.w * q1.z;
1410   	q3.x = tempvec1.x + tempvec2.x + tempvec3.x;
1411   	q3.y = tempvec1.y + tempvec2.y + tempvec3.y;
1412   	q3.z = tempvec1.z + tempvec2.z + tempvec3.z;
1413   	return NormaliseQuat(q3);
1414 }
1415 
1416 // returns the normal vector to a plane defined by three vertices
GetNormal(VECTOR vertex1,VECTOR vertex2,VECTOR vertex3)1417 VECTOR GetNormal(VECTOR vertex1, VECTOR vertex2, VECTOR vertex3)
1418 {
1419     	double ux, uy, uz, vx, vy, vz;
1420       	VECTOR temp_vertex;
1421 
1422       	ux = vertex1.x - vertex2.x;
1423       	uy = vertex1.y - vertex2.y;
1424       	uz = vertex1.z - vertex2.z;
1425       	vx = vertex3.x - vertex2.x;
1426       	vy = vertex3.y - vertex2.y;
1427       	vz = vertex3.z - vertex2.z;
1428       	temp_vertex.x = (uy*vz)-(vy*uz);
1429       	temp_vertex.y = (uz*vx)-(vz*ux);
1430       	temp_vertex.z = (ux*vy)-(vx*uy);
1431       	return temp_vertex;
1432 }
1433 
GetNorm(double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3)1434 VERTEX GetNorm(double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3)
1435 {
1436     	double ux;
1437     	double uy;
1438     	double uz;
1439     	double vx;
1440     	double vy;
1441     	double vz;
1442       	VERTEX temp_vertex;
1443       	ux = x1 - x2;
1444       	uy = y1 - y2;
1445       	uz = z1 - z2;
1446       	vx = x3 - x2;
1447       	vy = y3 - y2;
1448       	vz = z3 - z2;
1449       	temp_vertex.normal.x = (uy*vz)-(vy*uz);
1450       	temp_vertex.normal.y = (uz*vx)-(vz*ux);
1451       	temp_vertex.normal.z = (ux*vy)-(vx*uy);
1452       	return temp_vertex;
1453 }
1454 
MagnitudeVector(VECTOR vec1)1455 double MagnitudeVector(VECTOR vec1)
1456 {
1457   return(sqrt(vec1.x*vec1.x+vec1.y*vec1.y+vec1.z*vec1.z));
1458 }
1459 
GetUnitVector(VECTOR vector)1460 VECTOR GetUnitVector(VECTOR vector)
1461 {
1462 	// Reduces a normal vector specified as a set of three coordinates,
1463 	// to a unit normal vector of length one.
1464 
1465 	// Calculate the length of the vector
1466 	double length = (double) sqrt(( vector.x * vector.x) +
1467 						        ( vector.y * vector.y) +
1468 						        ( vector.z * vector.z) );
1469 
1470 	// Keep the program from blowing up by providing an exceptable
1471 	// value for vectors that may calculated too close to zero.
1472 	if(length == 0.0f)
1473 		length = 1.0f;
1474 
1475 	// Dividing each element by the length will result in a
1476 	// unit normal vector.
1477 	vector.x /= length;
1478 	vector.y /= length;
1479 	vector.z /= length;
1480         return vector;
1481 }
1482 
GetEdgeVector(VECTOR point1,VECTOR point2)1483 VECTOR GetEdgeVector(VECTOR point1, VECTOR point2)
1484 {
1485   VECTOR temp_vector;
1486   temp_vector.x = point1.x - point2.x;
1487   temp_vector.y = point1.y - point2.y;
1488   temp_vector.z = point1.z - point2.z;
1489   return temp_vector;
1490 }
1491 
1492 
1493 
CheckPointInTriangle(VECTOR point,VECTOR a,VECTOR b,VECTOR c)1494 bool CheckPointInTriangle(VECTOR point, VECTOR a, VECTOR b, VECTOR c)
1495 {
1496 
1497   double total_angles = 0.0f;
1498 
1499   // make the 3 vectors
1500   VECTOR TempVect;
1501   TempVect.x = point.x - a.x;
1502   TempVect.y = point.y - a.y;
1503   TempVect.z = point.z - a.z;
1504   VECTOR v1 = TempVect;
1505   TempVect.x = point.x - b.x;
1506   TempVect.y = point.y - b.y;
1507   TempVect.z = point.z - b.z;
1508   VECTOR v2 = TempVect;
1509   TempVect.x = point.x - c.x;
1510   TempVect.y = point.y - c.y;
1511   TempVect.z = point.z - c.z;
1512   VECTOR v3 = TempVect;
1513 
1514   v1 = GetUnitVector(v1);
1515   v2 = GetUnitVector(v2);
1516   v3 = GetUnitVector(v3);
1517   double Dot1 = DotProduct(v1,v2);
1518   if (Dot1 < -1)
1519     Dot1 = -1;
1520   if (Dot1 > 1)
1521     Dot1 = 1;
1522   total_angles += acos(Dot1);
1523   double Dot2 = DotProduct(v2,v3);
1524   if (Dot2 < -1)
1525     Dot2 = -1;
1526   if (Dot2 > 1)
1527     Dot2 = 1;
1528   total_angles += acos(Dot2);
1529   double Dot3 = DotProduct(v3,v1);
1530   if (Dot3 < -1)
1531     Dot3 = -1;
1532   if (Dot3 > 1)
1533     Dot3 = 1;
1534   total_angles += acos(Dot3);
1535 
1536   if (fabs(total_angles-2*pi) <= 0.005)
1537    return (true);
1538 
1539   return(false);
1540 }
1541 
1542 // Check if the point is in the quad
1543 // All the vertices of the quad MUST be in the same plane
CheckPointInQuad(VECTOR * point,QUAD * quad)1544 bool CheckPointInQuad(VECTOR* point, QUAD* quad)
1545 {
1546 	if (CheckPointInTriangle(*point,quad->Vertex[0].coords,quad->Vertex[1].coords,quad->Vertex[2].coords))
1547 		return true;
1548 
1549 	if (CheckPointInTriangle(*point,quad->Vertex[0].coords,quad->Vertex[3].coords,quad->Vertex[2].coords))
1550 		return true;
1551 
1552 	return false;
1553 }
1554 
1555 
1556 
line_plane_collision(VECTOR * a,VECTOR * b,QUAD * plane)1557 VECTOR line_plane_collision(VECTOR *a, VECTOR *b, QUAD *plane)
1558 {
1559     double Distance = -(plane->Vertex[0].coords.x * plane->Vertex[0].normal.x + plane->Vertex[0].coords.y * plane->Vertex[0].normal.y + plane->Vertex[0].coords.z * plane->Vertex[0].normal.z);
1560     double final_x,final_y,final_z,final_t;
1561     double t,i;
1562     VECTOR temp;
1563 
1564     t=0; i=0;
1565     i+=(plane->Vertex[0].normal.x*b->x)+(plane->Vertex[0].normal.y*b->y)+(plane->Vertex[0].normal.z*b->z)+(Distance);
1566     t+=(plane->Vertex[0].normal.x*(b->x*-1))+(plane->Vertex[0].normal.y*(b->y*-1))+(plane->Vertex[0].normal.z*(b->z*-1));
1567     t+=(plane->Vertex[0].normal.x*a->x)+(plane->Vertex[0].normal.y*a->y)+(plane->Vertex[0].normal.z*a->z);
1568 
1569     // Be wary of possible divide-by-zeros here (i.e. if t==0)
1570     if (t)
1571 		final_t = (-i)/t;
1572     else
1573 		final_t = 0;
1574     // Vertical Line Segment
1575     if ((a->x == b->x)&&(a->z == b->z))
1576     { // vertical line segment
1577 
1578         temp.x = a->x;
1579         temp.y = (-((plane->Vertex[0].normal.x*a->x)+(plane->Vertex[0].normal.z*a->z)+(Distance)))/(plane->Vertex[0].normal.y);
1580         temp.z = a->z;
1581 
1582         return(temp);
1583     }
1584 
1585     final_x = (((a->x)*(final_t))+((b->x)*(1-final_t)));
1586     final_y = (((a->y)*(final_t))+((b->y)*(1-final_t)));
1587     final_z = (((a->z)*(final_t))+((b->z)*(1-final_t)));
1588 
1589     temp.x = final_x;
1590     temp.y = final_y;
1591     temp.z = final_z;
1592 
1593     return(temp);
1594 }
1595 
1596 
1597 
1598 
1599 
CheckClipPlanes(CAMERA Camera,VECTOR Vect)1600 bool CheckClipPlanes(CAMERA Camera, VECTOR Vect)
1601 {
1602   double ProjectionMatrix[16];
1603   double ModelViewMatrix[16];
1604   double A, B, C, Distance;
1605   int Counter = 0;
1606   glGetDoublev(GL_PROJECTION_MATRIX, ProjectionMatrix);
1607   glGetDoublev(GL_MODELVIEW_MATRIX, ModelViewMatrix);
1608 
1609   MultMatrix(ProjectionMatrix, ModelViewMatrix);
1610 
1611   //right clipping plane
1612   A = ProjectionMatrix[0] - ProjectionMatrix[3];
1613   B = ProjectionMatrix[4] - ProjectionMatrix[7];
1614   C = ProjectionMatrix[8] - ProjectionMatrix[11];
1615 
1616   Distance = -1 * (A * (-Camera.Position.x + Vect.x) + B * (-Camera.Position.y + Vect.y) + C * (-Camera.Position.z + Vect.z));
1617   if (Distance > 0)
1618     Counter++;
1619 
1620   //left clipping plane
1621   A = ProjectionMatrix[0] + ProjectionMatrix[3];
1622   B = ProjectionMatrix[4] + ProjectionMatrix[7];
1623   C = ProjectionMatrix[8] + ProjectionMatrix[11];
1624 
1625   Distance = A * (-Camera.Position.x + Vect.x) + B * (-Camera.Position.y + Vect.y) + C * (-Camera.Position.z + Vect.z);
1626   if (Distance > 0)
1627     Counter++;
1628 
1629   //top clipping plane
1630   A = ProjectionMatrix[1] - ProjectionMatrix[3];
1631   B = ProjectionMatrix[5] - ProjectionMatrix[7];
1632   C = ProjectionMatrix[9] - ProjectionMatrix[11];
1633 
1634   Distance = -1 * (A * (-Camera.Position.x + Vect.x) + B * (-Camera.Position.y + Vect.y) + C * (-Camera.Position.z + Vect.z));
1635   if (Distance > 0)
1636     Counter++;
1637 
1638   //bottom clipping plane
1639   A = ProjectionMatrix[1] + ProjectionMatrix[3];
1640   B = ProjectionMatrix[5] + ProjectionMatrix[7];
1641   C = ProjectionMatrix[9] + ProjectionMatrix[11];
1642 
1643   Distance = A * (-Camera.Position.x + Vect.x) + B * (-Camera.Position.y + Vect.y) + C * (-Camera.Position.z + Vect.z);
1644   if (Distance > 0)
1645     Counter++;
1646 
1647   // near clipping plane (might not be necessary when the near frustrum plane is close, but just in case)
1648   A = ProjectionMatrix[2] - ProjectionMatrix[3];
1649   B = ProjectionMatrix[6] - ProjectionMatrix[7];
1650   C = ProjectionMatrix[10] - ProjectionMatrix[11];
1651 
1652   Distance = A * (-Camera.Position.x + Vect.x) + B * (-Camera.Position.y + Vect.y) + C * (-Camera.Position.z + Vect.z);
1653   if (Distance > 0)
1654     Counter++;
1655 
1656   // far clipping plane (the equation didn't work for the far plane, so I'll just use a distance test)
1657   VECTOR Vect2;
1658   Vect2.x = Vect.x - Camera.Position.x;
1659   Vect2.y = Vect.y - Camera.Position.y;
1660   Vect2.z = Vect.z - Camera.Position.z;
1661   if (MagnitudeVector(Vect2) < 200)
1662     Counter++;
1663 
1664   if (Counter == 6)
1665     return 1;
1666   else
1667     return 0;
1668 }
1669 
1670 
1671 
1672 
1673 
1674 // SelectPolygon() creates a ray from the current camera coordinates to a point half a unit into the screen that is coincident
1675 // with the mouse coordinates passed in and checks for an intersection along the line of this ray with the quad passed in.
1676 // The 4 vertices of the quad MUST be in the same plane
1677 // If there is an intersection then the function returns true.
1678 // MouseX and MouseY are relative to the top-left corner of the window.
SelectPolygon(double cameraX,double cameraY,double cameraZ,double Vertices[12],double Intersec[3],int MouseX,int MouseY)1679 bool SelectPolygon(double cameraX, double cameraY, double cameraZ,
1680 				   double Vertices[12], double Intersec[3],
1681 				   int MouseX, int MouseY)
1682 {
1683 	VECTOR CollisionPoint;
1684 	VECTOR WorldPos;
1685 	VECTOR posCam=VECTOR(cameraX,cameraY,cameraZ);
1686 	VECTOR* posCamera=&posCam;
1687 	QUAD aquad;
1688 
1689 	aquad.Vertex[0].coords.x=Vertices[0];
1690 	aquad.Vertex[0].coords.y=Vertices[1];
1691 	aquad.Vertex[0].coords.z=Vertices[2];
1692 
1693 	aquad.Vertex[1].coords.x=Vertices[3];
1694 	aquad.Vertex[1].coords.y=Vertices[4];
1695 	aquad.Vertex[1].coords.z=Vertices[5];
1696 
1697 	aquad.Vertex[2].coords.x=Vertices[6];
1698 	aquad.Vertex[2].coords.y=Vertices[7];
1699 	aquad.Vertex[2].coords.z=Vertices[8];
1700 
1701 	aquad.Vertex[3].coords.x=Vertices[9];
1702 	aquad.Vertex[3].coords.y=Vertices[10];
1703 	aquad.Vertex[3].coords.z=Vertices[11];
1704 
1705 	aquad.numVertices=4;
1706 	aquad.SetNormal();
1707 	QUAD* quad=&aquad;
1708 
1709 	// Get the modelview and projection matrices
1710 	double WorldPosX, WorldPosY, WorldPosZ, MousePosX, MousePosY, MousePosZ;
1711 	double ModelMatrix[16];
1712 	glGetDoublev(GL_MODELVIEW_MATRIX, ModelMatrix);
1713 	double ProjMatrix[16];
1714 	glGetDoublev(GL_PROJECTION_MATRIX, ProjMatrix);
1715 
1716 	// Get the current viewport
1717 	int Viewport[4];
1718 	glGetIntegerv(GL_VIEWPORT, Viewport);
1719 
1720 	if (MouseX >= Viewport[0] && MouseX <= Viewport[2] && MouseY >= Viewport[1] && MouseY <= Viewport[3])
1721 	{
1722 		// Set the end point of ray in windows coordinates
1723 		MousePosX = MouseX;
1724 		MousePosY = Viewport[3] - MouseY; // invert mouse Y coordinate
1725 		MousePosZ = 1.0; // near clip plane depth
1726 
1727 		// Get unprojected end point
1728 		gluUnProject
1729 		(
1730 			MousePosX,
1731 			MousePosY,
1732 			MousePosZ,
1733 			ModelMatrix,
1734 			ProjMatrix,
1735 			Viewport,
1736 		    &WorldPosX,
1737 		    &WorldPosY,
1738 		    &WorldPosZ
1739 		);
1740 		WorldPos.x = WorldPosX;
1741 		WorldPos.y = WorldPosY;
1742 		WorldPos.z = WorldPosZ;
1743 
1744         CollisionPoint = line_plane_collision(posCamera, &WorldPos, quad);
1745 		Intersec[0]=CollisionPoint.x;
1746 		Intersec[1]=CollisionPoint.y;
1747 		Intersec[2]=CollisionPoint.z;
1748 
1749         if (CheckPointInQuad(&CollisionPoint, quad))
1750 			return true;
1751 
1752 	}
1753 
1754 
1755 	return false;
1756 }
1757 
1758 
1759 
1760 
1761 
1762 
1763 
1764 
1765 
1766 
1767 
1768 
1769 
1770 
1771 
1772 
1773 
1774 
1775