Courtney Goeltzenleuchter | 4825f6a | 2014-10-28 10:27:47 -0600 | [diff] [blame^] | 1 | /* |
| 2 | DO WHAT THE **** YOU WANT TO PUBLIC LICENSE |
| 3 | Version 2, December 2004 |
| 4 | |
| 5 | Copyright (C) 2013 Wolfgang 'datenwolf' Draxinger <code@datenwolf.net> |
| 6 | |
| 7 | Everyone is permitted to copy and distribute verbatim or modified |
| 8 | copies of this license document, and changing it is allowed as long |
| 9 | as the name is changed. |
| 10 | |
| 11 | DO WHAT THE **** YOU WANT TO PUBLIC LICENSE |
| 12 | TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION |
| 13 | |
| 14 | 0. You just DO WHAT THE FUCK YOU WANT TO. |
| 15 | */ |
| 16 | |
| 17 | #ifndef LINMATH_H |
| 18 | #define LINMATH_H |
| 19 | |
| 20 | #define __USE_BSD |
| 21 | #include <math.h> |
| 22 | |
| 23 | // Converts degrees to radians. |
| 24 | #define degreesToRadians(angleDegrees) (angleDegrees * M_PI / 180.0) |
| 25 | |
| 26 | // Converts radians to degrees. |
| 27 | #define radiansToDegrees(angleRadians) (angleRadians * 180.0 / M_PI) |
| 28 | |
| 29 | typedef float vec3[3]; |
| 30 | static inline void vec3_add(vec3 r, vec3 const a, vec3 const b) |
| 31 | { |
| 32 | int i; |
| 33 | for(i=0; i<3; ++i) |
| 34 | r[i] = a[i] + b[i]; |
| 35 | } |
| 36 | static inline void vec3_sub(vec3 r, vec3 const a, vec3 const b) |
| 37 | { |
| 38 | int i; |
| 39 | for(i=0; i<3; ++i) |
| 40 | r[i] = a[i] - b[i]; |
| 41 | } |
| 42 | static inline void vec3_scale(vec3 r, vec3 const v, float const s) |
| 43 | { |
| 44 | int i; |
| 45 | for(i=0; i<3; ++i) |
| 46 | r[i] = v[i] * s; |
| 47 | } |
| 48 | static inline float vec3_mul_inner(vec3 const a, vec3 const b) |
| 49 | { |
| 50 | float p = 0.f; |
| 51 | int i; |
| 52 | for(i=0; i<3; ++i) |
| 53 | p += b[i]*a[i]; |
| 54 | return p; |
| 55 | } |
| 56 | static inline void vec3_mul_cross(vec3 r, vec3 const a, vec3 const b) |
| 57 | { |
| 58 | r[0] = a[1]*b[2] - a[2]*b[1]; |
| 59 | r[1] = a[2]*b[0] - a[0]*b[2]; |
| 60 | r[2] = a[0]*b[1] - a[1]*b[0]; |
| 61 | } |
| 62 | static inline float vec3_len(vec3 const v) |
| 63 | { |
| 64 | return sqrtf(vec3_mul_inner(v, v)); |
| 65 | } |
| 66 | static inline void vec3_norm(vec3 r, vec3 const v) |
| 67 | { |
| 68 | float k = 1.f / vec3_len(v); |
| 69 | vec3_scale(r, v, k); |
| 70 | } |
| 71 | static inline void vec3_reflect(vec3 r, vec3 const v, vec3 const n) |
| 72 | { |
| 73 | float p = 2.f*vec3_mul_inner(v, n); |
| 74 | int i; |
| 75 | for(i=0;i<3;++i) |
| 76 | r[i] = v[i] - p*n[i]; |
| 77 | } |
| 78 | |
| 79 | typedef float vec4[4]; |
| 80 | static inline void vec4_add(vec4 r, vec4 const a, vec4 const b) |
| 81 | { |
| 82 | int i; |
| 83 | for(i=0; i<4; ++i) |
| 84 | r[i] = a[i] + b[i]; |
| 85 | } |
| 86 | static inline void vec4_sub(vec4 r, vec4 const a, vec4 const b) |
| 87 | { |
| 88 | int i; |
| 89 | for(i=0; i<4; ++i) |
| 90 | r[i] = a[i] - b[i]; |
| 91 | } |
| 92 | static inline void vec4_scale(vec4 r, vec4 v, float s) |
| 93 | { |
| 94 | int i; |
| 95 | for(i=0; i<4; ++i) |
| 96 | r[i] = v[i] * s; |
| 97 | } |
| 98 | static inline float vec4_mul_inner(vec4 a, vec4 b) |
| 99 | { |
| 100 | float p = 0.f; |
| 101 | int i; |
| 102 | for(i=0; i<4; ++i) |
| 103 | p += b[i]*a[i]; |
| 104 | return p; |
| 105 | } |
| 106 | static inline void vec4_mul_cross(vec4 r, vec4 a, vec4 b) |
| 107 | { |
| 108 | r[0] = a[1]*b[2] - a[2]*b[1]; |
| 109 | r[1] = a[2]*b[0] - a[0]*b[2]; |
| 110 | r[2] = a[0]*b[1] - a[1]*b[0]; |
| 111 | r[3] = 1.f; |
| 112 | } |
| 113 | static inline float vec4_len(vec4 v) |
| 114 | { |
| 115 | return sqrtf(vec4_mul_inner(v, v)); |
| 116 | } |
| 117 | static inline void vec4_norm(vec4 r, vec4 v) |
| 118 | { |
| 119 | float k = 1.f / vec4_len(v); |
| 120 | vec4_scale(r, v, k); |
| 121 | } |
| 122 | static inline void vec4_reflect(vec4 r, vec4 v, vec4 n) |
| 123 | { |
| 124 | float p = 2.f*vec4_mul_inner(v, n); |
| 125 | int i; |
| 126 | for(i=0;i<4;++i) |
| 127 | r[i] = v[i] - p*n[i]; |
| 128 | } |
| 129 | |
| 130 | typedef vec4 mat4x4[4]; |
| 131 | static inline void mat4x4_identity(mat4x4 M) |
| 132 | { |
| 133 | int i, j; |
| 134 | for(i=0; i<4; ++i) |
| 135 | for(j=0; j<4; ++j) |
| 136 | M[i][j] = i==j ? 1.f : 0.f; |
| 137 | } |
| 138 | static inline void mat4x4_dup(mat4x4 M, mat4x4 N) |
| 139 | { |
| 140 | int i, j; |
| 141 | for(i=0; i<4; ++i) |
| 142 | for(j=0; j<4; ++j) |
| 143 | M[i][j] = N[i][j]; |
| 144 | } |
| 145 | static inline void mat4x4_row(vec4 r, mat4x4 M, int i) |
| 146 | { |
| 147 | int k; |
| 148 | for(k=0; k<4; ++k) |
| 149 | r[k] = M[k][i]; |
| 150 | } |
| 151 | static inline void mat4x4_col(vec4 r, mat4x4 M, int i) |
| 152 | { |
| 153 | int k; |
| 154 | for(k=0; k<4; ++k) |
| 155 | r[k] = M[i][k]; |
| 156 | } |
| 157 | static inline void mat4x4_transpose(mat4x4 M, mat4x4 N) |
| 158 | { |
| 159 | int i, j; |
| 160 | for(j=0; j<4; ++j) |
| 161 | for(i=0; i<4; ++i) |
| 162 | M[i][j] = N[j][i]; |
| 163 | } |
| 164 | static inline void mat4x4_add(mat4x4 M, mat4x4 a, mat4x4 b) |
| 165 | { |
| 166 | int i; |
| 167 | for(i=0; i<4; ++i) |
| 168 | vec4_add(M[i], a[i], b[i]); |
| 169 | } |
| 170 | static inline void mat4x4_sub(mat4x4 M, mat4x4 a, mat4x4 b) |
| 171 | { |
| 172 | int i; |
| 173 | for(i=0; i<4; ++i) |
| 174 | vec4_sub(M[i], a[i], b[i]); |
| 175 | } |
| 176 | static inline void mat4x4_scale(mat4x4 M, mat4x4 a, float k) |
| 177 | { |
| 178 | int i; |
| 179 | for(i=0; i<4; ++i) |
| 180 | vec4_scale(M[i], a[i], k); |
| 181 | } |
| 182 | static inline void mat4x4_scale_aniso(mat4x4 M, mat4x4 a, float x, float y, float z) |
| 183 | { |
| 184 | int i; |
| 185 | vec4_scale(M[0], a[0], x); |
| 186 | vec4_scale(M[1], a[1], y); |
| 187 | vec4_scale(M[2], a[2], z); |
| 188 | for(i = 0; i < 4; ++i) { |
| 189 | M[3][i] = a[3][i]; |
| 190 | } |
| 191 | } |
| 192 | static inline void mat4x4_mul(mat4x4 M, mat4x4 a, mat4x4 b) |
| 193 | { |
| 194 | int k, r, c; |
| 195 | for(c=0; c<4; ++c) for(r=0; r<4; ++r) { |
| 196 | M[c][r] = 0.f; |
| 197 | for(k=0; k<4; ++k) |
| 198 | M[c][r] += a[k][r] * b[c][k]; |
| 199 | } |
| 200 | } |
| 201 | static inline void mat4x4_mul_vec4(vec4 r, mat4x4 M, vec4 v) |
| 202 | { |
| 203 | int i, j; |
| 204 | for(j=0; j<4; ++j) { |
| 205 | r[j] = 0.f; |
| 206 | for(i=0; i<4; ++i) |
| 207 | r[j] += M[i][j] * v[i]; |
| 208 | } |
| 209 | } |
| 210 | static inline void mat4x4_translate(mat4x4 T, float x, float y, float z) |
| 211 | { |
| 212 | mat4x4_identity(T); |
| 213 | T[3][0] = x; |
| 214 | T[3][1] = y; |
| 215 | T[3][2] = z; |
| 216 | } |
| 217 | static inline void mat4x4_translate_in_place(mat4x4 M, float x, float y, float z) |
| 218 | { |
| 219 | vec4 t = {x, y, z, 0}; |
| 220 | vec4 r; |
| 221 | int i; |
| 222 | for (i = 0; i < 4; ++i) { |
| 223 | mat4x4_row(r, M, i); |
| 224 | M[3][i] += vec4_mul_inner(r, t); |
| 225 | } |
| 226 | } |
| 227 | static inline void mat4x4_from_vec3_mul_outer(mat4x4 M, vec3 a, vec3 b) |
| 228 | { |
| 229 | int i, j; |
| 230 | for(i=0; i<4; ++i) for(j=0; j<4; ++j) |
| 231 | M[i][j] = i<3 && j<3 ? a[i] * b[j] : 0.f; |
| 232 | } |
| 233 | static inline void mat4x4_rotate(mat4x4 R, mat4x4 M, float x, float y, float z, float angle) |
| 234 | { |
| 235 | float s = sinf(angle); |
| 236 | float c = cosf(angle); |
| 237 | vec3 u = {x, y, z}; |
| 238 | |
| 239 | if(vec3_len(u) > 1e-4) { |
| 240 | vec3_norm(u, u); |
| 241 | mat4x4 T; |
| 242 | mat4x4_from_vec3_mul_outer(T, u, u); |
| 243 | |
| 244 | mat4x4 S = { |
| 245 | { 0, u[2], -u[1], 0}, |
| 246 | {-u[2], 0, u[0], 0}, |
| 247 | { u[1], -u[0], 0, 0}, |
| 248 | { 0, 0, 0, 0} |
| 249 | }; |
| 250 | mat4x4_scale(S, S, s); |
| 251 | |
| 252 | mat4x4 C; |
| 253 | mat4x4_identity(C); |
| 254 | mat4x4_sub(C, C, T); |
| 255 | |
| 256 | mat4x4_scale(C, C, c); |
| 257 | |
| 258 | mat4x4_add(T, T, C); |
| 259 | mat4x4_add(T, T, S); |
| 260 | |
| 261 | T[3][3] = 1.; |
| 262 | mat4x4_mul(R, M, T); |
| 263 | } else { |
| 264 | mat4x4_dup(R, M); |
| 265 | } |
| 266 | } |
| 267 | static inline void mat4x4_rotate_X(mat4x4 Q, mat4x4 M, float angle) |
| 268 | { |
| 269 | float s = sinf(angle); |
| 270 | float c = cosf(angle); |
| 271 | mat4x4 R = { |
| 272 | {1.f, 0.f, 0.f, 0.f}, |
| 273 | {0.f, c, s, 0.f}, |
| 274 | {0.f, -s, c, 0.f}, |
| 275 | {0.f, 0.f, 0.f, 1.f} |
| 276 | }; |
| 277 | mat4x4_mul(Q, M, R); |
| 278 | } |
| 279 | static inline void mat4x4_rotate_Y(mat4x4 Q, mat4x4 M, float angle) |
| 280 | { |
| 281 | float s = sinf(angle); |
| 282 | float c = cosf(angle); |
| 283 | mat4x4 R = { |
| 284 | { c, 0.f, s, 0.f}, |
| 285 | { 0.f, 1.f, 0.f, 0.f}, |
| 286 | { -s, 0.f, c, 0.f}, |
| 287 | { 0.f, 0.f, 0.f, 1.f} |
| 288 | }; |
| 289 | mat4x4_mul(Q, M, R); |
| 290 | } |
| 291 | static inline void mat4x4_rotate_Z(mat4x4 Q, mat4x4 M, float angle) |
| 292 | { |
| 293 | float s = sinf(angle); |
| 294 | float c = cosf(angle); |
| 295 | mat4x4 R = { |
| 296 | { c, s, 0.f, 0.f}, |
| 297 | { -s, c, 0.f, 0.f}, |
| 298 | { 0.f, 0.f, 1.f, 0.f}, |
| 299 | { 0.f, 0.f, 0.f, 1.f} |
| 300 | }; |
| 301 | mat4x4_mul(Q, M, R); |
| 302 | } |
| 303 | static inline void mat4x4_invert(mat4x4 T, mat4x4 M) |
| 304 | { |
| 305 | float s[6]; |
| 306 | float c[6]; |
| 307 | s[0] = M[0][0]*M[1][1] - M[1][0]*M[0][1]; |
| 308 | s[1] = M[0][0]*M[1][2] - M[1][0]*M[0][2]; |
| 309 | s[2] = M[0][0]*M[1][3] - M[1][0]*M[0][3]; |
| 310 | s[3] = M[0][1]*M[1][2] - M[1][1]*M[0][2]; |
| 311 | s[4] = M[0][1]*M[1][3] - M[1][1]*M[0][3]; |
| 312 | s[5] = M[0][2]*M[1][3] - M[1][2]*M[0][3]; |
| 313 | |
| 314 | c[0] = M[2][0]*M[3][1] - M[3][0]*M[2][1]; |
| 315 | c[1] = M[2][0]*M[3][2] - M[3][0]*M[2][2]; |
| 316 | c[2] = M[2][0]*M[3][3] - M[3][0]*M[2][3]; |
| 317 | c[3] = M[2][1]*M[3][2] - M[3][1]*M[2][2]; |
| 318 | c[4] = M[2][1]*M[3][3] - M[3][1]*M[2][3]; |
| 319 | c[5] = M[2][2]*M[3][3] - M[3][2]*M[2][3]; |
| 320 | |
| 321 | /* Assumes it is invertible */ |
| 322 | float idet = 1.0f/( s[0]*c[5]-s[1]*c[4]+s[2]*c[3]+s[3]*c[2]-s[4]*c[1]+s[5]*c[0] ); |
| 323 | |
| 324 | T[0][0] = ( M[1][1] * c[5] - M[1][2] * c[4] + M[1][3] * c[3]) * idet; |
| 325 | T[0][1] = (-M[0][1] * c[5] + M[0][2] * c[4] - M[0][3] * c[3]) * idet; |
| 326 | T[0][2] = ( M[3][1] * s[5] - M[3][2] * s[4] + M[3][3] * s[3]) * idet; |
| 327 | T[0][3] = (-M[2][1] * s[5] + M[2][2] * s[4] - M[2][3] * s[3]) * idet; |
| 328 | |
| 329 | T[1][0] = (-M[1][0] * c[5] + M[1][2] * c[2] - M[1][3] * c[1]) * idet; |
| 330 | T[1][1] = ( M[0][0] * c[5] - M[0][2] * c[2] + M[0][3] * c[1]) * idet; |
| 331 | T[1][2] = (-M[3][0] * s[5] + M[3][2] * s[2] - M[3][3] * s[1]) * idet; |
| 332 | T[1][3] = ( M[2][0] * s[5] - M[2][2] * s[2] + M[2][3] * s[1]) * idet; |
| 333 | |
| 334 | T[2][0] = ( M[1][0] * c[4] - M[1][1] * c[2] + M[1][3] * c[0]) * idet; |
| 335 | T[2][1] = (-M[0][0] * c[4] + M[0][1] * c[2] - M[0][3] * c[0]) * idet; |
| 336 | T[2][2] = ( M[3][0] * s[4] - M[3][1] * s[2] + M[3][3] * s[0]) * idet; |
| 337 | T[2][3] = (-M[2][0] * s[4] + M[2][1] * s[2] - M[2][3] * s[0]) * idet; |
| 338 | |
| 339 | T[3][0] = (-M[1][0] * c[3] + M[1][1] * c[1] - M[1][2] * c[0]) * idet; |
| 340 | T[3][1] = ( M[0][0] * c[3] - M[0][1] * c[1] + M[0][2] * c[0]) * idet; |
| 341 | T[3][2] = (-M[3][0] * s[3] + M[3][1] * s[1] - M[3][2] * s[0]) * idet; |
| 342 | T[3][3] = ( M[2][0] * s[3] - M[2][1] * s[1] + M[2][2] * s[0]) * idet; |
| 343 | } |
| 344 | static inline void mat4x4_orthonormalize(mat4x4 R, mat4x4 M) |
| 345 | { |
| 346 | mat4x4_dup(R, M); |
| 347 | float s = 1.; |
| 348 | vec3 h; |
| 349 | |
| 350 | vec3_norm(R[2], R[2]); |
| 351 | |
| 352 | s = vec3_mul_inner(R[1], R[2]); |
| 353 | vec3_scale(h, R[2], s); |
| 354 | vec3_sub(R[1], R[1], h); |
| 355 | vec3_norm(R[2], R[2]); |
| 356 | |
| 357 | s = vec3_mul_inner(R[1], R[2]); |
| 358 | vec3_scale(h, R[2], s); |
| 359 | vec3_sub(R[1], R[1], h); |
| 360 | vec3_norm(R[1], R[1]); |
| 361 | |
| 362 | s = vec3_mul_inner(R[0], R[1]); |
| 363 | vec3_scale(h, R[1], s); |
| 364 | vec3_sub(R[0], R[0], h); |
| 365 | vec3_norm(R[0], R[0]); |
| 366 | } |
| 367 | |
| 368 | static inline void mat4x4_frustum(mat4x4 M, float l, float r, float b, float t, float n, float f) |
| 369 | { |
| 370 | M[0][0] = 2.f*n/(r-l); |
| 371 | M[0][1] = M[0][2] = M[0][3] = 0.f; |
| 372 | |
| 373 | M[1][1] = 2.*n/(t-b); |
| 374 | M[1][0] = M[1][2] = M[1][3] = 0.f; |
| 375 | |
| 376 | M[2][0] = (r+l)/(r-l); |
| 377 | M[2][1] = (t+b)/(t-b); |
| 378 | M[2][2] = -(f+n)/(f-n); |
| 379 | M[2][3] = -1.f; |
| 380 | |
| 381 | M[3][2] = -2.f*(f*n)/(f-n); |
| 382 | M[3][0] = M[3][1] = M[3][3] = 0.f; |
| 383 | } |
| 384 | static inline void mat4x4_ortho(mat4x4 M, float l, float r, float b, float t, float n, float f) |
| 385 | { |
| 386 | M[0][0] = 2.f/(r-l); |
| 387 | M[0][1] = M[0][2] = M[0][3] = 0.f; |
| 388 | |
| 389 | M[1][1] = 2.f/(t-b); |
| 390 | M[1][0] = M[1][2] = M[1][3] = 0.f; |
| 391 | |
| 392 | M[2][2] = -2.f/(f-n); |
| 393 | M[2][0] = M[2][1] = M[2][3] = 0.f; |
| 394 | |
| 395 | M[3][0] = -(r+l)/(r-l); |
| 396 | M[3][1] = -(t+b)/(t-b); |
| 397 | M[3][2] = -(f+n)/(f-n); |
| 398 | M[3][3] = 1.f; |
| 399 | } |
| 400 | static inline void mat4x4_perspective(mat4x4 m, float y_fov, float aspect, float n, float f) |
| 401 | { |
| 402 | /* NOTE: Degrees are an unhandy unit to work with. |
| 403 | * linmath.h uses radians for everything! */ |
| 404 | float const a = 1.f / tan(y_fov / 2.f); |
| 405 | |
| 406 | m[0][0] = a / aspect; |
| 407 | m[0][1] = 0.f; |
| 408 | m[0][2] = 0.f; |
| 409 | m[0][3] = 0.f; |
| 410 | |
| 411 | m[1][0] = 0.f; |
| 412 | m[1][1] = a; |
| 413 | m[1][2] = 0.f; |
| 414 | m[1][3] = 0.f; |
| 415 | |
| 416 | m[2][0] = 0.f; |
| 417 | m[2][1] = 0.f; |
| 418 | m[2][2] = -((f + n) / (f - n)); |
| 419 | m[2][3] = -1.f; |
| 420 | |
| 421 | m[3][0] = 0.f; |
| 422 | m[3][1] = 0.f; |
| 423 | m[3][2] = -((2.f * f * n) / (f - n)); |
| 424 | m[3][3] = 0.f; |
| 425 | } |
| 426 | static inline void mat4x4_look_at(mat4x4 m, vec3 eye, vec3 center, vec3 up) |
| 427 | { |
| 428 | /* Adapted from Android's OpenGL Matrix.java. */ |
| 429 | /* See the OpenGL GLUT documentation for gluLookAt for a description */ |
| 430 | /* of the algorithm. We implement it in a straightforward way: */ |
| 431 | |
| 432 | /* TODO: The negation of of can be spared by swapping the order of |
| 433 | * operands in the following cross products in the right way. */ |
| 434 | vec3 f; |
| 435 | vec3_sub(f, center, eye); |
| 436 | vec3_norm(f, f); |
| 437 | |
| 438 | vec3 s; |
| 439 | vec3_mul_cross(s, f, up); |
| 440 | vec3_norm(s, s); |
| 441 | |
| 442 | vec3 t; |
| 443 | vec3_mul_cross(t, s, f); |
| 444 | |
| 445 | m[0][0] = s[0]; |
| 446 | m[0][1] = t[0]; |
| 447 | m[0][2] = -f[0]; |
| 448 | m[0][3] = 0.f; |
| 449 | |
| 450 | m[1][0] = s[1]; |
| 451 | m[1][1] = t[1]; |
| 452 | m[1][2] = -f[1]; |
| 453 | m[1][3] = 0.f; |
| 454 | |
| 455 | m[2][0] = s[2]; |
| 456 | m[2][1] = t[2]; |
| 457 | m[2][2] = -f[2]; |
| 458 | m[2][3] = 0.f; |
| 459 | |
| 460 | m[3][0] = 0.f; |
| 461 | m[3][1] = 0.f; |
| 462 | m[3][2] = 0.f; |
| 463 | m[3][3] = 1.f; |
| 464 | |
| 465 | mat4x4_translate_in_place(m, -eye[0], -eye[1], -eye[2]); |
| 466 | } |
| 467 | |
| 468 | typedef float quat[4]; |
| 469 | static inline void quat_identity(quat q) |
| 470 | { |
| 471 | q[0] = q[1] = q[2] = 0.f; |
| 472 | q[3] = 1.f; |
| 473 | } |
| 474 | static inline void quat_add(quat r, quat a, quat b) |
| 475 | { |
| 476 | int i; |
| 477 | for(i=0; i<4; ++i) |
| 478 | r[i] = a[i] + b[i]; |
| 479 | } |
| 480 | static inline void quat_sub(quat r, quat a, quat b) |
| 481 | { |
| 482 | int i; |
| 483 | for(i=0; i<4; ++i) |
| 484 | r[i] = a[i] - b[i]; |
| 485 | } |
| 486 | static inline void quat_mul(quat r, quat p, quat q) |
| 487 | { |
| 488 | vec3 w; |
| 489 | vec3_mul_cross(r, p, q); |
| 490 | vec3_scale(w, p, q[3]); |
| 491 | vec3_add(r, r, w); |
| 492 | vec3_scale(w, q, p[3]); |
| 493 | vec3_add(r, r, w); |
| 494 | r[3] = p[3]*q[3] - vec3_mul_inner(p, q); |
| 495 | } |
| 496 | static inline void quat_scale(quat r, quat v, float s) |
| 497 | { |
| 498 | int i; |
| 499 | for(i=0; i<4; ++i) |
| 500 | r[i] = v[i] * s; |
| 501 | } |
| 502 | static inline float quat_inner_product(quat a, quat b) |
| 503 | { |
| 504 | float p = 0.f; |
| 505 | int i; |
| 506 | for(i=0; i<4; ++i) |
| 507 | p += b[i]*a[i]; |
| 508 | return p; |
| 509 | } |
| 510 | static inline void quat_conj(quat r, quat q) |
| 511 | { |
| 512 | int i; |
| 513 | for(i=0; i<3; ++i) |
| 514 | r[i] = -q[i]; |
| 515 | r[3] = q[3]; |
| 516 | } |
| 517 | #define quat_norm vec4_norm |
| 518 | static inline void quat_mul_vec3(vec3 r, quat q, vec3 v) |
| 519 | { |
| 520 | quat v_ = {v[0], v[1], v[2], 0.f}; |
| 521 | |
| 522 | quat_conj(r, q); |
| 523 | quat_norm(r, r); |
| 524 | quat_mul(r, v_, r); |
| 525 | quat_mul(r, q, r); |
| 526 | } |
| 527 | static inline void mat4x4_from_quat(mat4x4 M, quat q) |
| 528 | { |
| 529 | float a = q[3]; |
| 530 | float b = q[0]; |
| 531 | float c = q[1]; |
| 532 | float d = q[2]; |
| 533 | float a2 = a*a; |
| 534 | float b2 = b*b; |
| 535 | float c2 = c*c; |
| 536 | float d2 = d*d; |
| 537 | |
| 538 | M[0][0] = a2 + b2 - c2 - d2; |
| 539 | M[0][1] = 2.f*(b*c + a*d); |
| 540 | M[0][2] = 2.f*(b*d - a*c); |
| 541 | M[0][3] = 0.f; |
| 542 | |
| 543 | M[1][0] = 2*(b*c - a*d); |
| 544 | M[1][1] = a2 - b2 + c2 - d2; |
| 545 | M[1][2] = 2.f*(c*d + a*b); |
| 546 | M[1][3] = 0.f; |
| 547 | |
| 548 | M[2][0] = 2.f*(b*d + a*c); |
| 549 | M[2][1] = 2.f*(c*d - a*b); |
| 550 | M[2][2] = a2 - b2 - c2 + d2; |
| 551 | M[2][3] = 0.f; |
| 552 | |
| 553 | M[3][0] = M[3][1] = M[3][2] = 0.f; |
| 554 | M[3][3] = 1.f; |
| 555 | } |
| 556 | |
| 557 | static inline void mat4x4o_mul_quat(mat4x4 R, mat4x4 M, quat q) |
| 558 | { |
| 559 | /* XXX: The way this is written only works for othogonal matrices. */ |
| 560 | /* TODO: Take care of non-orthogonal case. */ |
| 561 | quat_mul_vec3(R[0], q, M[0]); |
| 562 | quat_mul_vec3(R[1], q, M[1]); |
| 563 | quat_mul_vec3(R[2], q, M[2]); |
| 564 | |
| 565 | R[3][0] = R[3][1] = R[3][2] = 0.f; |
| 566 | R[3][3] = 1.f; |
| 567 | } |
| 568 | static inline void quat_from_mat4x4(quat q, mat4x4 M) |
| 569 | { |
| 570 | float r=0.f; |
| 571 | int i; |
| 572 | |
| 573 | int perm[] = { 0, 1, 2, 0, 1 }; |
| 574 | int *p = perm; |
| 575 | |
| 576 | for(i = 0; i<3; i++) { |
| 577 | float m = M[i][i]; |
| 578 | if( m < r ) |
| 579 | continue; |
| 580 | m = r; |
| 581 | p = &perm[i]; |
| 582 | } |
| 583 | |
| 584 | r = sqrtf(1.f + M[p[0]][p[0]] - M[p[1]][p[1]] - M[p[2]][p[2]] ); |
| 585 | |
| 586 | if(r < 1e-6) { |
| 587 | q[0] = 1.f; |
| 588 | q[1] = q[2] = q[3] = 0.f; |
| 589 | return; |
| 590 | } |
| 591 | |
| 592 | q[0] = r/2.f; |
| 593 | q[1] = (M[p[0]][p[1]] - M[p[1]][p[0]])/(2.f*r); |
| 594 | q[2] = (M[p[2]][p[0]] - M[p[0]][p[2]])/(2.f*r); |
| 595 | q[3] = (M[p[2]][p[1]] - M[p[1]][p[2]])/(2.f*r); |
| 596 | } |
| 597 | |
| 598 | #endif |