Merged ball-test branch.
[neverball] / share / vec3.c
1 /*
2  * Copyright (C) 2003 Robert Kooima
3  *
4  * NEVERBALL is  free software; you can redistribute  it and/or modify
5  * it under the  terms of the GNU General  Public License as published
6  * by the Free  Software Foundation; either version 2  of the License,
7  * or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful, but
10  * WITHOUT  ANY  WARRANTY;  without   even  the  implied  warranty  of
11  * MERCHANTABILITY or  FITNESS FOR A PARTICULAR PURPOSE.   See the GNU
12  * General Public License for more details.
13  */
14
15 #include <stdio.h>
16 #include <math.h>
17
18 #include "vec3.h"
19
20 #define A 10
21 #define B 11
22 #define C 12
23 #define D 13
24 #define E 14
25 #define F 15
26
27 #define TINY 1e-5
28
29 /*---------------------------------------------------------------------------*/
30
31 void v_nrm(float *n, const float *v)
32 {
33     float d = v_len(v);
34
35     n[0] = v[0] / d;
36     n[1] = v[1] / d;
37     n[2] = v[2] / d;
38 }
39
40 void v_crs(float *u, const float *v, const float *w)
41 {
42     u[0] = v[1] * w[2] - v[2] * w[1];
43     u[1] = v[2] * w[0] - v[0] * w[2];
44     u[2] = v[0] * w[1] - v[1] * w[0];
45 }
46
47 /*---------------------------------------------------------------------------*/
48
49 void m_cpy(float *M, const float *N)
50 {
51     M[0] = N[0]; M[1] = N[1]; M[2] = N[2]; M[3] = N[3];
52     M[4] = N[4]; M[5] = N[5]; M[6] = N[6]; M[7] = N[7];
53     M[8] = N[8]; M[9] = N[9]; M[A] = N[A]; M[B] = N[B];
54     M[C] = N[C]; M[D] = N[D]; M[E] = N[E]; M[F] = N[F];
55 }
56
57 void m_xps(float *M, const float *N)
58 {
59     M[0] = N[0]; M[1] = N[4]; M[2] = N[8]; M[3] = N[C];
60     M[4] = N[1]; M[5] = N[5]; M[6] = N[9]; M[7] = N[D];
61     M[8] = N[2]; M[9] = N[6]; M[A] = N[A]; M[B] = N[E];
62     M[C] = N[3]; M[D] = N[7]; M[E] = N[B]; M[F] = N[F];
63 }
64
65 int  m_inv(float *I, const float *M)
66 {
67     float T[16];
68     float d;
69
70     T[0] = +(M[5] * (M[A] * M[F] - M[B] * M[E]) -
71              M[9] * (M[6] * M[F] - M[7] * M[E]) +
72              M[D] * (M[6] * M[B] - M[7] * M[A]));
73     T[1] = -(M[4] * (M[A] * M[F] - M[B] * M[E]) -
74              M[8] * (M[6] * M[F] - M[7] * M[E]) +
75              M[C] * (M[6] * M[B] - M[7] * M[A]));
76     T[2] = +(M[4] * (M[9] * M[F] - M[B] * M[D]) -
77              M[8] * (M[5] * M[F] - M[7] * M[D]) +
78              M[C] * (M[5] * M[B] - M[7] * M[9]));
79     T[3] = -(M[4] * (M[9] * M[E] - M[A] * M[D]) -
80              M[8] * (M[5] * M[E] - M[6] * M[D]) +
81              M[C] * (M[5] * M[A] - M[6] * M[9]));
82
83     T[4] = -(M[1] * (M[A] * M[F] - M[B] * M[E]) -
84              M[9] * (M[2] * M[F] - M[3] * M[E]) +
85              M[D] * (M[2] * M[B] - M[3] * M[A]));
86     T[5] = +(M[0] * (M[A] * M[F] - M[B] * M[E]) -
87              M[8] * (M[2] * M[F] - M[3] * M[E]) +
88              M[C] * (M[2] * M[B] - M[3] * M[A]));
89     T[6] = -(M[0] * (M[9] * M[F] - M[B] * M[D]) -
90              M[8] * (M[1] * M[F] - M[3] * M[D]) +
91              M[C] * (M[1] * M[B] - M[3] * M[9]));
92     T[7] = +(M[0] * (M[9] * M[E] - M[A] * M[D]) -
93              M[8] * (M[1] * M[E] - M[2] * M[D]) +
94              M[C] * (M[1] * M[A] - M[2] * M[9]));
95
96     T[8] = +(M[1] * (M[6] * M[F] - M[7] * M[E]) -
97              M[5] * (M[2] * M[F] - M[3] * M[E]) +
98              M[D] * (M[2] * M[7] - M[3] * M[6]));
99     T[9] = -(M[0] * (M[6] * M[F] - M[7] * M[E]) -
100              M[4] * (M[2] * M[F] - M[3] * M[E]) +
101              M[C] * (M[2] * M[7] - M[3] * M[6]));
102     T[A] = +(M[0] * (M[5] * M[F] - M[7] * M[D]) -
103              M[4] * (M[1] * M[F] - M[3] * M[D]) +
104              M[C] * (M[1] * M[7] - M[3] * M[5]));
105     T[B] = -(M[0] * (M[5] * M[E] - M[6] * M[D]) -
106              M[4] * (M[1] * M[E] - M[2] * M[D]) +
107              M[C] * (M[1] * M[6] - M[2] * M[5]));
108
109     T[C] = -(M[1] * (M[6] * M[B] - M[7] * M[A]) -
110              M[5] * (M[2] * M[B] - M[3] * M[A]) +
111              M[9] * (M[2] * M[7] - M[3] * M[6]));
112     T[D] = +(M[0] * (M[6] * M[B] - M[7] * M[A]) -
113              M[4] * (M[2] * M[B] - M[3] * M[A]) +
114              M[8] * (M[2] * M[7] - M[3] * M[6]));
115     T[E] = -(M[0] * (M[5] * M[B] - M[7] * M[9]) -
116              M[4] * (M[1] * M[B] - M[3] * M[9]) +
117              M[8] * (M[1] * M[7] - M[3] * M[5]));
118     T[F] = +(M[0] * (M[5] * M[A] - M[6] * M[9]) -
119              M[4] * (M[1] * M[A] - M[2] * M[9]) +
120              M[8] * (M[1] * M[6] - M[2] * M[5]));
121
122     d = M[0] * T[0] + M[4] * T[4] + M[8] * T[8] + M[C] * T[C];
123
124     if (fabs(d) > TINY)
125     {
126         I[0] = T[0] / d;
127         I[1] = T[4] / d;
128         I[2] = T[8] / d;
129         I[3] = T[C] / d;
130         I[4] = T[1] / d;
131         I[5] = T[5] / d;
132         I[6] = T[9] / d;
133         I[7] = T[D] / d;
134         I[8] = T[2] / d;
135         I[9] = T[6] / d;
136         I[A] = T[A] / d;
137         I[B] = T[E] / d;
138         I[C] = T[3] / d;
139         I[D] = T[7] / d;
140         I[E] = T[B] / d;
141         I[F] = T[F] / d;
142
143         return 1;
144     }
145     return 0;
146 }
147
148 /*---------------------------------------------------------------------------*/
149
150 void m_ident(float *M)
151 {
152     M[0] = 1.f; M[4] = 0.f; M[8] = 0.f; M[C] = 0.f;
153     M[1] = 0.f; M[5] = 1.f; M[9] = 0.f; M[D] = 0.f;
154     M[2] = 0.f; M[6] = 0.f; M[A] = 1.f; M[E] = 0.f;
155     M[3] = 0.f; M[7] = 0.f; M[B] = 0.f; M[F] = 1.f;
156 }
157
158 void m_basis(float *M,
159              const float e0[3],
160              const float e1[3],
161              const float e2[3])
162 {
163     M[0] = e0[0]; M[4] = e1[0]; M[8] = e2[0]; M[C] = 0.f;
164     M[1] = e0[1]; M[5] = e1[1]; M[9] = e2[1]; M[D] = 0.f;
165     M[2] = e0[2]; M[6] = e1[2]; M[A] = e2[2]; M[E] = 0.f;
166     M[3] =   0.f; M[7] =   0.f; M[B] =   0.f; M[F] = 1.f;
167 }
168
169 /*---------------------------------------------------------------------------*/
170
171 void m_xlt(float *M, const float *v)
172 {
173     M[0] = 1.f; M[4] = 0.f; M[8] = 0.f; M[C] = v[0];
174     M[1] = 0.f; M[5] = 1.f; M[9] = 0.f; M[D] = v[1];
175     M[2] = 0.f; M[6] = 0.f; M[A] = 1.f; M[E] = v[2];
176     M[3] = 0.f; M[7] = 0.f; M[B] = 0.f; M[F] =  1.f;
177 }
178
179 void m_scl(float *M, const float *v)
180 {
181     M[0] = v[0]; M[4] =  0.f; M[8] =  0.f; M[C] = 0.f;
182     M[1] =  0.f; M[5] = v[1]; M[9] =  0.f; M[D] = 0.f;
183     M[2] =  0.f; M[6] =  0.f; M[A] = v[2]; M[E] = 0.f;
184     M[3] =  0.f; M[7] =  0.f; M[B] =  0.f; M[F] = 1.f;
185 }
186
187 void m_rot(float *M, const float *v, float a)
188 {
189     float u[3];
190     float U[16];
191     float S[16];
192
193     const float s = fsinf(a);
194     const float c = fcosf(a);
195
196     v_nrm(u, v);
197
198     U[0] = u[0] * u[0]; U[4] = u[0] * u[1]; U[8] = u[0] * u[2];
199     U[1] = u[1] * u[0]; U[5] = u[1] * u[1]; U[9] = u[1] * u[2];
200     U[2] = u[2] * u[0]; U[6] = u[2] * u[1]; U[A] = u[2] * u[2];
201
202     S[0] =   0.f; S[4] = -u[2]; S[8] =  u[1];
203     S[1] =  u[2]; S[5] =   0.f; S[9] = -u[0];
204     S[2] = -u[1]; S[6] =  u[0]; S[A] =   0.f;
205
206     M[0] = U[0] + c * (1.f - U[0]) + s * S[0];
207     M[1] = U[1] + c * (0.f - U[1]) + s * S[1];
208     M[2] = U[2] + c * (0.f - U[2]) + s * S[2];
209     M[3] = 0.f;
210     M[4] = U[4] + c * (0.f - U[4]) + s * S[4];
211     M[5] = U[5] + c * (1.f - U[5]) + s * S[5];
212     M[6] = U[6] + c * (0.f - U[6]) + s * S[6];
213     M[7] = 0.f;
214     M[8] = U[8] + c * (0.f - U[8]) + s * S[8];
215     M[9] = U[9] + c * (0.f - U[9]) + s * S[9];
216     M[A] = U[A] + c * (1.f - U[A]) + s * S[A];
217     M[B] = 0.f;
218     M[C] = 0.f;
219     M[D] = 0.f;
220     M[E] = 0.f;
221     M[F] = 1.f;
222 }
223
224 /*---------------------------------------------------------------------------*/
225
226 void m_mult(float *M, const float *N, const float *O)
227 {
228     M[0] = N[0] * O[0] + N[4] * O[1] + N[8] * O[2] + N[C] * O[3];
229     M[1] = N[1] * O[0] + N[5] * O[1] + N[9] * O[2] + N[D] * O[3];
230     M[2] = N[2] * O[0] + N[6] * O[1] + N[A] * O[2] + N[E] * O[3];
231     M[3] = N[3] * O[0] + N[7] * O[1] + N[B] * O[2] + N[F] * O[3];
232
233     M[4] = N[0] * O[4] + N[4] * O[5] + N[8] * O[6] + N[C] * O[7];
234     M[5] = N[1] * O[4] + N[5] * O[5] + N[9] * O[6] + N[D] * O[7];
235     M[6] = N[2] * O[4] + N[6] * O[5] + N[A] * O[6] + N[E] * O[7];
236     M[7] = N[3] * O[4] + N[7] * O[5] + N[B] * O[6] + N[F] * O[7];
237
238     M[8] = N[0] * O[8] + N[4] * O[9] + N[8] * O[A] + N[C] * O[B];
239     M[9] = N[1] * O[8] + N[5] * O[9] + N[9] * O[A] + N[D] * O[B];
240     M[A] = N[2] * O[8] + N[6] * O[9] + N[A] * O[A] + N[E] * O[B];
241     M[B] = N[3] * O[8] + N[7] * O[9] + N[B] * O[A] + N[F] * O[B];
242
243     M[C] = N[0] * O[C] + N[4] * O[D] + N[8] * O[E] + N[C] * O[F];
244     M[D] = N[1] * O[C] + N[5] * O[D] + N[9] * O[E] + N[D] * O[F];
245     M[E] = N[2] * O[C] + N[6] * O[D] + N[A] * O[E] + N[E] * O[F];
246     M[F] = N[3] * O[C] + N[7] * O[D] + N[B] * O[E] + N[F] * O[F];
247 }
248
249 void m_pxfm(float *v, const float *M, const float *w)
250 {
251     float d = w[0] * M[3] + w[1] * M[7] + w[2] * M[B] + M[F];
252
253     v[0] = (w[0] * M[0] + w[1] * M[4] + w[2] * M[8] + M[C]) / d;
254     v[1] = (w[0] * M[1] + w[1] * M[5] + w[2] * M[9] + M[D]) / d;
255     v[2] = (w[0] * M[2] + w[1] * M[6] + w[2] * M[A] + M[E]) / d;
256 }
257
258 void m_vxfm(float *v, const float *M, const float *w)
259 {
260     v[0] = (w[0] * M[0] + w[1] * M[4] + w[2] * M[8]);
261     v[1] = (w[0] * M[1] + w[1] * M[5] + w[2] * M[9]);
262     v[2] = (w[0] * M[2] + w[1] * M[6] + w[2] * M[A]);
263 }
264
265 /*---------------------------------------------------------------------------*/
266
267 void m_view(float *M,
268             const float *c,
269             const float *p,
270             const float *u)
271 {
272     float x[3];
273     float y[3];
274     float z[3];
275
276     v_sub(z, p, c);
277     v_nrm(z, z);
278     v_crs(x, u, z);
279     v_nrm(x, x);
280     v_crs(y, z, x);
281
282     m_basis(M, x, y, z);
283 }
284
285 /*---------------------------------------------------------------------------*/