Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: downloads/ogre/Tools/3dsmaxExport/LEXIExporter/SharedUtilities/Sources/MathMatrix4x4.cpp @ 45

Last change on this file since 45 was 6, checked in by anonymous, 18 years ago

=…
File size: 16.9 KB
Line 
1/*
2-----------------------------------------------------------------------------
3This source file is part of LEXIExporter
4
5Copyright 2006 NDS Limited
6
7Author(s):
8Bo Krohn
9
10This program is free software; you can redistribute it and/or modify it under
11the terms of the GNU Lesser General Public License as published by the Free Software
12Foundation; either version 2 of the License, or (at your option) any later
13version.
14
15This program is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
18
19You should have received a copy of the GNU Lesser General Public License along with
20this program; if not, write to the Free Software Foundation, Inc., 59 Temple
21Place - Suite 330, Boston, MA 02111-1307, USA, or go to
22http://www.gnu.org/copyleft/lesser.txt.
23-----------------------------------------------------------------------------
24*/
25
26/////////////////////////////////////////////////////
27//
28//  Matrix 4x4 class
29//
30/////////////////////////////////////////////////////
31
32#include "stdafx.h"
33#include "MathMatrix4x4.h"
34
35//
36
37static float mZero[16] = {              0.0f, 0.0f, 0.0f, 0.0f,
38                                                                0.0f, 0.0f, 0.0f, 0.0f,
39                                                                0.0f, 0.0f, 0.0f, 0.0f,
40                                                                0.0f, 0.0f, 0.0f, 0.0f };
41
42const CMatrix CMatrix::_zero(mZero);
43
44//
45
46static float mIdentity[16] = {  1.0f, 0.0f, 0.0f, 0.0f,
47                                                                0.0f, 1.0f, 0.0f, 0.0f,
48                                                                0.0f, 0.0f, 1.0f, 0.0f,
49                                                                0.0f, 0.0f, 0.0f, 1.0f };
50
51const CMatrix CMatrix::_identity(mIdentity);
52
53//
54
55inline float det3x3(    float a1,float a2,float a3,
56                                                float b1,float b2,float b3,
57                                                float c1,float c2,float c3) 
58{ 
59        return a1*(b2*c3-b3*c2)-b1*(a2*c3-a3*c2)+c1*(a2*b3-a3*b2);
60}
61
62//
63
64void CMatrix::makeLookAt(const CVec3& eye, const CVec3& point, const CVec3& up)
65{
66        CVec3 f;
67        f.subtract(eye, point);                 // view vector  (maps to z)
68        const float flen = f.length2();
69        if(F_Min < flen)
70                f.scale( 1.0f / sqrtf(flen) );
71       
72        CVec3 upprime = up;
73        const float ulen = upprime.length2();
74        if(F_Min < ulen)
75                upprime.scale( 1.0f / sqrtf(ulen) );
76       
77        CVec3 s; 
78        s.cross(upprime, f);                    // s = up X f  (maps to x)
79        const float slen = s.length2();
80        if(F_Min < slen)
81                s.scale( 1.0f / sqrtf(slen) );
82
83        CVec3 u;
84        u.cross(f, s);                          // u = f X s;   (maps to y)
85        // s and f are normalized and orthogonal, so u is
86               
87        // this matrix maps to the eye point, we want to map the geometry so we
88        // need the inverse. since it's an orthonormal matrix by construction,
89        // we can simply transpose it.
90        // [1 0 0 0] [[  s ]]
91        // [0 1 0 0] [[  u ]]
92        // [0 0 1 0] [[  f ]]
93        // [0 0 0 1]
94       
95        mat[0][0] =  s.x; mat[0][1] =  u.x; mat[0][2] = f.x; mat[0][3] =  0.0f;
96        mat[1][0] =  s.y; mat[1][1] =  u.y; mat[1][2] = f.y; mat[1][3] =  0.0f;
97        mat[2][0] =  s.z; mat[2][1] =  u.z; mat[2][2] = f.z; mat[2][3] =  0.0f;
98
99        // translate eye to origin
100        mat[3][0] = -( mat[0][0] * eye.x + mat[1][0] * eye.y + mat[2][0] * eye.z );
101        mat[3][1] = -( mat[0][1] * eye.x + mat[1][1] * eye.y + mat[2][1] * eye.z );
102        mat[3][2] = -( mat[0][2] * eye.x + mat[1][2] * eye.y + mat[2][2] * eye.z );
103        mat[3][3] = 1.0f - ( mat[0][3] * eye.x + mat[1][3] * eye.y + mat[2][3] * eye.z );
104}
105
106void CMatrix::makeLookAtDirection(const CVec3& eye, const CVec3& dir, const CVec3& up)
107{
108        CVec3 f;
109        f.negate(dir);                           // view vector  (maps to z)
110        const float flen = f.length2();
111        if(F_Min < flen)
112                f.scale( 1.0f / sqrtf(flen) );
113
114        CVec3 upprime = up;
115        const float ulen = upprime.length2();
116        if(F_Min < ulen)
117                upprime.scale( 1.0f / sqrtf(ulen) );
118
119        CVec3 s; 
120        s.cross(upprime, f);                    // s = up X f  (maps to x)
121        const float slen = s.length2();
122        if(F_Min < slen)
123                s.scale( 1.0f / sqrtf(slen) );
124
125        CVec3 u;
126        u.cross(f, s);                          // u = f X s;   (maps to y)
127        // s and f are normalized and orthogonal, so u is
128               
129        // this matrix maps to the eye point, we want to map the geometry so we
130        // need the inverse. since it's an orthonormal matrix by construction,
131        // we can simply transpose it.
132        // [1 0 0 0] [[  s ]]
133        // [0 1 0 0] [[  u ]]
134        // [0 0 1 0] [[  f ]]
135        // [0 0 0 1]
136       
137        mat[0][0] =  s.x; mat[0][1] =  u.x; mat[0][2] = f.x; mat[0][3] =  0.0f;
138        mat[1][0] =  s.y; mat[1][1] =  u.y; mat[1][2] = f.y; mat[1][3] =  0.0f;
139        mat[2][0] =  s.z; mat[2][1] =  u.z; mat[2][2] = f.z; mat[2][3] =  0.0f;
140
141        // translate eye to origin
142        mat[3][0] = -( mat[0][0] * eye.x + mat[1][0] * eye.y + mat[2][0] * eye.z );
143        mat[3][1] = -( mat[0][1] * eye.x + mat[1][1] * eye.y + mat[2][1] * eye.z );
144        mat[3][2] = -( mat[0][2] * eye.x + mat[1][2] * eye.y + mat[2][2] * eye.z );
145        mat[3][3] = 1.0f - ( mat[0][3] * eye.x + mat[1][3] * eye.y + mat[2][3] * eye.z );
146}
147
148//
149
150void CMatrix::makePerspective(float left, float right, float bottom, float top, float znear, float zfar)
151{
152        float temp = 1.0f / ( right - left );
153
154        mat[0][0]= ( 2.0f * znear ) * temp;
155        mat[1][0]= 0.0f;       
156        mat[2][0]= ( right + left ) * temp;     // for asymmetric views
157        mat[3][0]= 0.0f;
158       
159        temp = 1.0f / ( top - bottom );
160
161        mat[0][1]= 0.0f;
162        mat[1][1]= ( 2.0f * znear ) * temp;
163        mat[2][1]= ( top + bottom ) * temp;     // for asymmetric views
164        mat[3][1]= 0.0f;
165       
166        temp = 1.0f / ( zfar - znear );
167
168        mat[0][2]= 0.0f;
169        mat[1][2]= 0.0f;
170        mat[2][2]= -( zfar + znear ) * temp;
171        mat[3][2]=      ( -2.0f * zfar * znear ) * temp;
172
173        mat[0][3]= 0.0f;
174        mat[1][3]= 0.0f;
175        mat[2][3]= -1.0f;
176        mat[3][3]= 0.0f;
177}
178
179void CMatrix::makePerspectiveFOV(float hfov, float vfov, float aspect, float znear, float zfar)
180{
181        float hfovRad=UtilDegToRad(hfov);
182        float vfovRad=UtilDegToRad(vfov);
183
184        float l,r,t,b;
185        float h,v;
186
187        if( hfovRad < 0.0f ) 
188        {
189                h = (znear * tanf( vfovRad * 0.5f )) * aspect;
190                h = atanf( h / znear ) * 2.0f;
191                v = vfovRad;
192        }
193        else 
194        {
195                v = (znear * tanf( hfovRad * 0.5f )) / aspect;
196                v = atanf( v / znear ) * 2.0f;
197                h = hfovRad;
198        }
199
200        l  = znear * -tanf( h * 0.5f );
201        r  = - l;
202        b  = znear * -tanf( v * 0.5f );
203        t  = - b;
204
205        makePerspective( l, r, b, t, znear, zfar );
206}
207
208void CMatrix::makeOrthogonalPerspective(float left, float right, float bottom, float top, float znear, float zfar)
209{
210        float temp = 1.0f / ( right - left );
211
212        mat[0][0]= 2.0f * temp;
213        mat[1][0]= 0.0f;
214        mat[2][0]= 0.0f;
215        mat[3][0]= - ( right + left ) * temp;
216
217        temp = 1.0f / ( top - bottom );
218       
219        mat[0][1]= 0.0f;
220        mat[1][1]= 2.0f * temp;
221        mat[2][1]= 0.0f;
222        mat[3][1]= -( top + bottom ) * temp;
223
224        temp = 1.0f / ( zfar - znear );
225       
226        mat[0][2]= 0.0f;
227        mat[1][2]= 0.0f;
228        mat[2][2]= -2.0f * temp;
229        mat[3][2]= -( zfar + znear ) * temp;
230       
231        mat[0][3]= 0.0f;
232        mat[1][3]= 0.0f;
233        mat[2][3]= 0.0f;
234        mat[3][3]= 1.0f;
235}
236
237//
238
239void CMatrix::setRotationRadians(float angle, const CVec3& axis)
240{
241        if(fabs(angle) < 0.0000005f) 
242        {
243                mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
244                mat[0][1] = mat[0][2] = mat[1][0] = mat[1][2] = mat[2][0] = mat[2][1] = 0.f;
245        } 
246        else 
247        {
248                float sine=sinf(angle);
249                float cosine=cosf(angle);
250
251                CVec3 sineAxis;
252                sineAxis.scale(sine, axis);
253
254                float t = 1.0f - cosine;
255                float tx = t * axis.x;
256                mat[0][0] = tx * axis.x + cosine;
257                mat[0][1] = tx * axis.y + sineAxis.z;
258                mat[0][2] = tx * axis.z - sineAxis.y;
259
260                float ty = t * axis.y;
261                mat[1][0] = ty * axis.x - sineAxis.z;
262                mat[1][1] = ty * axis.y + cosine;
263                mat[1][2] = ty * axis.z + sineAxis.x;
264
265                float tz = t * axis.z;
266                mat[2][0] = tz * axis.x + sineAxis.y;
267                mat[2][1] = tz * axis.y - sineAxis.x;
268                mat[2][2] = tz * axis.z + cosine;
269        }
270}
271
272void CMatrix::setRotationRadians(float anglex, float angley, float anglez)
273{
274        float sx, sy, sz, cx, cy, cz;
275        float sxsy, cxsz, cxcz;
276
277        if( anglex != 0.0f ) { sincos( anglex, sx, cx ); }
278        else { sx = 0.0f; cx = 1.0f; sxsy = 0.0f; }
279
280        if( angley != 0.0f ) { sincos( angley, sy, cy ); sxsy = sx * sy; }
281        else { sy = 0.0f; cy = 1.0f; sxsy = 0.0f; }
282
283        if( anglez != 0.0f ) { sincos( anglez, sz, cz ); cxsz = cx * sz; }
284        else { sz = 0.0f; cz = 1.0f; cxsz = 0.0f; }
285
286        cxcz = cx * cz;
287
288        mat[0][0] = cy * cz;
289        mat[0][1] = cy * sz;
290        mat[0][2] = -sy;
291
292        mat[1][0] = sxsy * cz - cxsz;
293        mat[1][1] = sxsy * sz + cx * cz;
294        mat[1][2] = sx * cy;
295
296        mat[2][0] = cxcz * sy + sx * sz;
297        mat[2][1] = cxsz * sy - sx * cz;
298        mat[2][2] = cx * cy;
299}
300
301void CMatrix::getRotationRadians(float& anglex, float& angley, float& anglez)
302{
303        CVec3 temp;
304        CVec3 row0( (float*)mat[0] );
305        CVec3 row1( (float*)mat[1] );
306        CVec3 row2( (float*)mat[2] );
307
308        if( mat[3][3] != 1.0f )
309        {
310                float global_scale_inverse = 1.0f / mat[3][3];
311                row0.scale( global_scale_inverse );
312                row1.scale( global_scale_inverse );
313                row2.scale( global_scale_inverse );
314        }
315
316        // possible scale or shearing must be removed...
317        row0.normalize();
318
319        // Compute XY shear factor and make 2nd row orthogonal to 1st.
320        float shearXY = row0.dot( row1 );
321        row1.addScaled( -shearXY, row0 );
322
323        // Now, normalize 2nd row.
324        row1.normalize();
325
326        // Compute XZ and YZ shears, orthogonalize 3rd row.
327        float shearXZ = row0.dot( row2 );
328        row2.addScaled( -shearXZ, row0 );
329        float shearYZ = row1.dot( row2 );
330        row2.addScaled( -shearYZ, row1 );
331
332        // Next, normalize 3rd row.
333        row2.normalize();
334
335        // Check for a coordinate system flip.  If the determinant is -1, then negate the rows.
336        temp.cross( row1, row2 );
337        if( row0.dot( temp ) < 0.0f )
338        {
339                row0.negate();
340                row1.negate();
341                row2.negate();
342        }
343
344        angley = asin( -row0.z );
345        if( cosf( angley ) != 0.0f ) 
346        {
347                anglex = atan2f( row1.z, row2.z );
348                anglez = atan2f( row0.y, row0.x );
349        } 
350        else 
351        {
352                anglex = atan2f( row1.x, row1.y );
353                anglez = 0.0f;
354        }
355}
356
357//
358
359void CMatrix::multiply(const CMatrix& m)
360{
361        /*register */unsigned int i;
362        float m2[4][4];
363
364        for( i=0; i<4; i++ ) {
365
366                m2[0][i] = ( m.mat[0][0] * mat[0][i] +
367                                                  m.mat[0][1] * mat[1][i] +
368                                                  m.mat[0][2] * mat[2][i] +
369                                                  m.mat[0][3] * mat[3][i] );
370
371                m2[1][i] = ( m.mat[1][0] * mat[0][i] +
372                                                  m.mat[1][1] * mat[1][i] +
373                                                  m.mat[1][2] * mat[2][i] +
374                                                  m.mat[1][3] * mat[3][i] );
375
376                m2[2][i] = ( m.mat[2][0] * mat[0][i] +
377                                                  m.mat[2][1] * mat[1][i] +
378                                                  m.mat[2][2] * mat[2][i] +
379                                                  m.mat[2][3] * mat[3][i] );
380
381                m2[3][i] = ( m.mat[3][0] * mat[0][i] +
382                                                  m.mat[3][1] * mat[1][i] +
383                                                  m.mat[3][2] * mat[2][i] +
384                                                  m.mat[3][3] * mat[3][i] );
385        }
386
387        memcpy(mat, m2, 16*sizeof(float));
388}
389
390void CMatrix::multiply(const CMatrix& m1, const CMatrix& m2)
391{
392        register unsigned int i;
393
394        if( this != &m1 && this != &m2 ) {
395
396                for( i=0; i<4; i++ ) {
397
398                        mat[0][i] = ( m1.mat[0][0] * m2.mat[0][i] +
399                                                  m1.mat[0][1] * m2.mat[1][i] +
400                                                  m1.mat[0][2] * m2.mat[2][i] +
401                                                  m1.mat[0][3] * m2.mat[3][i] );
402
403                        mat[1][i] = ( m1.mat[1][0] * m2.mat[0][i] +
404                                                  m1.mat[1][1] * m2.mat[1][i] +
405                                                  m1.mat[1][2] * m2.mat[2][i] +
406                                                  m1.mat[1][3] * m2.mat[3][i] );
407
408                        mat[2][i] = ( m1.mat[2][0] * m2.mat[0][i] +
409                                                  m1.mat[2][1] * m2.mat[1][i] +
410                                                  m1.mat[2][2] * m2.mat[2][i] +
411                                                  m1.mat[2][3] * m2.mat[3][i] );
412
413                        mat[3][i] = ( m1.mat[3][0] * m2.mat[0][i] +
414                                                  m1.mat[3][1] * m2.mat[1][i] +
415                                                  m1.mat[3][2] * m2.mat[2][i] +
416                                                  m1.mat[3][3] * m2.mat[3][i] );
417                }
418
419        }
420        else {
421               
422                float m3[4][4];
423
424                for( i=0; i<4; i++ ) {
425
426                        m3[0][i] = (    m1.mat[0][0] * m2.mat[0][i] +
427                                                        m1.mat[0][1] * m2.mat[1][i] +
428                                                        m1.mat[0][2] * m2.mat[2][i] +
429                                                        m1.mat[0][3] * m2.mat[3][i] );
430
431                        m3[1][i] = (    m1.mat[1][0] * m2.mat[0][i] +
432                                                        m1.mat[1][1] * m2.mat[1][i] +
433                                                        m1.mat[1][2] * m2.mat[2][i] +
434                                                        m1.mat[1][3] * m2.mat[3][i] );
435
436                        m3[2][i] = ( m1.mat[2][0] * m2.mat[0][i] +
437                                                          m1.mat[2][1] * m2.mat[1][i] +
438                                                          m1.mat[2][2] * m2.mat[2][i] +
439                                                          m1.mat[2][3] * m2.mat[3][i] );
440
441                        m3[3][i] = ( m1.mat[3][0] * m2.mat[0][i] +
442                                                          m1.mat[3][1] * m2.mat[1][i] +
443                                                          m1.mat[3][2] * m2.mat[2][i] +
444                                                          m1.mat[3][3] * m2.mat[3][i] );
445                }
446
447                memcpy(mat, m3, 16*sizeof(float));
448        }
449}
450
451//
452
453void CMatrix::transformPoints(const CVec3* from, CVec3* to, unsigned int iCount) const
454{
455        register float t0, t1, t2, w;
456
457        for( unsigned int i = 0; i < iCount; i++ )
458        {
459                t0 = from[i].x;
460                t1 = from[i].y;
461                t2 = from[i].z;
462
463                to[i].x = (t0 * mat[0][0] + t1 * mat[1][0] + t2 * mat[2][0] + mat[3][0]);
464                to[i].y = (t0 * mat[0][1] + t1 * mat[1][1] + t2 * mat[2][1] + mat[3][1]);
465                to[i].z = (t0 * mat[0][2] + t1 * mat[1][2] + t2 * mat[2][2] + mat[3][2]);
466                w = (t0 * mat[0][3] + t1 * mat[1][3] + t2 * mat[2][3] + mat[3][3]);
467
468                if( w != 1.0f ) {
469                        if( fabs( w ) < F_MinValue )
470                                w = F_MinValue;
471                        w = 1.0f / w;
472                        to[i].x *= w;
473                        to[i].y *= w;
474                        to[i].z *= w;
475                }
476        }
477}
478
479void CMatrix::transformPoints(const CVec4* from, CVec4* to, unsigned int iCount) const
480{
481        register float t0, t1, t2, t3;
482
483        for( unsigned int i = 0; i < iCount; i++ )
484        {
485                t0 = from[i].x;
486                t1 = from[i].y;
487                t2 = from[i].z;
488                t3 = from[i].w;
489
490                to[i].x = (t0 * mat[0][0] + t1 * mat[1][0] + t2 * mat[2][0] + t3 * mat[3][0]);
491                to[i].y = (t0 * mat[0][1] + t1 * mat[1][1] + t2 * mat[2][1] + t3 * mat[3][1]);
492                to[i].z = (t0 * mat[0][2] + t1 * mat[1][2] + t2 * mat[2][2] + t3 * mat[3][2]);
493                to[i].w = (t0 * mat[0][3] + t1 * mat[1][3] + t2 * mat[2][3] + t3 * mat[3][3]);
494        }
495}
496
497void CMatrix::transformVectors(const CVec3* from, CVec3* to, unsigned int iCount) const
498{
499        /*register */float t0, t1, t2;
500
501        for(unsigned int i = 0; i < iCount; i++, from++, to++)
502        {
503                t0 = from->x;
504                t1 = from->y;
505                t2 = from->z;
506
507                to->x = (t0 * mat[0][0] + t1 * mat[1][0] + t2 * mat[2][0]);
508                to->y = (t0 * mat[0][1] + t1 * mat[1][1] + t2 * mat[2][1]);
509                to->z = (t0 * mat[0][2] + t1 * mat[1][2] + t2 * mat[2][2]);
510        }
511}
512
513void CMatrix::transformVectors(const CVec4* from, CVec4* to, unsigned int iCount) const
514{
515        /*register */float t0, t1, t2;
516
517        for(unsigned int i = 0; i < iCount; i++, from++, to++)
518        {
519                t0 = from->x;
520                t1 = from->y;
521                t2 = from->z;
522
523                to->x = (t0 * mat[0][0] + t1 * mat[1][0] + t2 * mat[2][0]);
524                to->y = (t0 * mat[0][1] + t1 * mat[1][1] + t2 * mat[2][1]);
525                to->z = (t0 * mat[0][2] + t1 * mat[1][2] + t2 * mat[2][2]);
526                to->w = from->w;
527        }
528}
529
530//
531
532void CMatrix::invert()
533{
534        float det, idet;
535        CMatrix local_matrix;
536
537        const CMatrix& matrix=*this;
538
539        // calculate the adjoint matrix
540        adjoint( matrix, local_matrix );
541        // calculate the 4x4 determinant if the determinant is zero,
542        // then the inverse matrix is not unique.
543        det = matrix.determinant();
544       
545        // This test is only made to avoid crash
546        //   it is not a test of matrix inversibility
547        if( fabs( det ) < F_Min )
548                throw;
549
550        // scale the adjoint matrix to get the inverse
551        idet = 1.0f / det;
552        for(unsigned int i=0; i<4; ++i)
553                for(unsigned int j=0; j<4; ++j)
554                        mat[i][j] = local_matrix.mat[i][j] * idet;
555}
556
557//
558
559void CMatrix::transpose()
560{
561        float m[4][4];
562
563        m[0][0] = mat[0][0];
564        m[0][1] = mat[1][0];
565        m[0][2] = mat[2][0];
566        m[0][3] = mat[3][0];
567
568        m[1][0] = mat[0][1];
569        m[1][1] = mat[1][1];
570        m[1][2] = mat[2][1];
571        m[1][3] = mat[3][1];
572
573        m[2][0] = mat[0][2];
574        m[2][1] = mat[1][2];
575        m[2][2] = mat[2][2];
576        m[2][3] = mat[3][2];
577
578        m[3][0] = mat[0][3];
579        m[3][1] = mat[1][3];
580        m[3][2] = mat[2][3];
581        m[3][3] = mat[3][3];
582
583        memcpy(mat, m, 16*sizeof(float));
584}
585
586//
587
588float CMatrix::determinant() const
589{
590        float ans;
591        float a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;
592
593        /* assign to individual variable names to aid selecting */
594        /* correct elements */
595
596        a1 = mat[0][0]; b1 = mat[0][1]; 
597        c1 = mat[0][2]; d1 = mat[0][3];
598
599        a2 = mat[1][0]; b2 = mat[1][1]; 
600        c2 = mat[1][2]; d2 = mat[1][3];
601
602        a3 = mat[2][0]; b3 = mat[2][1]; 
603        c3 = mat[2][2]; d3 = mat[2][3];
604
605        a4 = mat[3][0]; b4 = mat[3][1]; 
606        c4 = mat[3][2]; d4 = mat[3][3];
607
608        ans = a1 * det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4)
609                - b1 * det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4)
610                + c1 * det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4)
611                - d1 * det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4);
612
613        return ans;
614}
615
616//
617
618void CMatrix::adjoint(const CMatrix& in, CMatrix& out)
619{
620    float a1, a2, a3, a4, b1, b2, b3, b4;
621    float c1, c2, c3, c4, d1, d2, d3, d4;
622
623    /* assign to individual variable names to aid  */
624    /* selecting correct values  */
625
626        a1 = in.mat[0][0]; b1 = in.mat[0][1]; 
627        c1 = in.mat[0][2]; d1 = in.mat[0][3];
628
629        a2 = in.mat[1][0]; b2 = in.mat[1][1]; 
630        c2 = in.mat[1][2]; d2 = in.mat[1][3];
631
632        a3 = in.mat[2][0]; b3 = in.mat[2][1];
633        c3 = in.mat[2][2]; d3 = in.mat[2][3];
634
635        a4 = in.mat[3][0]; b4 = in.mat[3][1]; 
636        c4 = in.mat[3][2]; d4 = in.mat[3][3];
637
638
639    /* row column labeling reversed since we transpose rows & columns */
640
641    out.mat[0][0] =   det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4);
642    out.mat[1][0] = - det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4);
643    out.mat[2][0] =   det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4);
644    out.mat[3][0] = - det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4);
645       
646    out.mat[0][1] = - det3x3( b1, b3, b4, c1, c3, c4, d1, d3, d4);
647    out.mat[1][1] =   det3x3( a1, a3, a4, c1, c3, c4, d1, d3, d4);
648    out.mat[2][1] = - det3x3( a1, a3, a4, b1, b3, b4, d1, d3, d4);
649    out.mat[3][1] =   det3x3( a1, a3, a4, b1, b3, b4, c1, c3, c4);
650       
651    out.mat[0][2] =   det3x3( b1, b2, b4, c1, c2, c4, d1, d2, d4);
652    out.mat[1][2] = - det3x3( a1, a2, a4, c1, c2, c4, d1, d2, d4);
653    out.mat[2][2] =   det3x3( a1, a2, a4, b1, b2, b4, d1, d2, d4);
654    out.mat[3][2] = - det3x3( a1, a2, a4, b1, b2, b4, c1, c2, c4);
655       
656    out.mat[0][3] = - det3x3( b1, b2, b3, c1, c2, c3, d1, d2, d3);
657    out.mat[1][3] =   det3x3( a1, a2, a3, c1, c2, c3, d1, d2, d3);
658    out.mat[2][3] = - det3x3( a1, a2, a3, b1, b2, b3, d1, d2, d3);
659    out.mat[3][3] =   det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3);
660}
661
662//
663
Note: See TracBrowser for help on using the repository browser.