3dphotoshoot

annotate libs/vmath/vmath.c @ 13:8cfe061ddd5f

find changesets by author, revision, files, or words in the commit message
fixed the assets directories in the android makefile
author John Tsiombikas <nuclear@member.fsf.org>
date Sun, 07 Jun 2015 02:29:10 +0300
parents
children
rev   line source
nuclear@10 1 /*
nuclear@10 2 libvmath - a vector math library
nuclear@10 3 Copyright (C) 2004-2015 John Tsiombikas <nuclear@member.fsf.org>
nuclear@10 4
nuclear@10 5 This program is free software: you can redistribute it and/or modify
nuclear@10 6 it under the terms of the GNU Lesser General Public License as published
nuclear@10 7 by the Free Software Foundation, either version 3 of the License, or
nuclear@10 8 (at your option) any later version.
nuclear@10 9
nuclear@10 10 This program is distributed in the hope that it will be useful,
nuclear@10 11 but WITHOUT ANY WARRANTY; without even the implied warranty of
nuclear@10 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
nuclear@10 13 GNU Lesser General Public License for more details.
nuclear@10 14
nuclear@10 15 You should have received a copy of the GNU Lesser General Public License
nuclear@10 16 along with this program. If not, see <http://www.gnu.org/licenses/>.
nuclear@10 17 */
nuclear@10 18 #include <stdlib.h>
nuclear@10 19 #include <math.h>
nuclear@10 20 #include "vmath.h"
nuclear@10 21
nuclear@10 22 #if defined(__APPLE__) && !defined(TARGET_IPHONE)
nuclear@10 23 #include <xmmintrin.h>
nuclear@10 24
nuclear@10 25 void enable_fpexcept(void)
nuclear@10 26 {
nuclear@10 27 unsigned int bits;
nuclear@10 28 bits = _MM_MASK_INVALID | _MM_MASK_DIV_ZERO | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW;
nuclear@10 29 _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~bits);
nuclear@10 30 }
nuclear@10 31
nuclear@10 32 void disable_fpexcept(void)
nuclear@10 33 {
nuclear@10 34 unsigned int bits;
nuclear@10 35 bits = _MM_MASK_INVALID | _MM_MASK_DIV_ZERO | _MM_MASK_OVERFLOW | _MM_MASK_UNDERFLOW;
nuclear@10 36 _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() | bits);
nuclear@10 37 }
nuclear@10 38
nuclear@10 39 #elif defined(__GNUC__) && !defined(TARGET_IPHONE) && !defined(__MINGW32__)
nuclear@10 40 #define __USE_GNU
nuclear@10 41 #include <fenv.h>
nuclear@10 42
nuclear@10 43 void enable_fpexcept(void)
nuclear@10 44 {
nuclear@10 45 feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW);
nuclear@10 46 }
nuclear@10 47
nuclear@10 48 void disable_fpexcept(void)
nuclear@10 49 {
nuclear@10 50 fedisableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW);
nuclear@10 51 }
nuclear@10 52
nuclear@10 53 #elif defined(_MSC_VER) || defined(__MINGW32__)
nuclear@10 54 #include <float.h>
nuclear@10 55
nuclear@10 56 #if defined(__MINGW32__) && !defined(_EM_OVERFLOW)
nuclear@10 57 /* if gcc's float.h gets precedence, the mingw MSVC includes won't be declared */
nuclear@10 58 #define _MCW_EM 0x8001f
nuclear@10 59 #define _EM_INVALID 0x10
nuclear@10 60 #define _EM_ZERODIVIDE 0x08
nuclear@10 61 #define _EM_OVERFLOW 0x04
nuclear@10 62 unsigned int __cdecl _clearfp(void);
nuclear@10 63 unsigned int __cdecl _controlfp(unsigned int, unsigned int);
nuclear@10 64 #endif
nuclear@10 65
nuclear@10 66 void enable_fpexcept(void)
nuclear@10 67 {
nuclear@10 68 _clearfp();
nuclear@10 69 _controlfp(_controlfp(0, 0) & ~(_EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW), _MCW_EM);
nuclear@10 70 }
nuclear@10 71
nuclear@10 72 void disable_fpexcept(void)
nuclear@10 73 {
nuclear@10 74 _clearfp();
nuclear@10 75 _controlfp(_controlfp(0, 0) | (_EM_INVALID | _EM_ZERODIVIDE | _EM_OVERFLOW), _MCW_EM);
nuclear@10 76 }
nuclear@10 77 #else
nuclear@10 78 void enable_fpexcept(void) {}
nuclear@10 79 void disable_fpexcept(void) {}
nuclear@10 80 #endif
nuclear@10 81
nuclear@10 82
nuclear@10 83 /** Numerical calculation of integrals using simpson's rule */
nuclear@10 84 scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples)
nuclear@10 85 {
nuclear@10 86 int i;
nuclear@10 87 scalar_t h = (high - low) / (scalar_t)samples;
nuclear@10 88 scalar_t sum = 0.0;
nuclear@10 89
nuclear@10 90 for(i=0; i<samples+1; i++) {
nuclear@10 91 scalar_t y = f((scalar_t)i * h + low);
nuclear@10 92 sum += ((!i || i == samples) ? y : ((i % 2) ? 4.0 * y : 2.0 * y)) * (h / 3.0);
nuclear@10 93 }
nuclear@10 94 return sum;
nuclear@10 95 }
nuclear@10 96
nuclear@10 97 /** Gaussuan function */
nuclear@10 98 scalar_t gaussian(scalar_t x, scalar_t mean, scalar_t sdev)
nuclear@10 99 {
nuclear@10 100 scalar_t exponent = -SQ(x - mean) / (2.0 * SQ(sdev));
nuclear@10 101 return 1.0 - -pow(M_E, exponent) / (sdev * sqrt(TWO_PI));
nuclear@10 102 }
nuclear@10 103
nuclear@10 104
nuclear@10 105 /** b-spline approximation */
nuclear@10 106 scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)
nuclear@10 107 {
nuclear@10 108 vec4_t tmp;
nuclear@10 109 scalar_t tsq = t * t;
nuclear@10 110
nuclear@10 111 static mat4_t bspline_mat = {
nuclear@10 112 {-1, 3, -3, 1},
nuclear@10 113 {3, -6, 3, 0},
nuclear@10 114 {-3, 0, 3, 0},
nuclear@10 115 {1, 4, 1, 0}
nuclear@10 116 };
nuclear@10 117
nuclear@10 118 tmp = v4_scale(v4_transform(v4_cons(a, b, c, d), bspline_mat), 1.0 / 6.0);
nuclear@10 119 return v4_dot(v4_cons(tsq * t, tsq, t, 1.0), tmp);
nuclear@10 120 }
nuclear@10 121
nuclear@10 122 /** Catmull-rom spline interpolation */
nuclear@10 123 scalar_t spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)
nuclear@10 124 {
nuclear@10 125 vec4_t tmp;
nuclear@10 126 scalar_t tsq = t * t;
nuclear@10 127
nuclear@10 128 static mat4_t crspline_mat = {
nuclear@10 129 {-1, 3, -3, 1},
nuclear@10 130 {2, -5, 4, -1},
nuclear@10 131 {-1, 0, 1, 0},
nuclear@10 132 {0, 2, 0, 0}
nuclear@10 133 };
nuclear@10 134
nuclear@10 135 tmp = v4_scale(v4_transform(v4_cons(a, b, c, d), crspline_mat), 0.5);
nuclear@10 136 return v4_dot(v4_cons(tsq * t, tsq, t, 1.0), tmp);
nuclear@10 137 }
nuclear@10 138
nuclear@10 139 /** Bezier interpolation */
nuclear@10 140 scalar_t bezier(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)
nuclear@10 141 {
nuclear@10 142 scalar_t omt, omt3, t3, f;
nuclear@10 143 t3 = t * t * t;
nuclear@10 144 omt = 1.0f - t;
nuclear@10 145 omt3 = omt * omt * omt;
nuclear@10 146 f = 3 * t * omt;
nuclear@10 147
nuclear@10 148 return (a * omt3) + (b * f * omt) + (c * f * t) + (d * t3);
nuclear@10 149 }
nuclear@10 150
nuclear@10 151 /* ---- Ken Perlin's implementation of noise ---- */
nuclear@10 152
nuclear@10 153 #define B 0x100
nuclear@10 154 #define BM 0xff
nuclear@10 155 #define N 0x1000
nuclear@10 156 #define NP 12 /* 2^N */
nuclear@10 157 #define NM 0xfff
nuclear@10 158
nuclear@10 159 #define s_curve(t) (t * t * (3.0f - 2.0f * t))
nuclear@10 160
nuclear@10 161 #define setup(elem, b0, b1, r0, r1) \
nuclear@10 162 do { \
nuclear@10 163 scalar_t t = elem + N; \
nuclear@10 164 b0 = ((int)t) & BM; \
nuclear@10 165 b1 = (b0 + 1) & BM; \
nuclear@10 166 r0 = t - (int)t; \
nuclear@10 167 r1 = r0 - 1.0f; \
nuclear@10 168 } while(0)
nuclear@10 169
nuclear@10 170
nuclear@10 171 static int perm[B + B + 2]; /* permuted index from g_n onto themselves */
nuclear@10 172 static vec3_t grad3[B + B + 2]; /* 3D random gradients */
nuclear@10 173 static vec2_t grad2[B + B + 2]; /* 2D random gradients */
nuclear@10 174 static scalar_t grad1[B + B + 2]; /* 1D random ... slopes */
nuclear@10 175 static int tables_valid;
nuclear@10 176
nuclear@10 177 static void init_noise()
nuclear@10 178 {
nuclear@10 179 int i;
nuclear@10 180
nuclear@10 181 /* calculate random gradients */
nuclear@10 182 for(i=0; i<B; i++) {
nuclear@10 183 perm[i] = i; /* .. and initialize permutation mapping to identity */
nuclear@10 184
nuclear@10 185 grad1[i] = (scalar_t)((rand() % (B + B)) - B) / B;
nuclear@10 186
nuclear@10 187 grad2[i].x = (scalar_t)((rand() % (B + B)) - B) / B;
nuclear@10 188 grad2[i].y = (scalar_t)((rand() % (B + B)) - B) / B;
nuclear@10 189 grad2[i] = v2_normalize(grad2[i]);
nuclear@10 190
nuclear@10 191 grad3[i].x = (scalar_t)((rand() % (B + B)) - B) / B;
nuclear@10 192 grad3[i].y = (scalar_t)((rand() % (B + B)) - B) / B;
nuclear@10 193 grad3[i].z = (scalar_t)((rand() % (B + B)) - B) / B;
nuclear@10 194 grad3[i] = v3_normalize(grad3[i]);
nuclear@10 195 }
nuclear@10 196
nuclear@10 197 /* permute indices by swapping them randomly */
nuclear@10 198 for(i=0; i<B; i++) {
nuclear@10 199 int rand_idx = rand() % B;
nuclear@10 200
nuclear@10 201 int tmp = perm[i];
nuclear@10 202 perm[i] = perm[rand_idx];
nuclear@10 203 perm[rand_idx] = tmp;
nuclear@10 204 }
nuclear@10 205
nuclear@10 206 /* fill up the rest of the arrays by duplicating the existing gradients */
nuclear@10 207 /* and permutations */
nuclear@10 208 for(i=0; i<B+2; i++) {
nuclear@10 209 perm[B + i] = perm[i];
nuclear@10 210 grad1[B + i] = grad1[i];
nuclear@10 211 grad2[B + i] = grad2[i];
nuclear@10 212 grad3[B + i] = grad3[i];
nuclear@10 213 }
nuclear@10 214 }
nuclear@10 215
nuclear@10 216 scalar_t noise1(scalar_t x)
nuclear@10 217 {
nuclear@10 218 int bx0, bx1;
nuclear@10 219 scalar_t rx0, rx1, sx, u, v;
nuclear@10 220
nuclear@10 221 if(!tables_valid) {
nuclear@10 222 init_noise();
nuclear@10 223 tables_valid = 1;
nuclear@10 224 }
nuclear@10 225
nuclear@10 226 setup(x, bx0, bx1, rx0, rx1);
nuclear@10 227 sx = s_curve(rx0);
nuclear@10 228 u = rx0 * grad1[perm[bx0]];
nuclear@10 229 v = rx1 * grad1[perm[bx1]];
nuclear@10 230
nuclear@10 231 return lerp(u, v, sx);
nuclear@10 232 }
nuclear@10 233
nuclear@10 234 scalar_t noise2(scalar_t x, scalar_t y)
nuclear@10 235 {
nuclear@10 236 int i, j, b00, b10, b01, b11;
nuclear@10 237 int bx0, bx1, by0, by1;
nuclear@10 238 scalar_t rx0, rx1, ry0, ry1;
nuclear@10 239 scalar_t sx, sy, u, v, a, b;
nuclear@10 240
nuclear@10 241 if(!tables_valid) {
nuclear@10 242 init_noise();
nuclear@10 243 tables_valid = 1;
nuclear@10 244 }
nuclear@10 245
nuclear@10 246 setup(x, bx0, bx1, rx0, rx1);
nuclear@10 247 setup(y, by0, by1, ry0, ry1);
nuclear@10 248
nuclear@10 249 i = perm[bx0];
nuclear@10 250 j = perm[bx1];
nuclear@10 251
nuclear@10 252 b00 = perm[i + by0];
nuclear@10 253 b10 = perm[j + by0];
nuclear@10 254 b01 = perm[i + by1];
nuclear@10 255 b11 = perm[j + by1];
nuclear@10 256
nuclear@10 257 /* calculate hermite inteprolating factors */
nuclear@10 258 sx = s_curve(rx0);
nuclear@10 259 sy = s_curve(ry0);
nuclear@10 260
nuclear@10 261 /* interpolate along the left edge */
nuclear@10 262 u = v2_dot(grad2[b00], v2_cons(rx0, ry0));
nuclear@10 263 v = v2_dot(grad2[b10], v2_cons(rx1, ry0));
nuclear@10 264 a = lerp(u, v, sx);
nuclear@10 265
nuclear@10 266 /* interpolate along the right edge */
nuclear@10 267 u = v2_dot(grad2[b01], v2_cons(rx0, ry1));
nuclear@10 268 v = v2_dot(grad2[b11], v2_cons(rx1, ry1));
nuclear@10 269 b = lerp(u, v, sx);
nuclear@10 270
nuclear@10 271 /* interpolate between them */
nuclear@10 272 return lerp(a, b, sy);
nuclear@10 273 }
nuclear@10 274
nuclear@10 275 scalar_t noise3(scalar_t x, scalar_t y, scalar_t z)
nuclear@10 276 {
nuclear@10 277 int i, j;
nuclear@10 278 int bx0, bx1, by0, by1, bz0, bz1;
nuclear@10 279 int b00, b10, b01, b11;
nuclear@10 280 scalar_t rx0, rx1, ry0, ry1, rz0, rz1;
nuclear@10 281 scalar_t sx, sy, sz;
nuclear@10 282 scalar_t u, v, a, b, c, d;
nuclear@10 283
nuclear@10 284 if(!tables_valid) {
nuclear@10 285 init_noise();
nuclear@10 286 tables_valid = 1;
nuclear@10 287 }
nuclear@10 288
nuclear@10 289 setup(x, bx0, bx1, rx0, rx1);
nuclear@10 290 setup(y, by0, by1, ry0, ry1);
nuclear@10 291 setup(z, bz0, bz1, rz0, rz1);
nuclear@10 292
nuclear@10 293 i = perm[bx0];
nuclear@10 294 j = perm[bx1];
nuclear@10 295
nuclear@10 296 b00 = perm[i + by0];
nuclear@10 297 b10 = perm[j + by0];
nuclear@10 298 b01 = perm[i + by1];
nuclear@10 299 b11 = perm[j + by1];
nuclear@10 300
nuclear@10 301 /* calculate hermite interpolating factors */
nuclear@10 302 sx = s_curve(rx0);
nuclear@10 303 sy = s_curve(ry0);
nuclear@10 304 sz = s_curve(rz0);
nuclear@10 305
nuclear@10 306 /* interpolate along the top slice of the cell */
nuclear@10 307 u = v3_dot(grad3[b00 + bz0], v3_cons(rx0, ry0, rz0));
nuclear@10 308 v = v3_dot(grad3[b10 + bz0], v3_cons(rx1, ry0, rz0));
nuclear@10 309 a = lerp(u, v, sx);
nuclear@10 310
nuclear@10 311 u = v3_dot(grad3[b01 + bz0], v3_cons(rx0, ry1, rz0));
nuclear@10 312 v = v3_dot(grad3[b11 + bz0], v3_cons(rx1, ry1, rz0));
nuclear@10 313 b = lerp(u, v, sx);
nuclear@10 314
nuclear@10 315 c = lerp(a, b, sy);
nuclear@10 316
nuclear@10 317 /* interpolate along the bottom slice of the cell */
nuclear@10 318 u = v3_dot(grad3[b00 + bz0], v3_cons(rx0, ry0, rz1));
nuclear@10 319 v = v3_dot(grad3[b10 + bz0], v3_cons(rx1, ry0, rz1));
nuclear@10 320 a = lerp(u, v, sx);
nuclear@10 321
nuclear@10 322 u = v3_dot(grad3[b01 + bz0], v3_cons(rx0, ry1, rz1));
nuclear@10 323 v = v3_dot(grad3[b11 + bz0], v3_cons(rx1, ry1, rz1));
nuclear@10 324 b = lerp(u, v, sx);
nuclear@10 325
nuclear@10 326 d = lerp(a, b, sy);
nuclear@10 327
nuclear@10 328 /* interpolate between slices */
nuclear@10 329 return lerp(c, d, sz);
nuclear@10 330 }
nuclear@10 331
nuclear@10 332 scalar_t fbm1(scalar_t x, int octaves)
nuclear@10 333 {
nuclear@10 334 int i;
nuclear@10 335 scalar_t res = 0.0f, freq = 1.0f;
nuclear@10 336 for(i=0; i<octaves; i++) {
nuclear@10 337 res += noise1(x * freq) / freq;
nuclear@10 338 freq *= 2.0f;
nuclear@10 339 }
nuclear@10 340 return res;
nuclear@10 341 }
nuclear@10 342
nuclear@10 343 scalar_t fbm2(scalar_t x, scalar_t y, int octaves)
nuclear@10 344 {
nuclear@10 345 int i;
nuclear@10 346 scalar_t res = 0.0f, freq = 1.0f;
nuclear@10 347 for(i=0; i<octaves; i++) {
nuclear@10 348 res += noise2(x * freq, y * freq) / freq;
nuclear@10 349 freq *= 2.0f;
nuclear@10 350 }
nuclear@10 351 return res;
nuclear@10 352 }
nuclear@10 353
nuclear@10 354 scalar_t fbm3(scalar_t x, scalar_t y, scalar_t z, int octaves)
nuclear@10 355 {
nuclear@10 356 int i;
nuclear@10 357 scalar_t res = 0.0f, freq = 1.0f;
nuclear@10 358 for(i=0; i<octaves; i++) {
nuclear@10 359 res += noise3(x * freq, y * freq, z * freq) / freq;
nuclear@10 360 freq *= 2.0f;
nuclear@10 361 }
nuclear@10 362 return res;
nuclear@10 363 }
nuclear@10 364
nuclear@10 365 scalar_t turbulence1(scalar_t x, int octaves)
nuclear@10 366 {
nuclear@10 367 int i;
nuclear@10 368 scalar_t res = 0.0f, freq = 1.0f;
nuclear@10 369 for(i=0; i<octaves; i++) {
nuclear@10 370 res += fabs(noise1(x * freq) / freq);
nuclear@10 371 freq *= 2.0f;
nuclear@10 372 }
nuclear@10 373 return res;
nuclear@10 374 }
nuclear@10 375
nuclear@10 376 scalar_t turbulence2(scalar_t x, scalar_t y, int octaves)
nuclear@10 377 {
nuclear@10 378 int i;
nuclear@10 379 scalar_t res = 0.0f, freq = 1.0f;
nuclear@10 380 for(i=0; i<octaves; i++) {
nuclear@10 381 res += fabs(noise2(x * freq, y * freq) / freq);
nuclear@10 382 freq *= 2.0f;
nuclear@10 383 }
nuclear@10 384 return res;
nuclear@10 385 }
nuclear@10 386
nuclear@10 387 scalar_t turbulence3(scalar_t x, scalar_t y, scalar_t z, int octaves)
nuclear@10 388 {
nuclear@10 389 int i;
nuclear@10 390 scalar_t res = 0.0f, freq = 1.0f;
nuclear@10 391 for(i=0; i<octaves; i++) {
nuclear@10 392 res += fabs(noise3(x * freq, y * freq, z * freq) / freq);
nuclear@10 393 freq *= 2.0f;
nuclear@10 394 }
nuclear@10 395 return res;
nuclear@10 396 }