« back to all changes in this revision
Viewing changes to algebra3.cpp
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Marcelo E. Magallon
- Date: 2001-09-23 12:57:28 UTC
- Revision ID: james.westby@ubuntu.com-20010923125728-qbts7xit2eg1ogo2
Tags: upstream-2.1.0.beta
Import upstream version 2.1.0.beta
- files added:
- msvc
added
removed
algebra3.cpp
1
/**************************************************************************
2
3
algebra3.cpp, algebra3.h - C++ Vector and Matrix Algebra routines
4
5
There are three vector classes and two matrix classes: vec2, vec3,
6
vec4, mat3, and mat4.
7
8
All the standard arithmetic operations are defined, with '*'
9
for dot product of two vectors and multiplication of two matrices,
10
and '^' for cross product of two vectors.
11
12
Additional functions include length(), normalize(), homogenize for
13
vectors, and print(), set(), apply() for all classes.
14
15
There is a function transpose() for matrices, but note that it
16
does not actually change the matrix,
17
18
When multiplied with a matrix, a vector is treated as a row vector
19
if it precedes the matrix (v*M), and as a column vector if it
20
follows the matrix (M*v).
21
22
Matrices are stored in row-major form.
23
24
A vector of one dimension (2d, 3d, or 4d) can be cast to a vector
25
of a higher or lower dimension. If casting to a higher dimension,
26
the new component is set by default to 1.0, unless a value is
27
specified:
28
vec3 a(1.0, 2.0, 3.0 );
29
vec4 b( a, 4.0 ); // now b == {1.0, 2.0, 3.0, 4.0};
30
When casting to a lower dimension, the vector is homogenized in
31
the lower dimension. E.g., if a 4d {X,Y,Z,W} is cast to 3d, the
32
resulting vector is {X/W, Y/W, Z/W}. It is up to the user to
33
insure the fourth component is not zero before casting.
34
35
There are also the following function for building matrices:
36
identity2D(), translation2D(), rotation2D(),
37
scaling2D(), identity3D(), translation3D(),
38
rotation3D(), rotation3Drad(), scaling3D(),
39
perspective3D()
40
41
42
---------------------------------------------------------------------
43
44
Author: Jean-Francois DOUEg
45
Revised: Paul Rademacher
46
Version 3.2 - Feb 1998
47
48
**************************************************************************/
49
50
#include <math.h>
51
#include "algebra3.h"
52
#include <ctype.h>
53
54
/****************************************************************
55
* *
56
* vec2 Member functions *
57
* *
58
****************************************************************/
59
60
/******************** vec2 CONSTRUCTORS ********************/
61
62
vec2::vec2(void)
63
{n[VX] = n[VY] = 0.0; }
64
65
vec2::vec2(const float x, const float y)
66
{ n[VX] = x; n[VY] = y; }
67
68
vec2::vec2(const float d)
69
{ n[VX] = n[VY] = d; }
70
71
vec2::vec2(const vec2& v)
72
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; }
73
74
vec2::vec2(const vec3& v) // it is up to caller to avoid divide-by-zero
75
{ n[VX] = v.n[VX]/v.n[VZ]; n[VY] = v.n[VY]/v.n[VZ]; };
76
77
vec2::vec2(const vec3& v, int dropAxis) {
78
switch (dropAxis) {
79
case VX: n[VX] = v.n[VY]; n[VY] = v.n[VZ]; break;
80
case VY: n[VX] = v.n[VX]; n[VY] = v.n[VZ]; break;
81
default: n[VX] = v.n[VX]; n[VY] = v.n[VY]; break;
82
}
83
}
84
85
86
/******************** vec2 ASSIGNMENT OPERATORS ******************/
87
88
vec2& vec2::operator = (const vec2& v)
89
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; return *this; }
90
91
vec2& vec2::operator += ( const vec2& v )
92
{ n[VX] += v.n[VX]; n[VY] += v.n[VY]; return *this; }
93
94
vec2& vec2::operator -= ( const vec2& v )
95
{ n[VX] -= v.n[VX]; n[VY] -= v.n[VY]; return *this; }
96
97
vec2& vec2::operator *= ( const float d )
98
{ n[VX] *= d; n[VY] *= d; return *this; }
99
100
vec2& vec2::operator /= ( const float d )
101
{ float d_inv = 1./d; n[VX] *= d_inv; n[VY] *= d_inv; return *this; }
102
103
float& vec2::operator [] ( int i) {
104
if (i < VX || i > VY)
105
//VEC_ERROR("vec2 [] operator: illegal access; index = " << i << '\n')
106
VEC_ERROR("vec2 [] operator: illegal access" );
107
return n[i];
108
}
109
110
111
/******************** vec2 SPECIAL FUNCTIONS ********************/
112
113
float vec2::length(void)
114
{ return sqrt(length2()); }
115
116
float vec2::length2(void)
117
{ return n[VX]*n[VX] + n[VY]*n[VY]; }
118
119
vec2& vec2::normalize(void) // it is up to caller to avoid divide-by-zero
120
{ *this /= length(); return *this; }
121
122
vec2& vec2::apply(V_FCT_PTR fct)
123
{ n[VX] = (*fct)(n[VX]); n[VY] = (*fct)(n[VY]); return *this; }
124
125
void vec2::set( float x, float y )
126
{ n[VX] = x; n[VY] = y; }
127
128
/******************** vec2 FRIENDS *****************************/
129
130
vec2 operator - (const vec2& a)
131
{ return vec2(-a.n[VX],-a.n[VY]); }
132
133
vec2 operator + (const vec2& a, const vec2& b)
134
{ return vec2(a.n[VX]+ b.n[VX], a.n[VY] + b.n[VY]); }
135
136
vec2 operator - (const vec2& a, const vec2& b)
137
{ return vec2(a.n[VX]-b.n[VX], a.n[VY]-b.n[VY]); }
138
139
vec2 operator * (const vec2& a, const float d)
140
{ return vec2(d*a.n[VX], d*a.n[VY]); }
141
142
vec2 operator * (const float d, const vec2& a)
143
{ return a*d; }
144
145
vec2 operator * (const mat3& a, const vec2& v) {
146
vec3 av;
147
148
av.n[VX] = a.v[0].n[VX]*v.n[VX] + a.v[0].n[VY]*v.n[VY] + a.v[0].n[VZ];
149
av.n[VY] = a.v[1].n[VX]*v.n[VX] + a.v[1].n[VY]*v.n[VY] + a.v[1].n[VZ];
150
av.n[VZ] = a.v[2].n[VX]*v.n[VX] + a.v[2].n[VY]*v.n[VY] + a.v[2].n[VZ];
151
return av;
152
}
153
154
vec2 operator * (const vec2& v, mat3& a)
155
{ return a.transpose() * v; }
156
157
vec3 operator * (const mat3& a, const vec3& v) {
158
vec3 av;
159
160
av.n[VX] =
161
a.v[0].n[VX]*v.n[VX] + a.v[0].n[VY]*v.n[VY] + a.v[0].n[VZ]*v.n[VZ];
162
av.n[VY] =
163
a.v[1].n[VX]*v.n[VX] + a.v[1].n[VY]*v.n[VY] + a.v[1].n[VZ]*v.n[VZ];
164
av.n[VZ] =
165
a.v[2].n[VX]*v.n[VX] + a.v[2].n[VY]*v.n[VY] + a.v[2].n[VZ]*v.n[VZ];
166
return av;
167
}
168
169
vec3 operator * (const vec3& v, mat3& a)
170
{ return a.transpose() * v; }
171
172
float operator * (const vec2& a, const vec2& b)
173
{ return (a.n[VX]*b.n[VX] + a.n[VY]*b.n[VY]); }
174
175
vec2 operator / (const vec2& a, const float d)
176
{ float d_inv = 1./d; return vec2(a.n[VX]*d_inv, a.n[VY]*d_inv); }
177
178
vec3 operator ^ (const vec2& a, const vec2& b)
179
{ return vec3(0.0, 0.0, a.n[VX] * b.n[VY] - b.n[VX] * a.n[VY]); }
180
181
int operator == (const vec2& a, const vec2& b)
182
{ return (a.n[VX] == b.n[VX]) && (a.n[VY] == b.n[VY]); }
183
184
int operator != (const vec2& a, const vec2& b)
185
{ return !(a == b); }
186
187
/*ostream& operator << (ostream& s, vec2& v)
188
{ return s << "| " << v.n[VX] << ' ' << v.n[VY] << " |"; }
189
*/
190
191
/*istream& operator >> (istream& s, vec2& v) {
192
vec2 v_tmp;
193
char c = ' ';
194
195
while (isspace(c))
196
s >> c;
197
// The vectors can be formatted either as x y or | x y |
198
if (c == '|') {
199
s >> v_tmp[VX] >> v_tmp[VY];
200
while (s >> c && isspace(c)) ;
201
if (c != '|')
202
;//s.set(_bad);
203
}
204
else {
205
s.putback(c);
206
s >> v_tmp[VX] >> v_tmp[VY];
207
}
208
if (s)
209
v = v_tmp;
210
return s;
211
}
212
*/
213
214
void swap(vec2& a, vec2& b)
215
{ vec2 tmp(a); a = b; b = tmp; }
216
217
vec2 min_vec(const vec2& a, const vec2& b)
218
{ return vec2(MIN(a.n[VX], b.n[VX]), MIN(a.n[VY], b.n[VY])); }
219
220
vec2 max_vec(const vec2& a, const vec2& b)
221
{ return vec2(MAX(a.n[VX], b.n[VX]), MAX(a.n[VY], b.n[VY])); }
222
223
vec2 prod(const vec2& a, const vec2& b)
224
{ return vec2(a.n[VX] * b.n[VX], a.n[VY] * b.n[VY]); }
225
226
/****************************************************************
227
* *
228
* vec3 Member functions *
229
* *
230
****************************************************************/
231
232
// CONSTRUCTORS
233
234
vec3::vec3(void)
235
{n[VX] = n[VY] = n[VZ] = 0.0;}
236
237
vec3::vec3(const float x, const float y, const float z)
238
{ n[VX] = x; n[VY] = y; n[VZ] = z; }
239
240
vec3::vec3(const float d)
241
{ n[VX] = n[VY] = n[VZ] = d; }
242
243
vec3::vec3(const vec3& v)
244
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; }
245
246
vec3::vec3(const vec2& v)
247
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = 1.0; }
248
249
vec3::vec3(const vec2& v, float d)
250
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = d; }
251
252
vec3::vec3(const vec4& v) // it is up to caller to avoid divide-by-zero
253
{ n[VX] = v.n[VX] / v.n[VW]; n[VY] = v.n[VY] / v.n[VW];
254
n[VZ] = v.n[VZ] / v.n[VW]; }
255
256
vec3::vec3(const vec4& v, int dropAxis) {
257
switch (dropAxis) {
258
case VX: n[VX] = v.n[VY]; n[VY] = v.n[VZ]; n[VZ] = v.n[VW]; break;
259
case VY: n[VX] = v.n[VX]; n[VY] = v.n[VZ]; n[VZ] = v.n[VW]; break;
260
case VZ: n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VW]; break;
261
default: n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; break;
262
}
263
}
264
265
266
// ASSIGNMENT OPERATORS
267
268
vec3& vec3::operator = (const vec3& v)
269
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; return *this; }
270
271
vec3& vec3::operator += ( const vec3& v )
272
{ n[VX] += v.n[VX]; n[VY] += v.n[VY]; n[VZ] += v.n[VZ]; return *this; }
273
274
vec3& vec3::operator -= ( const vec3& v )
275
{ n[VX] -= v.n[VX]; n[VY] -= v.n[VY]; n[VZ] -= v.n[VZ]; return *this; }
276
277
vec3& vec3::operator *= ( const float d )
278
{ n[VX] *= d; n[VY] *= d; n[VZ] *= d; return *this; }
279
280
vec3& vec3::operator /= ( const float d )
281
{ float d_inv = 1./d; n[VX] *= d_inv; n[VY] *= d_inv; n[VZ] *= d_inv;
282
return *this; }
283
284
float& vec3::operator [] ( int i) {
285
if (i < VX || i > VZ)
286
//VEC_ERROR("vec3 [] operator: illegal access; index = " << i << '\n')
287
VEC_ERROR("vec3 [] operator: illegal access" );
288
289
return n[i];
290
}
291
292
293
// SPECIAL FUNCTIONS
294
295
float vec3::length(void)
296
{ return sqrt(length2()); }
297
298
float vec3::length2(void)
299
{ return n[VX]*n[VX] + n[VY]*n[VY] + n[VZ]*n[VZ]; }
300
301
vec3& vec3::normalize(void) // it is up to caller to avoid divide-by-zero
302
{ *this /= length(); return *this; }
303
304
vec3& vec3::homogenize(void) // it is up to caller to avoid divide-by-zero
305
{ n[VX] /= n[VZ]; n[VY] /= n[VZ]; n[VZ] = 1.0; return *this; }
306
307
vec3& vec3::apply(V_FCT_PTR fct)
308
{ n[VX] = (*fct)(n[VX]); n[VY] = (*fct)(n[VY]); n[VZ] = (*fct)(n[VZ]);
309
return *this; }
310
311
void vec3::set( float x, float y, float z ) // set vector
312
{ n[VX] = x; n[VY] = y; n[VZ] = z; }
313
314
void vec3::print( FILE *file, char *name ) // print vector to a file
315
{
316
fprintf( file, "%s: <%f, %f, %f>\n", name, n[VX], n[VY], n[VZ] );
317
}
318
319
// FRIENDS
320
321
vec3 operator - (const vec3& a)
322
{ return vec3(-a.n[VX],-a.n[VY],-a.n[VZ]); }
323
324
vec3 operator + (const vec3& a, const vec3& b)
325
{ return vec3(a.n[VX]+ b.n[VX], a.n[VY] + b.n[VY], a.n[VZ] + b.n[VZ]); }
326
327
vec3 operator - (const vec3& a, const vec3& b)
328
{ return vec3(a.n[VX]-b.n[VX], a.n[VY]-b.n[VY], a.n[VZ]-b.n[VZ]); }
329
330
vec3 operator * (const vec3& a, const float d)
331
{ return vec3(d*a.n[VX], d*a.n[VY], d*a.n[VZ]); }
332
333
vec3 operator * (const float d, const vec3& a)
334
{ return a*d; }
335
336
vec3 operator * (const mat4& a, const vec3& v)
337
{ return a * vec4(v); }
338
339
vec3 operator * (const vec3& v, mat4& a)
340
{ return a.transpose() * v; }
341
342
float operator * (const vec3& a, const vec3& b)
343
{ return (a.n[VX]*b.n[VX] + a.n[VY]*b.n[VY] + a.n[VZ]*b.n[VZ]); }
344
345
vec3 operator / (const vec3& a, const float d)
346
{ float d_inv = 1./d; return vec3(a.n[VX]*d_inv, a.n[VY]*d_inv,
347
a.n[VZ]*d_inv); }
348
349
vec3 operator ^ (const vec3& a, const vec3& b) {
350
return vec3(a.n[VY]*b.n[VZ] - a.n[VZ]*b.n[VY],
351
a.n[VZ]*b.n[VX] - a.n[VX]*b.n[VZ],
352
a.n[VX]*b.n[VY] - a.n[VY]*b.n[VX]);
353
}
354
355
int operator == (const vec3& a, const vec3& b)
356
{ return (a.n[VX] == b.n[VX]) && (a.n[VY] == b.n[VY]) && (a.n[VZ] == b.n[VZ]);
357
}
358
359
int operator != (const vec3& a, const vec3& b)
360
{ return !(a == b); }
361
362
/*ostream& operator << (ostream& s, vec3& v)
363
{ return s << "| " << v.n[VX] << ' ' << v.n[VY] << ' ' << v.n[VZ] << " |"; }
364
365
istream& operator >> (istream& s, vec3& v) {
366
vec3 v_tmp;
367
char c = ' ';
368
369
while (isspace(c))
370
s >> c;
371
// The vectors can be formatted either as x y z or | x y z |
372
if (c == '|') {
373
s >> v_tmp[VX] >> v_tmp[VY] >> v_tmp[VZ];
374
while (s >> c && isspace(c)) ;
375
if (c != '|')
376
;//s.set(_bad);
377
}
378
else {
379
s.putback(c);
380
s >> v_tmp[VX] >> v_tmp[VY] >> v_tmp[VZ];
381
}
382
if (s)
383
v = v_tmp;
384
return s;
385
}
386
*/
387
388
void swap(vec3& a, vec3& b)
389
{ vec3 tmp(a); a = b; b = tmp; }
390
391
vec3 min_vec(const vec3& a, const vec3& b)
392
{ return vec3(MIN(a.n[VX], b.n[VX]), MIN(a.n[VY], b.n[VY]), MIN(a.n[VZ],
393
b.n[VZ])); }
394
395
vec3 max_vec(const vec3& a, const vec3& b)
396
{ return vec3(MAX(a.n[VX], b.n[VX]), MAX(a.n[VY], b.n[VY]), MAX(a.n[VZ],
397
b.n[VZ])); }
398
399
vec3 prod(const vec3& a, const vec3& b)
400
{ return vec3(a.n[VX] * b.n[VX], a.n[VY] * b.n[VY], a.n[VZ] * b.n[VZ]); }
401
402
/****************************************************************
403
* *
404
* vec4 Member functions *
405
* *
406
****************************************************************/
407
408
// CONSTRUCTORS
409
410
vec4::vec4(void)
411
{n[VX] = n[VY] = n[VZ] = 0.0; n[VW] = 1.0; }
412
413
vec4::vec4(const float x, const float y, const float z, const float w)
414
{ n[VX] = x; n[VY] = y; n[VZ] = z; n[VW] = w; }
415
416
vec4::vec4(const float d)
417
{ n[VX] = n[VY] = n[VZ] = n[VW] = d; }
418
419
vec4::vec4(const vec4& v)
420
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; n[VW] = v.n[VW]; }
421
422
vec4::vec4(const vec3& v)
423
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; n[VW] = 1.0; }
424
425
vec4::vec4(const vec3& v, const float d)
426
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; n[VW] = d; }
427
428
429
// ASSIGNMENT OPERATORS
430
431
vec4& vec4::operator = (const vec4& v)
432
{ n[VX] = v.n[VX]; n[VY] = v.n[VY]; n[VZ] = v.n[VZ]; n[VW] = v.n[VW];
433
return *this; }
434
435
vec4& vec4::operator += ( const vec4& v )
436
{ n[VX] += v.n[VX]; n[VY] += v.n[VY]; n[VZ] += v.n[VZ]; n[VW] += v.n[VW];
437
return *this; }
438
439
vec4& vec4::operator -= ( const vec4& v )
440
{ n[VX] -= v.n[VX]; n[VY] -= v.n[VY]; n[VZ] -= v.n[VZ]; n[VW] -= v.n[VW];
441
return *this; }
442
443
vec4& vec4::operator *= ( const float d )
444
{ n[VX] *= d; n[VY] *= d; n[VZ] *= d; n[VW] *= d; return *this; }
445
446
vec4& vec4::operator /= ( const float d )
447
{ float d_inv = 1./d; n[VX] *= d_inv; n[VY] *= d_inv; n[VZ] *= d_inv;
448
n[VW] *= d_inv; return *this; }
449
450
float& vec4::operator [] ( int i) {
451
if (i < VX || i > VW)
452
//VEC_ERROR("vec4 [] operator: illegal access; index = " << i << '\n')
453
VEC_ERROR("vec4 [] operator: illegal access" );
454
455
return n[i];
456
}
457
458
459
// SPECIAL FUNCTIONS
460
461
float vec4::length(void)
462
{ return sqrt(length2()); }
463
464
float vec4::length2(void)
465
{ return n[VX]*n[VX] + n[VY]*n[VY] + n[VZ]*n[VZ] + n[VW]*n[VW]; }
466
467
vec4& vec4::normalize(void) // it is up to caller to avoid divide-by-zero
468
{ *this /= length(); return *this; }
469
470
vec4& vec4::homogenize(void) // it is up to caller to avoid divide-by-zero
471
{ n[VX] /= n[VW]; n[VY] /= n[VW]; n[VZ] /= n[VW]; n[VW] = 1.0; return *this; }
472
473
vec4& vec4::apply(V_FCT_PTR fct)
474
{ n[VX] = (*fct)(n[VX]); n[VY] = (*fct)(n[VY]); n[VZ] = (*fct)(n[VZ]);
475
n[VW] = (*fct)(n[VW]); return *this; }
476
477
void vec4::print( FILE *file, char *name ) // print vector to a file
478
{
479
fprintf( file, "%s: <%f, %f, %f, %f>\n", name, n[VX], n[VY], n[VZ], n[VW] );
480
}
481
482
void vec4::set( float x, float y, float z, float a )
483
{
484
n[0] = x; n[1] = y; n[2] = z; n[3] = a;
485
}
486
487
488
// FRIENDS
489
490
vec4 operator - (const vec4& a)
491
{ return vec4(-a.n[VX],-a.n[VY],-a.n[VZ],-a.n[VW]); }
492
493
vec4 operator + (const vec4& a, const vec4& b)
494
{ return vec4(a.n[VX] + b.n[VX], a.n[VY] + b.n[VY], a.n[VZ] + b.n[VZ],
495
a.n[VW] + b.n[VW]); }
496
497
vec4 operator - (const vec4& a, const vec4& b)
498
{ return vec4(a.n[VX] - b.n[VX], a.n[VY] - b.n[VY], a.n[VZ] - b.n[VZ],
499
a.n[VW] - b.n[VW]); }
500
501
vec4 operator * (const vec4& a, const float d)
502
{ return vec4(d*a.n[VX], d*a.n[VY], d*a.n[VZ], d*a.n[VW] ); }
503
504
vec4 operator * (const float d, const vec4& a)
505
{ return a*d; }
506
507
vec4 operator * (const mat4& a, const vec4& v) {
508
#define ROWCOL(i) a.v[i].n[0]*v.n[VX] + a.v[i].n[1]*v.n[VY] \
509
+ a.v[i].n[2]*v.n[VZ] + a.v[i].n[3]*v.n[VW]
510
return vec4(ROWCOL(0), ROWCOL(1), ROWCOL(2), ROWCOL(3));
511
#undef ROWCOL
512
}
513
514
vec4 operator * (const vec4& v, mat4& a)
515
{ return a.transpose() * v; }
516
517
float operator * (const vec4& a, const vec4& b)
518
{ return (a.n[VX]*b.n[VX] + a.n[VY]*b.n[VY] + a.n[VZ]*b.n[VZ] +
519
a.n[VW]*b.n[VW]); }
520
521
vec4 operator / (const vec4& a, const float d)
522
{ float d_inv = 1./d; return vec4(a.n[VX]*d_inv, a.n[VY]*d_inv, a.n[VZ]*d_inv,
523
a.n[VW]*d_inv); }
524
525
int operator == (const vec4& a, const vec4& b)
526
{ return (a.n[VX] == b.n[VX]) && (a.n[VY] == b.n[VY]) && (a.n[VZ] == b.n[VZ])
527
&& (a.n[VW] == b.n[VW]); }
528
529
int operator != (const vec4& a, const vec4& b)
530
{ return !(a == b); }
531
532
/*ostream& operator << (ostream& s, vec4& v)
533
{ return s << "| " << v.n[VX] << ' ' << v.n[VY] << ' ' << v.n[VZ] << ' '
534
<< v.n[VW] << " |"; }
535
536
istream& operator >> (istream& s, vec4& v) {
537
vec4 v_tmp;
538
char c = ' ';
539
540
while (isspace(c))
541
s >> c;
542
// The vectors can be formatted either as x y z w or | x y z w |
543
if (c == '|') {
544
s >> v_tmp[VX] >> v_tmp[VY] >> v_tmp[VZ] >> v_tmp[VW];
545
while (s >> c && isspace(c)) ;
546
if (c != '|')
547
;//s.set(_bad);
548
}
549
else {
550
s.putback(c);
551
s >> v_tmp[VX] >> v_tmp[VY] >> v_tmp[VZ] >> v_tmp[VW];
552
}
553
if (s)
554
v = v_tmp;
555
return s;
556
}
557
*/
558
void swap(vec4& a, vec4& b)
559
{ vec4 tmp(a); a = b; b = tmp; }
560
561
vec4 min_vec(const vec4& a, const vec4& b)
562
{ return vec4(MIN(a.n[VX], b.n[VX]), MIN(a.n[VY], b.n[VY]), MIN(a.n[VZ],
563
b.n[VZ]), MIN(a.n[VW], b.n[VW])); }
564
565
vec4 max_vec(const vec4& a, const vec4& b)
566
{ return vec4(MAX(a.n[VX], b.n[VX]), MAX(a.n[VY], b.n[VY]), MAX(a.n[VZ],
567
b.n[VZ]), MAX(a.n[VW], b.n[VW])); }
568
569
vec4 prod(const vec4& a, const vec4& b)
570
{ return vec4(a.n[VX] * b.n[VX], a.n[VY] * b.n[VY], a.n[VZ] * b.n[VZ],
571
a.n[VW] * b.n[VW]); }
572
573
574
/****************************************************************
575
* *
576
* mat3 member functions *
577
* *
578
****************************************************************/
579
580
// CONSTRUCTORS
581
582
mat3::mat3(void) { *this = identity2D(); }
583
584
mat3::mat3(const vec3& v0, const vec3& v1, const vec3& v2)
585
{ this->set( v0, v1, v2 ); };
586
587
mat3::mat3(const float d)
588
{ v[0] = v[1] = v[2] = vec3(d); }
589
590
mat3::mat3(const mat3& m)
591
{ v[0] = m.v[0]; v[1] = m.v[1]; v[2] = m.v[2]; }
592
593
594
// ASSIGNMENT OPERATORS
595
596
mat3& mat3::operator = ( const mat3& m )
597
{ v[0] = m.v[0]; v[1] = m.v[1]; v[2] = m.v[2]; return *this; }
598
599
mat3& mat3::operator += ( const mat3& m )
600
{ v[0] += m.v[0]; v[1] += m.v[1]; v[2] += m.v[2]; return *this; }
601
602
mat3& mat3::operator -= ( const mat3& m )
603
{ v[0] -= m.v[0]; v[1] -= m.v[1]; v[2] -= m.v[2]; return *this; }
604
605
mat3& mat3::operator *= ( const float d )
606
{ v[0] *= d; v[1] *= d; v[2] *= d; return *this; }
607
608
mat3& mat3::operator /= ( const float d )
609
{ v[0] /= d; v[1] /= d; v[2] /= d; return *this; }
610
611
vec3& mat3::operator [] ( int i) {
612
if (i < VX || i > VZ)
613
//VEC_ERROR("mat3 [] operator: illegal access; index = " << i << '\n')
614
VEC_ERROR("mat3 [] operator: illegal access" );
615
return v[i];
616
}
617
618
void mat3::set( const vec3& v0, const vec3& v1, const vec3& v2 ) {
619
v[0] = v0; v[1] = v1; v[2] = v2;
620
}
621
622
// SPECIAL FUNCTIONS
623
624
mat3 mat3::transpose(void) {
625
return mat3(vec3(v[0][0], v[1][0], v[2][0]),
626
vec3(v[0][1], v[1][1], v[2][1]),
627
vec3(v[0][2], v[1][2], v[2][2]));
628
}
629
630
mat3 mat3::inverse(void) // Gauss-Jordan elimination with partial pivoting
631
{
632
mat3 a(*this), // As a evolves from original mat into identity
633
b(identity2D()); // b evolves from identity into inverse(a)
634
int i, j, i1;
635
636
// Loop over cols of a from left to right, eliminating above and below diag
637
for (j=0; j<3; j++) { // Find largest pivot in column j among rows j..2
638
i1 = j; // Row with largest pivot candidate
639
for (i=j+1; i<3; i++)
640
if (fabs(a.v[i].n[j]) > fabs(a.v[i1].n[j]))
641
i1 = i;
642
643
// Swap rows i1 and j in a and b to put pivot on diagonal
644
swap(a.v[i1], a.v[j]);
645
swap(b.v[i1], b.v[j]);
646
647
// Scale row j to have a unit diagonal
648
if (a.v[j].n[j]==0.)
649
VEC_ERROR("mat3::inverse: singular matrix; can't invert\n");
650
b.v[j] /= a.v[j].n[j];
651
a.v[j] /= a.v[j].n[j];
652
653
// Eliminate off-diagonal elems in col j of a, doing identical ops to b
654
for (i=0; i<3; i++)
655
if (i!=j) {
656
b.v[i] -= a.v[i].n[j]*b.v[j];
657
a.v[i] -= a.v[i].n[j]*a.v[j];
658
}
659
}
660
return b;
661
}
662
663
mat3& mat3::apply(V_FCT_PTR fct) {
664
v[VX].apply(fct);
665
v[VY].apply(fct);
666
v[VZ].apply(fct);
667
return *this;
668
}
669
670
671
// FRIENDS
672
673
mat3 operator - (const mat3& a)
674
{ return mat3(-a.v[0], -a.v[1], -a.v[2]); }
675
676
mat3 operator + (const mat3& a, const mat3& b)
677
{ return mat3(a.v[0] + b.v[0], a.v[1] + b.v[1], a.v[2] + b.v[2]); }
678
679
mat3 operator - (const mat3& a, const mat3& b)
680
{ return mat3(a.v[0] - b.v[0], a.v[1] - b.v[1], a.v[2] - b.v[2]); }
681
682
mat3 operator * (mat3& a, mat3& b) {
683
#define ROWCOL(i, j) \
684
a.v[i].n[0]*b.v[0][j] + a.v[i].n[1]*b.v[1][j] + a.v[i].n[2]*b.v[2][j]
685
return mat3(vec3(ROWCOL(0,0), ROWCOL(0,1), ROWCOL(0,2)),
686
vec3(ROWCOL(1,0), ROWCOL(1,1), ROWCOL(1,2)),
687
vec3(ROWCOL(2,0), ROWCOL(2,1), ROWCOL(2,2)));
688
#undef ROWCOL
689
}
690
691
mat3 operator * (const mat3& a, const float d)
692
{ return mat3(a.v[0] * d, a.v[1] * d, a.v[2] * d); }
693
694
mat3 operator * (const float d, const mat3& a)
695
{ return a*d; }
696
697
mat3 operator / (const mat3& a, const float d)
698
{ return mat3(a.v[0] / d, a.v[1] / d, a.v[2] / d); }
699
700
int operator == (const mat3& a, const mat3& b)
701
{ return (a.v[0] == b.v[0]) && (a.v[1] == b.v[1]) && (a.v[2] == b.v[2]); }
702
703
int operator != (const mat3& a, const mat3& b)
704
{ return !(a == b); }
705
706
/*ostream& operator << (ostream& s, mat3& m)
707
{ return s << m.v[VX] << '\n' << m.v[VY] << '\n' << m.v[VZ]; }
708
709
istream& operator >> (istream& s, mat3& m) {
710
mat3 m_tmp;
711
712
s >> m_tmp[VX] >> m_tmp[VY] >> m_tmp[VZ];
713
if (s)
714
m = m_tmp;
715
return s;
716
}
717
*/
718
719
void swap(mat3& a, mat3& b)
720
{ mat3 tmp(a); a = b; b = tmp; }
721
722
void mat3::print( FILE *file, char *name )
723
{
724
int i, j;
725
726
fprintf( stderr, "%s:\n", name );
727
728
for( i = 0; i < 3; i++ )
729
{
730
fprintf( stderr, " " );
731
for( j = 0; j < 3; j++ )
732
{
733
fprintf( stderr, "%f ", v[i][j] );
734
}
735
fprintf( stderr, "\n" );
736
}
737
}
738
739
740
741
/****************************************************************
742
* *
743
* mat4 member functions *
744
* *
745
****************************************************************/
746
747
// CONSTRUCTORS
748
749
mat4::mat4(void) { *this = identity3D();}
750
751
mat4::mat4(const vec4& v0, const vec4& v1, const vec4& v2, const vec4& v3)
752
{ v[0] = v0; v[1] = v1; v[2] = v2; v[3] = v3; }
753
754
mat4::mat4(const float d)
755
{ v[0] = v[1] = v[2] = v[3] = vec4(d); }
756
757
mat4::mat4(const mat4& m)
758
{ v[0] = m.v[0]; v[1] = m.v[1]; v[2] = m.v[2]; v[3] = m.v[3]; }
759
760
mat4::mat4(const float a00, const float a01, const float a02, const float a03,
761
const float a10, const float a11, const float a12, const float a13,
762
const float a20, const float a21, const float a22, const float a23,
763
const float a30, const float a31, const float a32, const float a33 )
764
{
765
v[0][0] = a00; v[0][1] = a01; v[0][2] = a02; v[0][3] = a03;
766
v[1][0] = a10; v[1][1] = a11; v[1][2] = a12; v[1][3] = a13;
767
v[2][0] = a20; v[2][1] = a21; v[2][2] = a22; v[2][3] = a23;
768
v[3][0] = a30; v[3][1] = a31; v[3][2] = a32; v[3][3] = a33;
769
}
770
771
// ASSIGNMENT OPERATORS
772
773
mat4& mat4::operator = ( const mat4& m )
774
{ v[0] = m.v[0]; v[1] = m.v[1]; v[2] = m.v[2]; v[3] = m.v[3];
775
return *this; }
776
777
mat4& mat4::operator += ( const mat4& m )
778
{ v[0] += m.v[0]; v[1] += m.v[1]; v[2] += m.v[2]; v[3] += m.v[3];
779
return *this; }
780
781
mat4& mat4::operator -= ( const mat4& m )
782
{ v[0] -= m.v[0]; v[1] -= m.v[1]; v[2] -= m.v[2]; v[3] -= m.v[3];
783
return *this; }
784
785
mat4& mat4::operator *= ( const float d )
786
{ v[0] *= d; v[1] *= d; v[2] *= d; v[3] *= d; return *this; }
787
788
mat4& mat4::operator /= ( const float d )
789
{ v[0] /= d; v[1] /= d; v[2] /= d; v[3] /= d; return *this; }
790
791
vec4& mat4::operator [] ( int i) {
792
if (i < VX || i > VW)
793
//VEC_ERROR("mat4 [] operator: illegal access; index = " << i << '\n')
794
VEC_ERROR("mat4 [] operator: illegal access" );
795
return v[i];
796
}
797
798
// SPECIAL FUNCTIONS;
799
800
mat4 mat4::transpose(void) {
801
return mat4(vec4(v[0][0], v[1][0], v[2][0], v[3][0]),
802
vec4(v[0][1], v[1][1], v[2][1], v[3][1]),
803
vec4(v[0][2], v[1][2], v[2][2], v[3][2]),
804
vec4(v[0][3], v[1][3], v[2][3], v[3][3]));
805
}
806
807
mat4 mat4::inverse(void) // Gauss-Jordan elimination with partial pivoting
808
{
809
mat4 a(*this), // As a evolves from original mat into identity
810
b(identity3D()); // b evolves from identity into inverse(a)
811
int i, j, i1;
812
813
// Loop over cols of a from left to right, eliminating above and below diag
814
for (j=0; j<4; j++) { // Find largest pivot in column j among rows j..3
815
i1 = j; // Row with largest pivot candidate
816
for (i=j+1; i<4; i++)
817
if (fabs(a.v[i].n[j]) > fabs(a.v[i1].n[j]))
818
i1 = i;
819
820
// Swap rows i1 and j in a and b to put pivot on diagonal
821
swap(a.v[i1], a.v[j]);
822
swap(b.v[i1], b.v[j]);
823
824
// Scale row j to have a unit diagonal
825
if (a.v[j].n[j]==0.)
826
VEC_ERROR("mat4::inverse: singular matrix; can't invert\n");
827
b.v[j] /= a.v[j].n[j];
828
a.v[j] /= a.v[j].n[j];
829
830
// Eliminate off-diagonal elems in col j of a, doing identical ops to b
831
for (i=0; i<4; i++)
832
if (i!=j) {
833
b.v[i] -= a.v[i].n[j]*b.v[j];
834
a.v[i] -= a.v[i].n[j]*a.v[j];
835
}
836
}
837
return b;
838
}
839
840
mat4& mat4::apply(V_FCT_PTR fct)
841
{ v[VX].apply(fct); v[VY].apply(fct); v[VZ].apply(fct); v[VW].apply(fct);
842
return *this; }
843
844
845
void mat4::print( FILE *file, char *name )
846
{
847
int i, j;
848
849
fprintf( stderr, "%s:\n", name );
850
851
for( i = 0; i < 4; i++ )
852
{
853
fprintf( stderr, " " );
854
for( j = 0; j < 4; j++ )
855
{
856
fprintf( stderr, "%f ", v[i][j] );
857
}
858
fprintf( stderr, "\n" );
859
}
860
}
861
862
void mat4::swap_rows( int i, int j )
863
{
864
vec4 t;
865
866
t = v[i];
867
v[i] = v[j];
868
v[j] = t;
869
}
870
871
void mat4::swap_cols( int i, int j )
872
{
873
float t;
874
int k;
875
876
for(k=0; k<4; k++ ) {
877
t = v[k][i];
878
v[k][i] = v[k][j];
879
v[k][j] = t;
880
}
881
}
882
883
884
// FRIENDS
885
886
mat4 operator - (const mat4& a)
887
{ return mat4(-a.v[0], -a.v[1], -a.v[2], -a.v[3]); }
888
889
mat4 operator + (const mat4& a, const mat4& b)
890
{ return mat4(a.v[0] + b.v[0], a.v[1] + b.v[1], a.v[2] + b.v[2],
891
a.v[3] + b.v[3]);
892
}
893
894
mat4 operator - (const mat4& a, const mat4& b)
895
{ return mat4(a.v[0] - b.v[0], a.v[1] - b.v[1], a.v[2] - b.v[2], a.v[3] - b.v[3]); }
896
897
mat4 operator * (mat4& a, mat4& b) {
898
#define ROWCOL(i, j) a.v[i].n[0]*b.v[0][j] + a.v[i].n[1]*b.v[1][j] + \
899
a.v[i].n[2]*b.v[2][j] + a.v[i].n[3]*b.v[3][j]
900
return mat4(
901
vec4(ROWCOL(0,0), ROWCOL(0,1), ROWCOL(0,2), ROWCOL(0,3)),
902
vec4(ROWCOL(1,0), ROWCOL(1,1), ROWCOL(1,2), ROWCOL(1,3)),
903
vec4(ROWCOL(2,0), ROWCOL(2,1), ROWCOL(2,2), ROWCOL(2,3)),
904
vec4(ROWCOL(3,0), ROWCOL(3,1), ROWCOL(3,2), ROWCOL(3,3))
905
);
906
}
907
908
mat4 operator * (const mat4& a, const float d)
909
{ return mat4(a.v[0] * d, a.v[1] * d, a.v[2] * d, a.v[3] * d); }
910
911
mat4 operator * (const float d, const mat4& a)
912
{ return a*d; }
913
914
mat4 operator / (const mat4& a, const float d)
915
{ return mat4(a.v[0] / d, a.v[1] / d, a.v[2] / d, a.v[3] / d); }
916
917
int operator == (const mat4& a, const mat4& b)
918
{ return ((a.v[0] == b.v[0]) && (a.v[1] == b.v[1]) && (a.v[2] == b.v[2]) &&
919
(a.v[3] == b.v[3])); }
920
921
int operator != (const mat4& a, const mat4& b)
922
{ return !(a == b); }
923
924
/*ostream& operator << (ostream& s, mat4& m)
925
{ return s << m.v[VX] << '\n' << m.v[VY] << '\n' << m.v[VZ] << '\n' << m.v[VW]; }
926
927
istream& operator >> (istream& s, mat4& m)
928
{
929
mat4 m_tmp;
930
931
s >> m_tmp[VX] >> m_tmp[VY] >> m_tmp[VZ] >> m_tmp[VW];
932
if (s)
933
m = m_tmp;
934
return s;
935
}
936
*/
937
void swap(mat4& a, mat4& b)
938
{ mat4 tmp(a); a = b; b = tmp; }
939
940
941
/****************************************************************
942
* *
943
* 2D functions and 3D functions *
944
* *
945
****************************************************************/
946
947
mat3 identity2D(void)
948
{ return mat3(vec3(1.0, 0.0, 0.0),
949
vec3(0.0, 1.0, 0.0),
950
vec3(0.0, 0.0, 1.0)); }
951
952
mat3 translation2D(vec2& v)
953
{ return mat3(vec3(1.0, 0.0, v[VX]),
954
vec3(0.0, 1.0, v[VY]),
955
vec3(0.0, 0.0, 1.0)); }
956
957
mat3 rotation2D(vec2& Center, const float angleDeg) {
958
float angleRad = angleDeg * M_PI / 180.0,
959
c = cos(angleRad),
960
s = sin(angleRad);
961
962
return mat3(vec3(c, -s, Center[VX] * (1.0-c) + Center[VY] * s),
963
vec3(s, c, Center[VY] * (1.0-c) - Center[VX] * s),
964
vec3(0.0, 0.0, 1.0));
965
}
966
967
mat3 scaling2D(vec2& scaleVector)
968
{ return mat3(vec3(scaleVector[VX], 0.0, 0.0),
969
vec3(0.0, scaleVector[VY], 0.0),
970
vec3(0.0, 0.0, 1.0)); }
971
972
mat4 identity3D(void)
973
{ return mat4(vec4(1.0, 0.0, 0.0, 0.0),
974
vec4(0.0, 1.0, 0.0, 0.0),
975
vec4(0.0, 0.0, 1.0, 0.0),
976
vec4(0.0, 0.0, 0.0, 1.0)); }
977
978
mat4 translation3D(vec3& v)
979
{ return mat4(vec4(1.0, 0.0, 0.0, v[VX]),
980
vec4(0.0, 1.0, 0.0, v[VY]),
981
vec4(0.0, 0.0, 1.0, v[VZ]),
982
vec4(0.0, 0.0, 0.0, 1.0)); }
983
984
mat4 rotation3D(vec3& Axis, const float angleDeg) {
985
float angleRad = angleDeg * M_PI / 180.0,
986
c = cos(angleRad),
987
s = sin(angleRad),
988
t = 1.0 - c;
989
990
Axis.normalize();
991
return mat4(vec4(t * Axis[VX] * Axis[VX] + c,
992
t * Axis[VX] * Axis[VY] - s * Axis[VZ],
993
t * Axis[VX] * Axis[VZ] + s * Axis[VY],
994
0.0),
995
vec4(t * Axis[VX] * Axis[VY] + s * Axis[VZ],
996
t * Axis[VY] * Axis[VY] + c,
997
t * Axis[VY] * Axis[VZ] - s * Axis[VX],
998
0.0),
999
vec4(t * Axis[VX] * Axis[VZ] - s * Axis[VY],
1000
t * Axis[VY] * Axis[VZ] + s * Axis[VX],
1001
t * Axis[VZ] * Axis[VZ] + c,
1002
0.0),
1003
vec4(0.0, 0.0, 0.0, 1.0));
1004
}
1005
1006
mat4 rotation3Drad(vec3& Axis, const float angleRad) {
1007
float c = cos(angleRad),
1008
s = sin(angleRad),
1009
t = 1.0 - c;
1010
1011
Axis.normalize();
1012
return mat4(vec4(t * Axis[VX] * Axis[VX] + c,
1013
t * Axis[VX] * Axis[VY] - s * Axis[VZ],
1014
t * Axis[VX] * Axis[VZ] + s * Axis[VY],
1015
0.0),
1016
vec4(t * Axis[VX] * Axis[VY] + s * Axis[VZ],
1017
t * Axis[VY] * Axis[VY] + c,
1018
t * Axis[VY] * Axis[VZ] - s * Axis[VX],
1019
0.0),
1020
vec4(t * Axis[VX] * Axis[VZ] - s * Axis[VY],
1021
t * Axis[VY] * Axis[VZ] + s * Axis[VX],
1022
t * Axis[VZ] * Axis[VZ] + c,
1023
0.0),
1024
vec4(0.0, 0.0, 0.0, 1.0));
1025
}
1026
1027
mat4 scaling3D(vec3& scaleVector)
1028
{ return mat4(vec4(scaleVector[VX], 0.0, 0.0, 0.0),
1029
vec4(0.0, scaleVector[VY], 0.0, 0.0),
1030
vec4(0.0, 0.0, scaleVector[VZ], 0.0),
1031
vec4(0.0, 0.0, 0.0, 1.0)); }
1032
1033
mat4 perspective3D(const float d)
1034
{ return mat4(vec4(1.0, 0.0, 0.0, 0.0),
1035
vec4(0.0, 1.0, 0.0, 0.0),
1036
vec4(0.0, 0.0, 1.0, 0.0),
1037
vec4(0.0, 0.0, 1.0/d, 0.0)); }
Loggerhead is a web-based interface for Breezy
Version: 2.0.1