source: trunk/packages/vizservers/nanovis/Mat4x4.cpp @ 1028

Last change on this file since 1028 was 1028, checked in by gah, 13 years ago

various cleanups

File size: 5.6 KB
Line 
1
2/*
3 * ----------------------------------------------------------------------
4 * Mat4x4.cpp: Mat4x4 class
5 *
6 * ======================================================================
7 *  AUTHOR:  Wei Qiao <qiaow@purdue.edu>
8 *           Purdue Rendering and Perceptualization Lab (PURPL)
9 *
10 *  Copyright (c) 2004-2006  Purdue Research Foundation
11 *
12 *  See the file "license.terms" for information on usage and
13 *  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
14 * ======================================================================
15 */
16#include "Mat4x4.h"
17#include <math.h>
18#include <stdio.h>
19#include <assert.h>
20
21Mat4x4::Mat4x4(float *vals) {
22    int i;
23    for (i=0;i<16;i++)
24        m[i] = vals[i];
25}
26
27void
28Mat4x4::print(){
29    fprintf(stderr, "\n%f %f %f %f\n", m[0], m[1], m[2], m[3]);
30    fprintf(stderr, "%f %f %f %f\n", m[4], m[5], m[6], m[7]);
31    fprintf(stderr, "%f %f %f %f\n", m[8], m[9], m[10], m[11]);
32    fprintf(stderr, "%f %f %f %f\n\n", m[12], m[13], m[14], m[15]);
33}
34
35Mat4x4
36Mat4x4::inverse(){
37    Mat4x4 p;
38       
39    float m00 = m[0], m01 = m[4], m02 = m[8], m03 = m[12];
40    float m10 = m[1], m11 = m[5], m12 = m[9], m13 = m[13];
41    float m20 = m[2], m21 = m[6], m22 = m[10], m23 = m[14];
42    float m30 = m[3], m31 = m[7], m32 = m[11], m33 = m[15];
43
44    float det0, inv0;
45
46#define det33(a1,a2,a3,b1,b2,b3,c1,c2,c3) \
47    (a1*(b2*c3-b3*c2) + b1*(c2*a3-a2*c3) + c1*(a2*b3-a3*b2))
48
49    float cof00 =  det33(m11, m12, m13, m21, m22, m23, m31, m32, m33);
50    float cof01 = -det33(m12, m13, m10, m22, m23, m20, m32, m33, m30);
51    float cof02 =  det33(m13, m10, m11, m23, m20, m21, m33, m30, m31);
52    float cof03 = -det33(m10, m11, m12, m20, m21, m22, m30, m31, m32);
53   
54    float cof10 = -det33(m21, m22, m23, m31, m32, m33, m01, m02, m03);
55    float cof11 =  det33(m22, m23, m20, m32, m33, m30, m02, m03, m00);
56    float cof12 = -det33(m23, m20, m21, m33, m30, m31, m03, m00, m01);
57    float cof13 =  det33(m20, m21, m22, m30, m31, m32, m00, m01, m02);
58   
59    float cof20 =  det33(m31, m32, m33, m01, m02, m03, m11, m12, m13);
60    float cof21 = -det33(m32, m33, m30, m02, m03, m00, m12, m13, m10);
61    float cof22 =  det33(m33, m30, m31, m03, m00, m01, m13, m10, m11);
62    float cof23 = -det33(m30, m31, m32, m00, m01, m02, m10, m11, m12);
63   
64    float cof30 = -det33(m01, m02, m03, m11, m12, m13, m21, m22, m23);
65    float cof31 =  det33(m02, m03, m00, m12, m13, m10, m22, m23, m20);
66    float cof32 = -det33(m03, m00, m01, m13, m10, m11, m23, m20, m21);
67    float cof33 =  det33(m00, m01, m02, m10, m11, m12, m20, m21, m22);
68
69#undef det33
70
71#if 0
72    fprintf(stderr, "Invert:\n");
73    fprintf(stderr, "   %12.9f %12.9f %12.9f %12.9f\n", m00, m01, m02, m03);
74    fprintf(stderr, "   %12.9f %12.9f %12.9f %12.9f\n", m10, m11, m12, m13);
75    fprintf(stderr, "   %12.9f %12.9f %12.9f %12.9f\n", m20, m21, m22, m23);
76    fprintf(stderr, "   %12.9f %12.9f %12.9f %12.9f\n", m30, m31, m32, m33);
77#endif
78
79#if 0
80    det0  = m00 * cof00 + m01 * cof01 + m02 * cof02 + m03 * cof03;
81#else
82    det0  = m30 * cof30 + m31 * cof31 + m32 * cof32 + m33 * cof33;
83#endif
84    inv0  = (det0 == 0.0f) ? 0.0f : 1.0f / det0;
85
86    p.m[0] = cof00 * inv0;
87    p.m[1] = cof01 * inv0;
88    p.m[2] = cof02 * inv0;
89    p.m[3] = cof03 * inv0;
90
91    p.m[4] = cof10 * inv0;
92    p.m[5] = cof11 * inv0;
93    p.m[6] = cof12 * inv0;
94    p.m[7] = cof13 * inv0;
95
96    p.m[8] = cof20 * inv0;
97    p.m[9] = cof21 * inv0;
98    p.m[10] = cof22 * inv0;
99    p.m[11] = cof23 * inv0;
100
101    p.m[12] = cof30 * inv0;
102    p.m[13] = cof31 * inv0;
103    p.m[14] = cof32 * inv0;
104    p.m[15]= cof33 * inv0;     
105       
106    return p;
107}
108
109Vector4
110Mat4x4::multiply_row_vector(Vector4 v){
111    Vector4 ret;
112    ret.x = (m[0]  * v.x) + (m[1]  * v.y) + (m[2]  * v.z) + (m[3]  * v.w);
113    ret.y = (m[4]  * v.x) + (m[5]  * v.y) + (m[6]  * v.z) + (m[7]  * v.w);
114    ret.z = (m[8]  * v.x) + (m[9]  * v.y) + (m[10] * v.z) + (m[11] * v.w);
115    ret.w = (m[12] * v.x) + (m[13] * v.y) + (m[14] * v.z) + (m[15] * v.w);
116    return ret;
117}
118
119Vector4
120Mat4x4::transform(Vector4 v){
121    Vector4 ret;
122    ret.x = v.x*m[0]+v.y*m[4]+v.z*m[8]+v.w*m[12];
123    ret.y = v.x*m[1]+v.y*m[5]+v.z*m[9]+v.w*m[13];
124    ret.z = v.x*m[2]+v.y*m[6]+v.z*m[10]+v.w*m[14];
125    ret.w = v.x*m[3]+v.y*m[7]+v.z*m[11]+v.w*m[15];
126    return ret;
127}
128
129Mat4x4
130Mat4x4::operator *(Mat4x4 mat){
131    Mat4x4 ret;
132   
133    ret.m[0] = m[0]*mat.m[0]+m[1]*mat.m[4]+m[2]*mat.m[8]+m[3]*mat.m[12];
134    ret.m[1] = m[0]*mat.m[1]+m[1]*mat.m[5]+m[2]*mat.m[9]+m[3]*mat.m[13];
135    ret.m[2] = m[0]*mat.m[2]+m[1]*mat.m[6]+m[2]*mat.m[10]+m[3]*mat.m[14];
136    ret.m[3] = m[0]*mat.m[3]+m[1]*mat.m[7]+m[2]*mat.m[11]+m[3]*mat.m[15];
137   
138    ret.m[4] = m[4]*mat.m[0]+m[5]*mat.m[4]+m[6]*mat.m[8]+m[7]*mat.m[12];
139    ret.m[5] = m[4]*mat.m[1]+m[5]*mat.m[5]+m[6]*mat.m[9]+m[7]*mat.m[13];
140    ret.m[6] = m[4]*mat.m[2]+m[5]*mat.m[6]+m[6]*mat.m[10]+m[7]*mat.m[14];
141    ret.m[7] = m[4]*mat.m[3]+m[5]*mat.m[7]+m[6]*mat.m[11]+m[7]*mat.m[15];
142   
143    ret.m[8] = m[8]*mat.m[0]+m[9]*mat.m[4]+m[10]*mat.m[8]+m[11]*mat.m[12];
144    ret.m[9] = m[8]*mat.m[1]+m[9]*mat.m[5]+m[10]*mat.m[9]+m[11]*mat.m[13];
145    ret.m[10] = m[8]*mat.m[2]+m[9]*mat.m[6]+m[10]*mat.m[10]+m[11]*mat.m[14];
146    ret.m[11] = m[8]*mat.m[3]+m[9]*mat.m[7]+m[10]*mat.m[11]+m[11]*mat.m[15];
147   
148    ret.m[12] = m[12]*mat.m[0]+m[13]*mat.m[4]+m[14]*mat.m[8]+m[15]*mat.m[12];
149    ret.m[13] = m[12]*mat.m[1]+m[13]*mat.m[5]+m[14]*mat.m[9]+m[15]*mat.m[13];
150    ret.m[14] = m[12]*mat.m[2]+m[13]*mat.m[6]+m[14]*mat.m[10]+m[15]*mat.m[14];
151    ret.m[15] = m[12]*mat.m[3]+m[13]*mat.m[7]+m[14]*mat.m[11]+m[15]*mat.m[15];
152   
153    return ret;
154}
155
156Mat4x4
157Mat4x4::transpose(){
158    Mat4x4 ret;
159    for(int i=0; i<4; i++){
160        for(int j=0; j<4; j++){
161            ret.m[j+4*i] = this->m[i+4*j];
162        }
163    }
164    return ret;
165}
166
167
Note: See TracBrowser for help on using the repository browser.