source: GTP/trunk/App/Demos/Vis/FriendlyCulling/src/Matrix4x4.cpp @ 2782

Revision 2782, 13.5 KB checked in by mattausch, 16 years ago (diff)
Line 
1#include "Matrix4x4.h"
2#include "Vector3.h"
3
4// standard headers
5#include <iomanip>
6using namespace std;
7
8namespace CHCDemoEngine
9{
10
11Matrix4x4::Matrix4x4()
12{}
13
14
15Matrix4x4::Matrix4x4(const Vector3 &a, const Vector3 &b, const Vector3 &c)
16{
17        // first index is column [x], the second is row [y]
18        x[0][0] = a.x;
19        x[1][0] = b.x;
20        x[2][0] = c.x;
21        x[3][0] = 0.0f;
22
23        x[0][1] = a.y;
24        x[1][1] = b.y;
25        x[2][1] = c.y;
26        x[3][1] = 0.0f;
27
28        x[0][2] = a.z;
29        x[1][2] = b.z;
30        x[2][2] = c.z;
31        x[3][2] = 0.0f;
32
33        x[0][3] = 0.0f;
34        x[1][3] = 0.0f;
35        x[2][3] = 0.0f;
36        x[3][3] = 1.0f;
37}
38
39
40void Matrix4x4::SetColumns(const Vector3 &a, const Vector3 &b, const Vector3 &c)
41{
42        // first index is column [x], the second is row [y]
43        x[0][0] = a.x;
44        x[1][0] = a.y;
45        x[2][0] = a.z;
46        x[3][0] = 0.0f;
47
48        x[0][1] = b.x;
49        x[1][1] = b.y;
50        x[2][1] = b.z;
51        x[3][1] = 0.0f;
52
53        x[0][2] = c.x;
54        x[1][2] = c.y;
55        x[2][2] = c.z;
56        x[3][2] = 0.0f;
57
58        x[0][3] = 0.0f;
59        x[1][3] = 0.0f;
60        x[2][3] = 0.0f;
61        x[3][3] = 1.0f;
62}
63
64// full constructor
65Matrix4x4::Matrix4x4(float x11, float x12, float x13, float x14,
66                                         float x21, float x22, float x23, float x24,
67                                         float x31, float x32, float x33, float x34,
68                                         float x41, float x42, float x43, float x44)
69{
70        // first index is column [x], the second is row [y]
71        x[0][0] = x11;
72        x[1][0] = x12;
73        x[2][0] = x13;
74        x[3][0] = x14;
75
76        x[0][1] = x21;
77        x[1][1] = x22;
78        x[2][1] = x23;
79        x[3][1] = x24;
80
81        x[0][2] = x31;
82        x[1][2] = x32;
83        x[2][2] = x33;
84        x[3][2] = x34;
85
86        x[0][3] = x41;
87        x[1][3] = x42;
88        x[2][3] = x43;
89        x[3][3] = x44;
90}
91
92
93float Matrix4x4::Det3x3() const
94{
95        return (x[0][0]*x[1][1]*x[2][2] + \
96                x[1][0]*x[2][1]*x[0][2] + \
97                x[2][0]*x[0][1]*x[1][2] - \
98                x[2][0]*x[1][1]*x[0][2] - \
99                x[0][0]*x[2][1]*x[1][2] - \
100                x[1][0]*x[0][1]*x[2][2]);
101}
102
103
104
105// inverse matrix computation gauss_jacobiho method .. from N.R. in C
106// if matrix is regular = computatation successfull = returns 0
107// in case of singular matrix returns 1
108int Matrix4x4::Invert()
109{
110        int indxc[4],indxr[4],ipiv[4];
111        int i,icol,irow,j,k,l,ll,n;
112        float big,dum,pivinv,temp;
113        // satisfy the compiler
114        icol = irow = 0;
115
116        // the size of the matrix
117        n = 4;
118
119        for ( j = 0 ; j < n ; j++) /* zero pivots */
120                ipiv[j] = 0;
121
122        for ( i = 0; i < n; i++)
123        {
124                big = 0.0;
125                for (j = 0 ; j < n ; j++)
126                        if (ipiv[j] != 1)
127                                for ( k = 0 ; k<n ; k++)
128                                {
129                                        if (ipiv[k] == 0)
130                                        {
131                                                if (fabs(x[k][j]) >= big)
132                                                {
133                                                        big = fabs(x[k][j]);
134                                                        irow = j;
135                                                        icol = k;
136                                                }
137                                        }
138                                        else
139                                                if (ipiv[k] > 1)
140                                                        return 1; /* singular matrix */
141                                }
142                                ++(ipiv[icol]);
143                                if (irow != icol)
144                                {
145                                        for ( l = 0 ; l<n ; l++)
146                                        {
147                                                temp = x[l][icol];
148                                                x[l][icol] = x[l][irow];
149                                                x[l][irow] = temp;
150                                        }
151                                }
152                                indxr[i] = irow;
153                                indxc[i] = icol;
154                                if (x[icol][icol] == 0.0)
155                                        return 1; /* singular matrix */
156
157                                pivinv = 1.0f / x[icol][icol];
158                                x[icol][icol] = 1.0f ;
159                                for ( l = 0 ; l<n ; l++)
160                                        x[l][icol] = x[l][icol] * pivinv ;
161
162                                for (ll = 0 ; ll < n ; ll++)
163                                        if (ll != icol)
164                                        {
165                                                dum = x[icol][ll];
166                                                x[icol][ll] = 0.0;
167                                                for ( l = 0 ; l<n ; l++)
168                                                        x[l][ll] = x[l][ll] - x[l][icol] * dum ;
169                                        }
170        }
171        for ( l = n; l--; )
172        {
173                if (indxr[l] != indxc[l])
174                        for ( k = 0; k<n ; k++)
175                        {
176                                temp = x[indxr[l]][k];
177                                x[indxr[l]][k] = x[indxc[l]][k];
178                                x[indxc[l]][k] = temp;
179                        }
180        }
181
182        return 0 ; // matrix is regular .. inversion has been succesfull
183}
184
185
186// Invert the given matrix using the above inversion routine.
187Matrix4x4 Invert(const Matrix4x4& M)
188{
189        Matrix4x4 InvertMe = M;
190        InvertMe.Invert();
191        return InvertMe;
192}
193
194
195// Transpose the matrix.
196void Matrix4x4::Transpose()
197{
198        for (int i = 0; i < 4; i++)
199                for (int j = i; j < 4; j++)
200                        if (i != j) {
201                                float temp = x[i][j];
202                                x[i][j] = x[j][i];
203                                x[j][i] = temp;
204                        }
205}
206
207
208// Transpose the given matrix using the transpose routine above.
209Matrix4x4 Transpose(const Matrix4x4& M)
210{
211        Matrix4x4 TransposeMe = M;
212        TransposeMe.Transpose();
213        return TransposeMe;
214}
215
216
217// Construct an identity matrix.
218Matrix4x4 IdentityMatrix()
219{
220        Matrix4x4 M;
221
222        for (int i = 0; i < 4; i++)
223                for (int j = 0; j < 4; j++)
224                        M.x[i][j] = (i == j) ? 1.0f : 0.0f;
225        return M;
226}
227
228
229// Construct a zero matrix.
230Matrix4x4 ZeroMatrix()
231{
232        Matrix4x4 M;
233        for (int i = 0; i < 4; i++)
234                for (int j = 0; j < 4; j++)
235                        M.x[i][j] = 0;
236        return M;
237}
238
239// Construct a translation matrix given the location to translate to.
240Matrix4x4
241TranslationMatrix(const Vector3& Location)
242{
243        Matrix4x4 M = IdentityMatrix();
244        M.x[3][0] = Location.x;
245        M.x[3][1] = Location.y;
246        M.x[3][2] = Location.z;
247        return M;
248}
249
250// Construct a rotation matrix.  Rotates Angle radians about the
251// X axis.
252Matrix4x4
253RotationXMatrix(float Angle)
254{
255        Matrix4x4 M = IdentityMatrix();
256        float Cosine = cos(Angle);
257        float Sine = sin(Angle);
258        M.x[1][1] = Cosine;
259        M.x[2][1] = -Sine;
260        M.x[1][2] = Sine;
261        M.x[2][2] = Cosine;
262        return M;
263}
264
265// Construct a rotation matrix.  Rotates Angle radians about the
266// Y axis.
267Matrix4x4 RotationYMatrix(float angle)
268{
269        Matrix4x4 m = IdentityMatrix();
270
271        float cosine = cos(angle);
272        float sine = sin(angle);
273
274        m.x[0][0] = cosine;
275        m.x[2][0] = -sine;
276        m.x[0][2] = sine;
277        m.x[2][2] = cosine;
278
279        return m;
280}
281
282
283// Construct a rotation matrix.  Rotates Angle radians about the
284// Z axis.
285Matrix4x4 RotationZMatrix(float angle)
286{
287        Matrix4x4 m = IdentityMatrix();
288
289        float cosine = cos(angle);
290        float sine = sin(angle);
291       
292        m.x[0][0] = cosine;
293        m.x[1][0] = -sine;
294        m.x[0][1] = sine;
295        m.x[1][1] = cosine;
296
297        return m;
298}
299
300
301// Construct a yaw-pitch-roll rotation matrix.  Rotate Yaw
302// radians about the XY axis, rotate Pitch radians in the
303// plane defined by the Yaw rotation, and rotate Roll radians
304// about the axis defined by the previous two angles.
305Matrix4x4 RotationYPRMatrix(float Yaw, float Pitch, float Roll)
306{
307        Matrix4x4 M;
308        float ch = cos(Yaw);
309        float sh = sin(Yaw);
310        float cp = cos(Pitch);
311        float sp = sin(Pitch);
312        float cr = cos(Roll);
313        float sr = sin(Roll);
314
315        M.x[0][0] = ch * cr + sh * sp * sr;
316        M.x[1][0] = -ch * sr + sh * sp * cr;
317        M.x[2][0] = sh * cp;
318        M.x[0][1] = sr * cp;
319        M.x[1][1] = cr * cp;
320        M.x[2][1] = -sp;
321        M.x[0][2] = -sh * cr - ch * sp * sr;
322        M.x[1][2] = sr * sh + ch * sp * cr;
323        M.x[2][2] = ch * cp;
324        for (int i = 0; i < 4; i++)
325                M.x[3][i] = M.x[i][3] = 0;
326        M.x[3][3] = 1;
327
328        return M;
329}
330
331// Construct a rotation of a given angle about a given axis.
332// Derived from Eric Haines's SPD (Standard Procedural
333// Database).
334Matrix4x4 RotationAxisMatrix(const Vector3& axis, float angle)
335{
336        Matrix4x4 M;
337
338        float cosine = cos(angle);
339        float sine = sin(angle);
340        float one_minus_cosine = 1 - cosine;
341
342        M.x[0][0] = axis.x * axis.x + (1.0f - axis.x * axis.x) * cosine;
343        M.x[0][1] = axis.x * axis.y * one_minus_cosine + axis.z * sine;
344        M.x[0][2] = axis.x * axis.z * one_minus_cosine - axis.y * sine;
345        M.x[0][3] = 0;
346
347        M.x[1][0] = axis.x * axis.y * one_minus_cosine - axis.z * sine;
348        M.x[1][1] = axis.y * axis.y + (1.0f - axis.y * axis.y) * cosine;
349        M.x[1][2] = axis.y * axis.z * one_minus_cosine + axis.x * sine;
350        M.x[1][3] = 0;
351
352        M.x[2][0] = axis.x * axis.z * one_minus_cosine + axis.y * sine;
353        M.x[2][1] = axis.y * axis.z * one_minus_cosine - axis.x * sine;
354        M.x[2][2] = axis.z * axis.z + (1.0f - axis.z * axis.z) * cosine;
355        M.x[2][3] = 0;
356
357        M.x[3][0] = 0;
358        M.x[3][1] = 0;
359        M.x[3][2] = 0;
360        M.x[3][3] = 1;
361
362        return M;
363}
364
365
366// Constructs the rotation matrix that rotates 'vec1' to 'vec2'
367Matrix4x4 RotationVectorsMatrix(const Vector3 &vecStart, const Vector3 &vecTo)
368{
369        Vector3 vec = CrossProd(vecStart, vecTo);
370
371        if (Magnitude(vec) > Limits::Small) {
372                // vector exist, compute angle
373                float angle = acos(DotProd(vecStart, vecTo));
374                // normalize for sure
375                vec.Normalize();
376                return RotationAxisMatrix(vec, angle);
377        }
378
379        // opposite or colinear vectors
380        Matrix4x4 ret = IdentityMatrix();
381        if (DotProd(vecStart, vecTo) < 0.0f)
382                ret *= -1.0f; // opposite vectors
383
384        return ret;
385}
386
387
388// Construct a scale matrix given the X, Y, and Z parameters
389// to scale by.  To scale uniformly, let X==Y==Z.
390Matrix4x4 ScaleMatrix(float X, float Y, float Z)
391{
392        Matrix4x4 M = IdentityMatrix();
393
394        M.x[0][0] = X;
395        M.x[1][1] = Y;
396        M.x[2][2] = Z;
397
398        return M;
399}
400
401
402// Construct a rotation matrix that makes the x, y, z axes
403// correspond to the vectors given.
404Matrix4x4 GenRotation(const Vector3 &x, const Vector3 &y, const Vector3 &z)
405{
406        Matrix4x4 M = IdentityMatrix();
407
408        //  x  y
409        M.x[0][0] = x.x;
410        M.x[1][0] = x.y;
411        M.x[2][0] = x.z;
412
413        M.x[0][1] = y.x;
414        M.x[1][1] = y.y;
415        M.x[2][1] = y.z;
416
417        M.x[0][2] = z.x;
418        M.x[1][2] = z.y;
419        M.x[2][2] = z.z;
420
421        return M;
422}
423
424// Construct a quadric matrix.  After Foley et al. pp. 528-529.
425Matrix4x4
426QuadricMatrix(float a, float b, float c, float d, float e,
427                          float f, float g, float h, float j, float k)
428{
429        Matrix4x4 M;
430
431        M.x[0][0] = a;  M.x[0][1] = d;  M.x[0][2] = f;  M.x[0][3] = g;
432        M.x[1][0] = d;  M.x[1][1] = b;  M.x[1][2] = e;  M.x[1][3] = h;
433        M.x[2][0] = f;  M.x[2][1] = e;  M.x[2][2] = c;  M.x[2][3] = j;
434        M.x[3][0] = g;  M.x[3][1] = h;  M.x[3][2] = j;  M.x[3][3] = k;
435
436        return M;
437}
438
439// Construct various "mirror" matrices, which flip coordinate
440// signs in the various axes specified.
441Matrix4x4
442MirrorX()
443{
444        Matrix4x4 M = IdentityMatrix();
445        M.x[0][0] = -1;
446        return M;
447}
448
449Matrix4x4
450MirrorY()
451{
452        Matrix4x4 M = IdentityMatrix();
453        M.x[1][1] = -1;
454        return M;
455}
456
457Matrix4x4
458MirrorZ()
459{
460        Matrix4x4 M = IdentityMatrix();
461        M.x[2][2] = -1;
462        return M;
463}
464
465Matrix4x4
466RotationOnly(const Matrix4x4& x)
467{
468        Matrix4x4 M = x;
469        M.x[3][0] = M.x[3][1] = M.x[3][2] = 0;
470        return M;
471}
472
473// Add corresponding elements of the two matrices.
474Matrix4x4&
475Matrix4x4::operator+= (const Matrix4x4& A)
476{
477        for (int i = 0; i < 4; i++)
478                for (int j = 0; j < 4; j++)
479                        x[i][j] += A.x[i][j];
480        return *this;
481}
482
483// Subtract corresponding elements of the matrices.
484Matrix4x4&
485Matrix4x4::operator-= (const Matrix4x4& A)
486{
487        for (int i = 0; i < 4; i++)
488                for (int j = 0; j < 4; j++)
489                        x[i][j] -= A.x[i][j];
490        return *this;
491}
492
493// Scale each element of the matrix by A.
494Matrix4x4&
495Matrix4x4::operator*= (float A)
496{
497        for (int i = 0; i < 4; i++)
498                for (int j = 0; j < 4; j++)
499                        x[i][j] *= A;
500        return *this;
501}
502
503// Multiply two matrices.
504Matrix4x4&
505Matrix4x4::operator*= (const Matrix4x4& A)
506{
507        Matrix4x4 ret = *this;
508
509        for (int i = 0; i < 4; i++)
510                for (int j = 0; j < 4; j++) {
511                        float subt = 0;
512                        for (int k = 0; k < 4; k++)
513                                subt += ret.x[i][k] * A.x[k][j];
514                        x[i][j] = subt;
515                }
516                return *this;
517}
518
519// Add corresponding elements of the matrices.
520Matrix4x4
521operator+ (const Matrix4x4& A, const Matrix4x4& B)
522{
523        Matrix4x4 ret;
524
525        for (int i = 0; i < 4; i++)
526                for (int j = 0; j < 4; j++)
527                        ret.x[i][j] = A.x[i][j] + B.x[i][j];
528        return ret;
529}
530
531// Subtract corresponding elements of the matrices.
532Matrix4x4
533operator- (const Matrix4x4& A, const Matrix4x4& B)
534{
535        Matrix4x4 ret;
536
537        for (int i = 0; i < 4; i++)
538                for (int j = 0; j < 4; j++)
539                        ret.x[i][j] = A.x[i][j] - B.x[i][j];
540        return ret;
541}
542
543// Multiply matrices.
544Matrix4x4
545operator* (const Matrix4x4& A, const Matrix4x4& B)
546{
547        Matrix4x4 ret;
548
549        for (int i = 0; i < 4; i++)
550                for (int j = 0; j < 4; j++) {
551                        float subt = 0;
552                        for (int k = 0; k < 4; k++)
553                                subt += A.x[i][k] * B.x[k][j];
554                        ret.x[i][j] = subt;
555                }
556                return ret;
557}
558
559// Transform a vector by a matrix.
560Vector3
561operator* (const Matrix4x4& M, const Vector3& v)
562{
563        Vector3 ret;
564        float denom;
565
566        ret.x = v.x * M.x[0][0] + v.y * M.x[1][0] + v.z * M.x[2][0] + M.x[3][0];
567        ret.y = v.x * M.x[0][1] + v.y * M.x[1][1] + v.z * M.x[2][1] + M.x[3][1];
568        ret.z = v.x * M.x[0][2] + v.y * M.x[1][2] + v.z * M.x[2][2] + M.x[3][2];
569        denom = M.x[0][3] + M.x[1][3] + M.x[2][3] + M.x[3][3];
570        if (denom != 1.0)
571                ret /= denom;
572        return ret;
573}
574
575// Apply the rotation portion of a matrix to a vector.
576Vector3
577RotateOnly(const Matrix4x4& M, const Vector3& v)
578{
579        Vector3 ret;
580        float denom;
581
582        ret.x = v.x * M.x[0][0] + v.y * M.x[1][0] + v.z * M.x[2][0];
583        ret.y = v.x * M.x[0][1] + v.y * M.x[1][1] + v.z * M.x[2][1];
584        ret.z = v.x * M.x[0][2] + v.y * M.x[1][2] + v.z * M.x[2][2];
585        denom = M.x[0][3] + M.x[1][3] + M.x[2][3] + M.x[3][3];
586        if (denom != 1.0)
587                ret /= denom;
588        return ret;
589}
590
591// Scale each element of the matrix by B.
592Matrix4x4
593operator* (const Matrix4x4& A, float B)
594{
595        Matrix4x4 ret;
596
597        for (int i = 0; i < 4; i++)
598                for (int j = 0; j < 4; j++)
599                        ret.x[i][j] = A.x[i][j] * B;
600        return ret;
601}
602
603// Overloaded << for C++-style output.
604ostream&
605operator<< (ostream& s, const Matrix4x4& M)
606{
607        for (int i = 0; i < 4; i++) { // y
608                for (int j = 0; j < 4; j++) { // x
609                        //  x  y
610                        s << setprecision(4) << setw(10) << M.x[j][i];
611                }
612                s << '\n';
613        }
614        return s;
615}
616
617
618Vector3 PlaneRotate(const Matrix4x4& tform, const Vector3& p)
619{
620        // I sure hope that matrix is invertible...
621        Matrix4x4 use = Transpose(Invert(tform));
622
623        return RotateOnly(use, p);
624}
625
626// Transform a normal
627Vector3 TransformNormal(const Matrix4x4& tform, const Vector3& n)
628{
629        Matrix4x4 use = NormalTransformMatrix(tform);
630
631        return RotateOnly(use, n);
632}
633
634Matrix4x4
635NormalTransformMatrix(const Matrix4x4 &tform)
636{
637        Matrix4x4 m = tform;
638        // for normal translation vector must be zero!
639        m.x[3][0] = m.x[3][1] = m.x[3][2] = 0.0;
640        // I sure hope that matrix is invertible...
641        return Transpose(Invert(m));
642}
643
644
645Vector3 GetTranslation(const Matrix4x4 &M)
646{
647        return Vector3(M.x[3][0], M.x[3][1], M.x[3][2]);
648}
649
650}
Note: See TracBrowser for help on using the repository browser.