dungeon_crawler

annotate prototype/vmath/matrix.cc @ 1:96de911d05d4

find changesets by author, revision, files, or words in the commit message
started a rough prototype
author John Tsiombikas <nuclear@mutantstargoat.com>
date Thu, 28 Jun 2012 06:05:50 +0300
parents
children
rev   line source
nuclear@1 1 #include <cstdio>
nuclear@1 2 #include <cmath>
nuclear@1 3 #include "matrix.h"
nuclear@1 4 #include "vector.h"
nuclear@1 5 #include "quat.h"
nuclear@1 6
nuclear@1 7 using namespace std;
nuclear@1 8
nuclear@1 9 // ----------- Matrix3x3 --------------
nuclear@1 10
nuclear@1 11 Matrix3x3 Matrix3x3::identity = Matrix3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
nuclear@1 12
nuclear@1 13 Matrix3x3::Matrix3x3()
nuclear@1 14 {
nuclear@1 15 *this = identity;
nuclear@1 16 }
nuclear@1 17
nuclear@1 18 Matrix3x3::Matrix3x3( scalar_t m11, scalar_t m12, scalar_t m13,
nuclear@1 19 scalar_t m21, scalar_t m22, scalar_t m23,
nuclear@1 20 scalar_t m31, scalar_t m32, scalar_t m33)
nuclear@1 21 {
nuclear@1 22 m[0][0] = m11; m[0][1] = m12; m[0][2] = m13;
nuclear@1 23 m[1][0] = m21; m[1][1] = m22; m[1][2] = m23;
nuclear@1 24 m[2][0] = m31; m[2][1] = m32; m[2][2] = m33;
nuclear@1 25 }
nuclear@1 26
nuclear@1 27 Matrix3x3::Matrix3x3(const Vector3 &ivec, const Vector3 &jvec, const Vector3 &kvec)
nuclear@1 28 {
nuclear@1 29 set_row_vector(ivec, 0);
nuclear@1 30 set_row_vector(jvec, 1);
nuclear@1 31 set_row_vector(kvec, 2);
nuclear@1 32 }
nuclear@1 33
nuclear@1 34 Matrix3x3::Matrix3x3(const mat3_t cmat)
nuclear@1 35 {
nuclear@1 36 memcpy(m, cmat, sizeof(mat3_t));
nuclear@1 37 }
nuclear@1 38
nuclear@1 39 Matrix3x3::Matrix3x3(const Matrix4x4 &mat4x4)
nuclear@1 40 {
nuclear@1 41 for(int i=0; i<3; i++) {
nuclear@1 42 for(int j=0; j<3; j++) {
nuclear@1 43 m[i][j] = mat4x4[i][j];
nuclear@1 44 }
nuclear@1 45 }
nuclear@1 46 }
nuclear@1 47
nuclear@1 48 Matrix3x3 operator +(const Matrix3x3 &m1, const Matrix3x3 &m2)
nuclear@1 49 {
nuclear@1 50 Matrix3x3 res;
nuclear@1 51 const scalar_t *op1 = m1.m[0], *op2 = m2.m[0];
nuclear@1 52 scalar_t *dest = res.m[0];
nuclear@1 53
nuclear@1 54 for(int i=0; i<9; i++) {
nuclear@1 55 *dest++ = *op1++ + *op2++;
nuclear@1 56 }
nuclear@1 57 return res;
nuclear@1 58 }
nuclear@1 59
nuclear@1 60 Matrix3x3 operator -(const Matrix3x3 &m1, const Matrix3x3 &m2)
nuclear@1 61 {
nuclear@1 62 Matrix3x3 res;
nuclear@1 63 const scalar_t *op1 = m1.m[0], *op2 = m2.m[0];
nuclear@1 64 scalar_t *dest = res.m[0];
nuclear@1 65
nuclear@1 66 for(int i=0; i<9; i++) {
nuclear@1 67 *dest++ = *op1++ - *op2++;
nuclear@1 68 }
nuclear@1 69 return res;
nuclear@1 70 }
nuclear@1 71
nuclear@1 72 Matrix3x3 operator *(const Matrix3x3 &m1, const Matrix3x3 &m2)
nuclear@1 73 {
nuclear@1 74 Matrix3x3 res;
nuclear@1 75 for(int i=0; i<3; i++) {
nuclear@1 76 for(int j=0; j<3; j++) {
nuclear@1 77 res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j];
nuclear@1 78 }
nuclear@1 79 }
nuclear@1 80 return res;
nuclear@1 81 }
nuclear@1 82
nuclear@1 83 void operator +=(Matrix3x3 &m1, const Matrix3x3 &m2)
nuclear@1 84 {
nuclear@1 85 scalar_t *op1 = m1.m[0];
nuclear@1 86 const scalar_t *op2 = m2.m[0];
nuclear@1 87
nuclear@1 88 for(int i=0; i<9; i++) {
nuclear@1 89 *op1++ += *op2++;
nuclear@1 90 }
nuclear@1 91 }
nuclear@1 92
nuclear@1 93 void operator -=(Matrix3x3 &m1, const Matrix3x3 &m2)
nuclear@1 94 {
nuclear@1 95 scalar_t *op1 = m1.m[0];
nuclear@1 96 const scalar_t *op2 = m2.m[0];
nuclear@1 97
nuclear@1 98 for(int i=0; i<9; i++) {
nuclear@1 99 *op1++ -= *op2++;
nuclear@1 100 }
nuclear@1 101 }
nuclear@1 102
nuclear@1 103 void operator *=(Matrix3x3 &m1, const Matrix3x3 &m2)
nuclear@1 104 {
nuclear@1 105 Matrix3x3 res;
nuclear@1 106 for(int i=0; i<3; i++) {
nuclear@1 107 for(int j=0; j<3; j++) {
nuclear@1 108 res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j];
nuclear@1 109 }
nuclear@1 110 }
nuclear@1 111 memcpy(m1.m, res.m, 9 * sizeof(scalar_t));
nuclear@1 112 }
nuclear@1 113
nuclear@1 114 Matrix3x3 operator *(const Matrix3x3 &mat, scalar_t scalar)
nuclear@1 115 {
nuclear@1 116 Matrix3x3 res;
nuclear@1 117 const scalar_t *mptr = mat.m[0];
nuclear@1 118 scalar_t *dptr = res.m[0];
nuclear@1 119
nuclear@1 120 for(int i=0; i<9; i++) {
nuclear@1 121 *dptr++ = *mptr++ * scalar;
nuclear@1 122 }
nuclear@1 123 return res;
nuclear@1 124 }
nuclear@1 125
nuclear@1 126 Matrix3x3 operator *(scalar_t scalar, const Matrix3x3 &mat)
nuclear@1 127 {
nuclear@1 128 Matrix3x3 res;
nuclear@1 129 const scalar_t *mptr = mat.m[0];
nuclear@1 130 scalar_t *dptr = res.m[0];
nuclear@1 131
nuclear@1 132 for(int i=0; i<9; i++) {
nuclear@1 133 *dptr++ = *mptr++ * scalar;
nuclear@1 134 }
nuclear@1 135 return res;
nuclear@1 136 }
nuclear@1 137
nuclear@1 138 void operator *=(Matrix3x3 &mat, scalar_t scalar)
nuclear@1 139 {
nuclear@1 140 scalar_t *mptr = mat.m[0];
nuclear@1 141
nuclear@1 142 for(int i=0; i<9; i++) {
nuclear@1 143 *mptr++ *= scalar;
nuclear@1 144 }
nuclear@1 145 }
nuclear@1 146
nuclear@1 147 void Matrix3x3::translate(const Vector2 &trans)
nuclear@1 148 {
nuclear@1 149 Matrix3x3 tmat(1, 0, trans.x, 0, 1, trans.y, 0, 0, 1);
nuclear@1 150 *this *= tmat;
nuclear@1 151 }
nuclear@1 152
nuclear@1 153 void Matrix3x3::set_translation(const Vector2 &trans)
nuclear@1 154 {
nuclear@1 155 *this = Matrix3x3(1, 0, trans.x, 0, 1, trans.y, 0, 0, 1);
nuclear@1 156 }
nuclear@1 157
nuclear@1 158 void Matrix3x3::rotate(scalar_t angle)
nuclear@1 159 {
nuclear@1 160 scalar_t cos_a = cos(angle);
nuclear@1 161 scalar_t sin_a = sin(angle);
nuclear@1 162 Matrix3x3 rmat( cos_a, -sin_a, 0,
nuclear@1 163 sin_a, cos_a, 0,
nuclear@1 164 0, 0, 1);
nuclear@1 165 *this *= rmat;
nuclear@1 166 }
nuclear@1 167
nuclear@1 168 void Matrix3x3::set_rotation(scalar_t angle)
nuclear@1 169 {
nuclear@1 170 scalar_t cos_a = cos(angle);
nuclear@1 171 scalar_t sin_a = sin(angle);
nuclear@1 172 *this = Matrix3x3(cos_a, -sin_a, 0, sin_a, cos_a, 0, 0, 0, 1);
nuclear@1 173 }
nuclear@1 174
nuclear@1 175 void Matrix3x3::rotate(const Vector3 &euler_angles)
nuclear@1 176 {
nuclear@1 177 Matrix3x3 xrot, yrot, zrot;
nuclear@1 178
nuclear@1 179 xrot = Matrix3x3( 1, 0, 0,
nuclear@1 180 0, cos(euler_angles.x), -sin(euler_angles.x),
nuclear@1 181 0, sin(euler_angles.x), cos(euler_angles.x));
nuclear@1 182
nuclear@1 183 yrot = Matrix3x3( cos(euler_angles.y), 0, sin(euler_angles.y),
nuclear@1 184 0, 1, 0,
nuclear@1 185 -sin(euler_angles.y), 0, cos(euler_angles.y));
nuclear@1 186
nuclear@1 187 zrot = Matrix3x3( cos(euler_angles.z), -sin(euler_angles.z), 0,
nuclear@1 188 sin(euler_angles.z), cos(euler_angles.z), 0,
nuclear@1 189 0, 0, 1);
nuclear@1 190
nuclear@1 191 *this *= xrot * yrot * zrot;
nuclear@1 192 }
nuclear@1 193
nuclear@1 194 void Matrix3x3::set_rotation(const Vector3 &euler_angles)
nuclear@1 195 {
nuclear@1 196 Matrix3x3 xrot, yrot, zrot;
nuclear@1 197
nuclear@1 198 xrot = Matrix3x3( 1, 0, 0,
nuclear@1 199 0, cos(euler_angles.x), -sin(euler_angles.x),
nuclear@1 200 0, sin(euler_angles.x), cos(euler_angles.x));
nuclear@1 201
nuclear@1 202 yrot = Matrix3x3( cos(euler_angles.y), 0, sin(euler_angles.y),
nuclear@1 203 0, 1, 0,
nuclear@1 204 -sin(euler_angles.y), 0, cos(euler_angles.y));
nuclear@1 205
nuclear@1 206 zrot = Matrix3x3( cos(euler_angles.z), -sin(euler_angles.z), 0,
nuclear@1 207 sin(euler_angles.z), cos(euler_angles.z), 0,
nuclear@1 208 0, 0, 1);
nuclear@1 209
nuclear@1 210 *this = xrot * yrot * zrot;
nuclear@1 211 }
nuclear@1 212
nuclear@1 213 void Matrix3x3::rotate(const Vector3 &axis, scalar_t angle)
nuclear@1 214 {
nuclear@1 215 scalar_t sina = (scalar_t)sin(angle);
nuclear@1 216 scalar_t cosa = (scalar_t)cos(angle);
nuclear@1 217 scalar_t invcosa = 1-cosa;
nuclear@1 218 scalar_t nxsq = axis.x * axis.x;
nuclear@1 219 scalar_t nysq = axis.y * axis.y;
nuclear@1 220 scalar_t nzsq = axis.z * axis.z;
nuclear@1 221
nuclear@1 222 Matrix3x3 xform;
nuclear@1 223 xform.m[0][0] = nxsq + (1-nxsq) * cosa;
nuclear@1 224 xform.m[0][1] = axis.x * axis.y * invcosa - axis.z * sina;
nuclear@1 225 xform.m[0][2] = axis.x * axis.z * invcosa + axis.y * sina;
nuclear@1 226
nuclear@1 227 xform.m[1][0] = axis.x * axis.y * invcosa + axis.z * sina;
nuclear@1 228 xform.m[1][1] = nysq + (1-nysq) * cosa;
nuclear@1 229 xform.m[1][2] = axis.y * axis.z * invcosa - axis.x * sina;
nuclear@1 230
nuclear@1 231 xform.m[2][0] = axis.x * axis.z * invcosa - axis.y * sina;
nuclear@1 232 xform.m[2][1] = axis.y * axis.z * invcosa + axis.x * sina;
nuclear@1 233 xform.m[2][2] = nzsq + (1-nzsq) * cosa;
nuclear@1 234
nuclear@1 235 *this *= xform;
nuclear@1 236 }
nuclear@1 237
nuclear@1 238 void Matrix3x3::set_rotation(const Vector3 &axis, scalar_t angle)
nuclear@1 239 {
nuclear@1 240 scalar_t sina = (scalar_t)sin(angle);
nuclear@1 241 scalar_t cosa = (scalar_t)cos(angle);
nuclear@1 242 scalar_t invcosa = 1-cosa;
nuclear@1 243 scalar_t nxsq = axis.x * axis.x;
nuclear@1 244 scalar_t nysq = axis.y * axis.y;
nuclear@1 245 scalar_t nzsq = axis.z * axis.z;
nuclear@1 246
nuclear@1 247 reset_identity();
nuclear@1 248 m[0][0] = nxsq + (1-nxsq) * cosa;
nuclear@1 249 m[0][1] = axis.x * axis.y * invcosa - axis.z * sina;
nuclear@1 250 m[0][2] = axis.x * axis.z * invcosa + axis.y * sina;
nuclear@1 251 m[1][0] = axis.x * axis.y * invcosa + axis.z * sina;
nuclear@1 252 m[1][1] = nysq + (1-nysq) * cosa;
nuclear@1 253 m[1][2] = axis.y * axis.z * invcosa - axis.x * sina;
nuclear@1 254 m[2][0] = axis.x * axis.z * invcosa - axis.y * sina;
nuclear@1 255 m[2][1] = axis.y * axis.z * invcosa + axis.x * sina;
nuclear@1 256 m[2][2] = nzsq + (1-nzsq) * cosa;
nuclear@1 257 }
nuclear@1 258
nuclear@1 259 void Matrix3x3::scale(const Vector3 &scale_vec)
nuclear@1 260 {
nuclear@1 261 Matrix3x3 smat( scale_vec.x, 0, 0,
nuclear@1 262 0, scale_vec.y, 0,
nuclear@1 263 0, 0, scale_vec.z);
nuclear@1 264 *this *= smat;
nuclear@1 265 }
nuclear@1 266
nuclear@1 267 void Matrix3x3::set_scaling(const Vector3 &scale_vec)
nuclear@1 268 {
nuclear@1 269 *this = Matrix3x3( scale_vec.x, 0, 0,
nuclear@1 270 0, scale_vec.y, 0,
nuclear@1 271 0, 0, scale_vec.z);
nuclear@1 272 }
nuclear@1 273
nuclear@1 274 void Matrix3x3::set_column_vector(const Vector3 &vec, unsigned int col_index)
nuclear@1 275 {
nuclear@1 276 m[0][col_index] = vec.x;
nuclear@1 277 m[1][col_index] = vec.y;
nuclear@1 278 m[2][col_index] = vec.z;
nuclear@1 279 }
nuclear@1 280
nuclear@1 281 void Matrix3x3::set_row_vector(const Vector3 &vec, unsigned int row_index)
nuclear@1 282 {
nuclear@1 283 m[row_index][0] = vec.x;
nuclear@1 284 m[row_index][1] = vec.y;
nuclear@1 285 m[row_index][2] = vec.z;
nuclear@1 286 }
nuclear@1 287
nuclear@1 288 Vector3 Matrix3x3::get_column_vector(unsigned int col_index) const
nuclear@1 289 {
nuclear@1 290 return Vector3(m[0][col_index], m[1][col_index], m[2][col_index]);
nuclear@1 291 }
nuclear@1 292
nuclear@1 293 Vector3 Matrix3x3::get_row_vector(unsigned int row_index) const
nuclear@1 294 {
nuclear@1 295 return Vector3(m[row_index][0], m[row_index][1], m[row_index][2]);
nuclear@1 296 }
nuclear@1 297
nuclear@1 298 void Matrix3x3::transpose()
nuclear@1 299 {
nuclear@1 300 Matrix3x3 tmp = *this;
nuclear@1 301 for(int i=0; i<3; i++) {
nuclear@1 302 for(int j=0; j<3; j++) {
nuclear@1 303 m[i][j] = tmp[j][i];
nuclear@1 304 }
nuclear@1 305 }
nuclear@1 306 }
nuclear@1 307
nuclear@1 308 Matrix3x3 Matrix3x3::transposed() const
nuclear@1 309 {
nuclear@1 310 Matrix3x3 res;
nuclear@1 311 for(int i=0; i<3; i++) {
nuclear@1 312 for(int j=0; j<3; j++) {
nuclear@1 313 res[i][j] = m[j][i];
nuclear@1 314 }
nuclear@1 315 }
nuclear@1 316 return res;
nuclear@1 317 }
nuclear@1 318
nuclear@1 319 scalar_t Matrix3x3::determinant() const
nuclear@1 320 {
nuclear@1 321 return m[0][0] * (m[1][1]*m[2][2] - m[1][2]*m[2][1]) -
nuclear@1 322 m[0][1] * (m[1][0]*m[2][2] - m[1][2]*m[2][0]) +
nuclear@1 323 m[0][2] * (m[1][0]*m[2][1] - m[1][1]*m[2][0]);
nuclear@1 324 }
nuclear@1 325
nuclear@1 326 Matrix3x3 Matrix3x3::inverse() const
nuclear@1 327 {
nuclear@1 328 // TODO: implement 3x3 inverse
nuclear@1 329 return *this;
nuclear@1 330 }
nuclear@1 331
nuclear@1 332 ostream &operator <<(ostream &out, const Matrix3x3 &mat)
nuclear@1 333 {
nuclear@1 334 for(int i=0; i<3; i++) {
nuclear@1 335 char str[100];
nuclear@1 336 sprintf(str, "[ %12.5f %12.5f %12.5f ]\n", (float)mat.m[i][0], (float)mat.m[i][1], (float)mat.m[i][2]);
nuclear@1 337 out << str;
nuclear@1 338 }
nuclear@1 339 return out;
nuclear@1 340 }
nuclear@1 341
nuclear@1 342
nuclear@1 343
nuclear@1 344 /* ----------------- Matrix4x4 implementation --------------- */
nuclear@1 345
nuclear@1 346 Matrix4x4 Matrix4x4::identity = Matrix4x4(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
nuclear@1 347
nuclear@1 348 Matrix4x4::Matrix4x4()
nuclear@1 349 {
nuclear@1 350 *this = identity;
nuclear@1 351 }
nuclear@1 352
nuclear@1 353 Matrix4x4::Matrix4x4( scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14,
nuclear@1 354 scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24,
nuclear@1 355 scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34,
nuclear@1 356 scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44)
nuclear@1 357 {
nuclear@1 358 m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; m[0][3] = m14;
nuclear@1 359 m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; m[1][3] = m24;
nuclear@1 360 m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; m[2][3] = m34;
nuclear@1 361 m[3][0] = m41; m[3][1] = m42; m[3][2] = m43; m[3][3] = m44;
nuclear@1 362 }
nuclear@1 363
nuclear@1 364 Matrix4x4::Matrix4x4(const mat4_t cmat)
nuclear@1 365 {
nuclear@1 366 memcpy(m, cmat, sizeof(mat4_t));
nuclear@1 367 }
nuclear@1 368
nuclear@1 369 Matrix4x4::Matrix4x4(const Matrix3x3 &mat3x3)
nuclear@1 370 {
nuclear@1 371 reset_identity();
nuclear@1 372 for(int i=0; i<3; i++) {
nuclear@1 373 for(int j=0; j<3; j++) {
nuclear@1 374 m[i][j] = mat3x3[i][j];
nuclear@1 375 }
nuclear@1 376 }
nuclear@1 377 }
nuclear@1 378
nuclear@1 379 Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2)
nuclear@1 380 {
nuclear@1 381 Matrix4x4 res;
nuclear@1 382 const scalar_t *op1 = m1.m[0], *op2 = m2.m[0];
nuclear@1 383 scalar_t *dest = res.m[0];
nuclear@1 384
nuclear@1 385 for(int i=0; i<16; i++) {
nuclear@1 386 *dest++ = *op1++ + *op2++;
nuclear@1 387 }
nuclear@1 388 return res;
nuclear@1 389 }
nuclear@1 390
nuclear@1 391 Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2)
nuclear@1 392 {
nuclear@1 393 Matrix4x4 res;
nuclear@1 394 const scalar_t *op1 = m1.m[0], *op2 = m2.m[0];
nuclear@1 395 scalar_t *dest = res.m[0];
nuclear@1 396
nuclear@1 397 for(int i=0; i<16; i++) {
nuclear@1 398 *dest++ = *op1++ - *op2++;
nuclear@1 399 }
nuclear@1 400 return res;
nuclear@1 401 }
nuclear@1 402
nuclear@1 403 /*
nuclear@1 404 Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2) {
nuclear@1 405 Matrix4x4 res;
nuclear@1 406
nuclear@1 407 for(int i=0; i<4; i++) {
nuclear@1 408 for(int j=0; j<4; j++) {
nuclear@1 409 res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j];
nuclear@1 410 }
nuclear@1 411 }
nuclear@1 412
nuclear@1 413 return res;
nuclear@1 414 }
nuclear@1 415 */
nuclear@1 416
nuclear@1 417 void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2)
nuclear@1 418 {
nuclear@1 419 scalar_t *op1 = m1.m[0];
nuclear@1 420 const scalar_t *op2 = m2.m[0];
nuclear@1 421
nuclear@1 422 for(int i=0; i<16; i++) {
nuclear@1 423 *op1++ += *op2++;
nuclear@1 424 }
nuclear@1 425 }
nuclear@1 426
nuclear@1 427 void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2)
nuclear@1 428 {
nuclear@1 429 scalar_t *op1 = m1.m[0];
nuclear@1 430 const scalar_t *op2 = m2.m[0];
nuclear@1 431
nuclear@1 432 for(int i=0; i<16; i++) {
nuclear@1 433 *op1++ -= *op2++;
nuclear@1 434 }
nuclear@1 435 }
nuclear@1 436
nuclear@1 437 void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2)
nuclear@1 438 {
nuclear@1 439 Matrix4x4 res;
nuclear@1 440 for(int i=0; i<4; i++) {
nuclear@1 441 for(int j=0; j<4; j++) {
nuclear@1 442 res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j];
nuclear@1 443 }
nuclear@1 444 }
nuclear@1 445 memcpy(m1.m, res.m, 16 * sizeof(scalar_t));
nuclear@1 446 }
nuclear@1 447
nuclear@1 448 Matrix4x4 operator *(const Matrix4x4 &mat, scalar_t scalar)
nuclear@1 449 {
nuclear@1 450 Matrix4x4 res;
nuclear@1 451 const scalar_t *mptr = mat.m[0];
nuclear@1 452 scalar_t *dptr = res.m[0];
nuclear@1 453
nuclear@1 454 for(int i=0; i<16; i++) {
nuclear@1 455 *dptr++ = *mptr++ * scalar;
nuclear@1 456 }
nuclear@1 457 return res;
nuclear@1 458 }
nuclear@1 459
nuclear@1 460 Matrix4x4 operator *(scalar_t scalar, const Matrix4x4 &mat)
nuclear@1 461 {
nuclear@1 462 Matrix4x4 res;
nuclear@1 463 const scalar_t *mptr = mat.m[0];
nuclear@1 464 scalar_t *dptr = res.m[0];
nuclear@1 465
nuclear@1 466 for(int i=0; i<16; i++) {
nuclear@1 467 *dptr++ = *mptr++ * scalar;
nuclear@1 468 }
nuclear@1 469 return res;
nuclear@1 470 }
nuclear@1 471
nuclear@1 472 void operator *=(Matrix4x4 &mat, scalar_t scalar)
nuclear@1 473 {
nuclear@1 474 scalar_t *mptr = mat.m[0];
nuclear@1 475
nuclear@1 476 for(int i=0; i<16; i++) {
nuclear@1 477 *mptr++ *= scalar;
nuclear@1 478 }
nuclear@1 479 }
nuclear@1 480
nuclear@1 481 void Matrix4x4::translate(const Vector3 &trans)
nuclear@1 482 {
nuclear@1 483 Matrix4x4 tmat(1, 0, 0, trans.x, 0, 1, 0, trans.y, 0, 0, 1, trans.z, 0, 0, 0, 1);
nuclear@1 484 *this *= tmat;
nuclear@1 485 }
nuclear@1 486
nuclear@1 487 void Matrix4x4::set_translation(const Vector3 &trans)
nuclear@1 488 {
nuclear@1 489 *this = Matrix4x4(1, 0, 0, trans.x, 0, 1, 0, trans.y, 0, 0, 1, trans.z, 0, 0, 0, 1);
nuclear@1 490 }
nuclear@1 491
nuclear@1 492 void Matrix4x4::rotate(const Vector3 &euler_angles)
nuclear@1 493 {
nuclear@1 494 Matrix3x3 xrot, yrot, zrot;
nuclear@1 495
nuclear@1 496 xrot = Matrix3x3( 1, 0, 0,
nuclear@1 497 0, cos(euler_angles.x), -sin(euler_angles.x),
nuclear@1 498 0, sin(euler_angles.x), cos(euler_angles.x));
nuclear@1 499
nuclear@1 500 yrot = Matrix3x3( cos(euler_angles.y), 0, sin(euler_angles.y),
nuclear@1 501 0, 1, 0,
nuclear@1 502 -sin(euler_angles.y), 0, cos(euler_angles.y));
nuclear@1 503
nuclear@1 504 zrot = Matrix3x3( cos(euler_angles.z), -sin(euler_angles.z), 0,
nuclear@1 505 sin(euler_angles.z), cos(euler_angles.z), 0,
nuclear@1 506 0, 0, 1);
nuclear@1 507
nuclear@1 508 *this *= Matrix4x4(xrot * yrot * zrot);
nuclear@1 509 }
nuclear@1 510
nuclear@1 511 void Matrix4x4::set_rotation(const Vector3 &euler_angles)
nuclear@1 512 {
nuclear@1 513 Matrix3x3 xrot, yrot, zrot;
nuclear@1 514
nuclear@1 515 xrot = Matrix3x3( 1, 0, 0,
nuclear@1 516 0, cos(euler_angles.x), -sin(euler_angles.x),
nuclear@1 517 0, sin(euler_angles.x), cos(euler_angles.x));
nuclear@1 518
nuclear@1 519 yrot = Matrix3x3( cos(euler_angles.y), 0, sin(euler_angles.y),
nuclear@1 520 0, 1, 0,
nuclear@1 521 -sin(euler_angles.y), 0, cos(euler_angles.y));
nuclear@1 522
nuclear@1 523 zrot = Matrix3x3( cos(euler_angles.z), -sin(euler_angles.z), 0,
nuclear@1 524 sin(euler_angles.z), cos(euler_angles.z), 0,
nuclear@1 525 0, 0, 1);
nuclear@1 526
nuclear@1 527 *this = Matrix4x4(xrot * yrot * zrot);
nuclear@1 528 }
nuclear@1 529
nuclear@1 530 void Matrix4x4::rotate(const Vector3 &axis, scalar_t angle)
nuclear@1 531 {
nuclear@1 532 scalar_t sina = (scalar_t)sin(angle);
nuclear@1 533 scalar_t cosa = (scalar_t)cos(angle);
nuclear@1 534 scalar_t invcosa = 1-cosa;
nuclear@1 535 scalar_t nxsq = axis.x * axis.x;
nuclear@1 536 scalar_t nysq = axis.y * axis.y;
nuclear@1 537 scalar_t nzsq = axis.z * axis.z;
nuclear@1 538
nuclear@1 539 Matrix3x3 xform;
nuclear@1 540 xform[0][0] = nxsq + (1-nxsq) * cosa;
nuclear@1 541 xform[0][1] = axis.x * axis.y * invcosa - axis.z * sina;
nuclear@1 542 xform[0][2] = axis.x * axis.z * invcosa + axis.y * sina;
nuclear@1 543 xform[1][0] = axis.x * axis.y * invcosa + axis.z * sina;
nuclear@1 544 xform[1][1] = nysq + (1-nysq) * cosa;
nuclear@1 545 xform[1][2] = axis.y * axis.z * invcosa - axis.x * sina;
nuclear@1 546 xform[2][0] = axis.x * axis.z * invcosa - axis.y * sina;
nuclear@1 547 xform[2][1] = axis.y * axis.z * invcosa + axis.x * sina;
nuclear@1 548 xform[2][2] = nzsq + (1-nzsq) * cosa;
nuclear@1 549
nuclear@1 550 *this *= Matrix4x4(xform);
nuclear@1 551 }
nuclear@1 552
nuclear@1 553 void Matrix4x4::set_rotation(const Vector3 &axis, scalar_t angle)
nuclear@1 554 {
nuclear@1 555 scalar_t sina = (scalar_t)sin(angle);
nuclear@1 556 scalar_t cosa = (scalar_t)cos(angle);
nuclear@1 557 scalar_t invcosa = 1-cosa;
nuclear@1 558 scalar_t nxsq = axis.x * axis.x;
nuclear@1 559 scalar_t nysq = axis.y * axis.y;
nuclear@1 560 scalar_t nzsq = axis.z * axis.z;
nuclear@1 561
nuclear@1 562 reset_identity();
nuclear@1 563 m[0][0] = nxsq + (1-nxsq) * cosa;
nuclear@1 564 m[0][1] = axis.x * axis.y * invcosa - axis.z * sina;
nuclear@1 565 m[0][2] = axis.x * axis.z * invcosa + axis.y * sina;
nuclear@1 566 m[1][0] = axis.x * axis.y * invcosa + axis.z * sina;
nuclear@1 567 m[1][1] = nysq + (1-nysq) * cosa;
nuclear@1 568 m[1][2] = axis.y * axis.z * invcosa - axis.x * sina;
nuclear@1 569 m[2][0] = axis.x * axis.z * invcosa - axis.y * sina;
nuclear@1 570 m[2][1] = axis.y * axis.z * invcosa + axis.x * sina;
nuclear@1 571 m[2][2] = nzsq + (1-nzsq) * cosa;
nuclear@1 572 }
nuclear@1 573
nuclear@1 574 void Matrix4x4::scale(const Vector4 &scale_vec)
nuclear@1 575 {
nuclear@1 576 Matrix4x4 smat( scale_vec.x, 0, 0, 0,
nuclear@1 577 0, scale_vec.y, 0, 0,
nuclear@1 578 0, 0, scale_vec.z, 0,
nuclear@1 579 0, 0, 0, scale_vec.w);
nuclear@1 580 *this *= smat;
nuclear@1 581 }
nuclear@1 582
nuclear@1 583 void Matrix4x4::set_scaling(const Vector4 &scale_vec)
nuclear@1 584 {
nuclear@1 585 *this = Matrix4x4( scale_vec.x, 0, 0, 0,
nuclear@1 586 0, scale_vec.y, 0, 0,
nuclear@1 587 0, 0, scale_vec.z, 0,
nuclear@1 588 0, 0, 0, scale_vec.w);
nuclear@1 589 }
nuclear@1 590
nuclear@1 591 void Matrix4x4::set_column_vector(const Vector4 &vec, unsigned int col_index)
nuclear@1 592 {
nuclear@1 593 m[0][col_index] = vec.x;
nuclear@1 594 m[1][col_index] = vec.y;
nuclear@1 595 m[2][col_index] = vec.z;
nuclear@1 596 m[3][col_index] = vec.w;
nuclear@1 597 }
nuclear@1 598
nuclear@1 599 void Matrix4x4::set_row_vector(const Vector4 &vec, unsigned int row_index)
nuclear@1 600 {
nuclear@1 601 m[row_index][0] = vec.x;
nuclear@1 602 m[row_index][1] = vec.y;
nuclear@1 603 m[row_index][2] = vec.z;
nuclear@1 604 m[row_index][3] = vec.w;
nuclear@1 605 }
nuclear@1 606
nuclear@1 607 Vector4 Matrix4x4::get_column_vector(unsigned int col_index) const
nuclear@1 608 {
nuclear@1 609 return Vector4(m[0][col_index], m[1][col_index], m[2][col_index], m[3][col_index]);
nuclear@1 610 }
nuclear@1 611
nuclear@1 612 Vector4 Matrix4x4::get_row_vector(unsigned int row_index) const
nuclear@1 613 {
nuclear@1 614 return Vector4(m[row_index][0], m[row_index][1], m[row_index][2], m[row_index][3]);
nuclear@1 615 }
nuclear@1 616
nuclear@1 617 void Matrix4x4::transpose()
nuclear@1 618 {
nuclear@1 619 Matrix4x4 tmp = *this;
nuclear@1 620 for(int i=0; i<4; i++) {
nuclear@1 621 for(int j=0; j<4; j++) {
nuclear@1 622 m[i][j] = tmp[j][i];
nuclear@1 623 }
nuclear@1 624 }
nuclear@1 625 }
nuclear@1 626
nuclear@1 627 Matrix4x4 Matrix4x4::transposed() const
nuclear@1 628 {
nuclear@1 629 Matrix4x4 res;
nuclear@1 630 for(int i=0; i<4; i++) {
nuclear@1 631 for(int j=0; j<4; j++) {
nuclear@1 632 res[i][j] = m[j][i];
nuclear@1 633 }
nuclear@1 634 }
nuclear@1 635 return res;
nuclear@1 636 }
nuclear@1 637
nuclear@1 638 scalar_t Matrix4x4::determinant() const
nuclear@1 639 {
nuclear@1 640 scalar_t det11 = (m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
nuclear@1 641 (m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
nuclear@1 642 (m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
nuclear@1 643
nuclear@1 644 scalar_t det12 = (m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
nuclear@1 645 (m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
nuclear@1 646 (m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
nuclear@1 647
nuclear@1 648 scalar_t det13 = (m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
nuclear@1 649 (m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
nuclear@1 650 (m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
nuclear@1 651
nuclear@1 652 scalar_t det14 = (m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
nuclear@1 653 (m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
nuclear@1 654 (m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
nuclear@1 655
nuclear@1 656 return m[0][0] * det11 - m[0][1] * det12 + m[0][2] * det13 - m[0][3] * det14;
nuclear@1 657 }
nuclear@1 658
nuclear@1 659
nuclear@1 660 Matrix4x4 Matrix4x4::adjoint() const
nuclear@1 661 {
nuclear@1 662 Matrix4x4 coef;
nuclear@1 663
nuclear@1 664 coef.m[0][0] = (m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
nuclear@1 665 (m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
nuclear@1 666 (m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
nuclear@1 667 coef.m[0][1] = (m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
nuclear@1 668 (m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
nuclear@1 669 (m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
nuclear@1 670 coef.m[0][2] = (m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
nuclear@1 671 (m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
nuclear@1 672 (m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
nuclear@1 673 coef.m[0][3] = (m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
nuclear@1 674 (m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
nuclear@1 675 (m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
nuclear@1 676
nuclear@1 677 coef.m[1][0] = (m[0][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
nuclear@1 678 (m[0][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
nuclear@1 679 (m[0][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
nuclear@1 680 coef.m[1][1] = (m[0][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
nuclear@1 681 (m[0][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
nuclear@1 682 (m[0][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
nuclear@1 683 coef.m[1][2] = (m[0][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
nuclear@1 684 (m[0][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
nuclear@1 685 (m[0][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
nuclear@1 686 coef.m[1][3] = (m[0][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
nuclear@1 687 (m[0][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
nuclear@1 688 (m[0][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
nuclear@1 689
nuclear@1 690 coef.m[2][0] = (m[0][1] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) -
nuclear@1 691 (m[0][2] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) +
nuclear@1 692 (m[0][3] * (m[1][1] * m[3][2] - m[3][1] * m[1][2]));
nuclear@1 693 coef.m[2][1] = (m[0][0] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) -
nuclear@1 694 (m[0][2] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) +
nuclear@1 695 (m[0][3] * (m[1][0] * m[3][2] - m[3][0] * m[1][2]));
nuclear@1 696 coef.m[2][2] = (m[0][0] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) -
nuclear@1 697 (m[0][1] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) +
nuclear@1 698 (m[0][3] * (m[1][0] * m[3][1] - m[3][0] * m[1][1]));
nuclear@1 699 coef.m[2][3] = (m[0][0] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])) -
nuclear@1 700 (m[0][1] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])) +
nuclear@1 701 (m[0][2] * (m[1][0] * m[3][1] - m[3][0] * m[1][1]));
nuclear@1 702
nuclear@1 703 coef.m[3][0] = (m[0][1] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) -
nuclear@1 704 (m[0][2] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) +
nuclear@1 705 (m[0][3] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]));
nuclear@1 706 coef.m[3][1] = (m[0][0] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) -
nuclear@1 707 (m[0][2] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) +
nuclear@1 708 (m[0][3] * (m[1][0] * m[2][2] - m[2][0] * m[1][2]));
nuclear@1 709 coef.m[3][2] = (m[0][0] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) -
nuclear@1 710 (m[0][1] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) +
nuclear@1 711 (m[0][3] * (m[1][0] * m[2][1] - m[2][0] * m[1][1]));
nuclear@1 712 coef.m[3][3] = (m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])) -
nuclear@1 713 (m[0][1] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])) +
nuclear@1 714 (m[0][2] * (m[1][0] * m[2][1] - m[2][0] * m[1][1]));
nuclear@1 715
nuclear@1 716 coef.transpose();
nuclear@1 717
nuclear@1 718 for(int i=0; i<4; i++) {
nuclear@1 719 for(int j=0; j<4; j++) {
nuclear@1 720 coef.m[i][j] = j%2 ? -coef.m[i][j] : coef.m[i][j];
nuclear@1 721 if(i%2) coef.m[i][j] = -coef.m[i][j];
nuclear@1 722 }
nuclear@1 723 }
nuclear@1 724
nuclear@1 725 return coef;
nuclear@1 726 }
nuclear@1 727
nuclear@1 728 Matrix4x4 Matrix4x4::inverse() const
nuclear@1 729 {
nuclear@1 730 Matrix4x4 adj = adjoint();
nuclear@1 731
nuclear@1 732 return adj * (1.0f / determinant());
nuclear@1 733 }
nuclear@1 734
nuclear@1 735 ostream &operator <<(ostream &out, const Matrix4x4 &mat)
nuclear@1 736 {
nuclear@1 737 for(int i=0; i<4; i++) {
nuclear@1 738 char str[100];
nuclear@1 739 sprintf(str, "[ %12.5f %12.5f %12.5f %12.5f ]\n", (float)mat.m[i][0], (float)mat.m[i][1], (float)mat.m[i][2], (float)mat.m[i][3]);
nuclear@1 740 out << str;
nuclear@1 741 }
nuclear@1 742 return out;
nuclear@1 743 }