git.mutantstargoat.com Git - dosdemo/blob - src/noise.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include "noise.h"
   4
   5 /* ---- Ken Perlin's implementation of noise ---- */
   6 #define B       0x100
   7 #define BM      0xff
   8 #define N       0x1000
   9 #define NP      12   /* 2^N */
  10 #define NM      0xfff
  11
  12 #define s_curve(t) (t * t * (3.0f - 2.0f * t))
  13
  14 #define setup(elem, b0, b1, r0, r1) \
  15         do {                                                    \
  16                 float t = elem + N;             \
  17                 b0 = ((int)t) & BM;                     \
  18                 b1 = (b0 + 1) & BM;                     \
  19                 r0 = t - (int)t;                        \
  20                 r1 = r0 - 1.0f;                         \
  21         } while(0)
  22
  23 #define setup_p(elem, b0, b1, r0, r1, p) \
  24         do {                                                    \
  25                 float t = elem + N;             \
  26                 b0 = (((int)t) & BM) % p;       \
  27                 b1 = ((b0 + 1) & BM) % p;       \
  28                 r0 = t - (int)t;                        \
  29                 r1 = r0 - 1.0f;                         \
  30         } while(0)
  31
  32
  33 static int perm[B + B + 2];                     /* permuted index from g_n onto themselves */
  34 static float grad3[B + B + 2][3];               /* 3D random gradients */
  35 static float grad2[B + B + 2][2];               /* 2D random gradients */
  36 static float grad1[B + B + 2];          /* 1D random ... slopes */
  37 static int tables_valid;
  38
  39 #define init_once()     if(!tables_valid) init_noise()
  40
  41 static void init_noise()
  42 {
  43         int i;
  44         float len;
  45
  46         /* calculate random gradients */
  47         for(i=0; i<B; i++) {
  48                 perm[i] = i;    /* .. and initialize permutation mapping to identity */
  49
  50                 grad1[i] = (float)((rand() % (B + B)) - B) / B;
  51
  52                 grad2[i][0] = (float)((rand() % (B + B)) - B) / B;
  53                 grad2[i][1] = (float)((rand() % (B + B)) - B) / B;
  54                 if((len = sqrt(grad2[i][0] * grad2[i][0] + grad2[i][1] * grad2[i][1])) != 0.0f) {
  55                         grad2[i][0] /= len;
  56                         grad2[i][1] /= len;
  57                 }
  58
  59                 grad3[i][0] = (float)((rand() % (B + B)) - B) / B;
  60                 grad3[i][1] = (float)((rand() % (B + B)) - B) / B;
  61                 grad3[i][2] = (float)((rand() % (B + B)) - B) / B;
  62                 if((len = sqrt(grad3[i][0] * grad3[i][0] + grad3[i][1] * grad3[i][1]) + grad3[i][2] * grad3[i][2]) != 0.0f) {
  63                         grad3[i][0] /= len;
  64                         grad3[i][1] /= len;
  65                         grad3[i][2] /= len;
  66                 }
  67         }
  68
  69         /* permute indices by swapping them randomly */
  70         for(i=0; i<B; i++) {
  71                 int rand_idx = rand() % B;
  72
  73                 int tmp = perm[i];
  74                 perm[i] = perm[rand_idx];
  75                 perm[rand_idx] = tmp;
  76         }
  77
  78         /* fill up the rest of the arrays by duplicating the existing gradients */
  79         /* and permutations */
  80         for(i=0; i<B+2; i++) {
  81                 perm[B + i] = perm[i];
  82                 grad1[B + i] = grad1[i];
  83                 grad2[B + i][0] = grad2[i][0];
  84                 grad2[B + i][1] = grad2[i][1];
  85                 grad3[B + i][0] = grad3[i][0];
  86                 grad3[B + i][1] = grad3[i][1];
  87                 grad3[B + i][2] = grad3[i][2];
  88         }
  89
  90         tables_valid = 1;
  91 }
  92
  93 #define lerp(a, b, t) ((a) + ((b) - (a)) * (t))
  94 #define dotgrad2(g, x, y) ((g)[0] * (x) + (g)[1] * (y))
  95 #define dotgrad3(g, x, y, z) ((g)[0] * (x) + (g)[1] * (y) + (g)[2] * (z))
  96
  97 float noise1(float x)
  98 {
  99         int bx0, bx1;
 100         float rx0, rx1, sx, u, v;
 101
 102         init_once();
 103
 104         setup(x, bx0, bx1, rx0, rx1);
 105         sx = s_curve(rx0);
 106         u = rx0 * grad1[perm[bx0]];
 107         v = rx1 * grad1[perm[bx1]];
 108         return lerp(u, v, sx);
 109 }
 110
 111 float noise2(float x, float y)
 112 {
 113         int i, j, b00, b10, b01, b11;
 114         int bx0, bx1, by0, by1;
 115         float rx0, rx1, ry0, ry1;
 116         float sx, sy, u, v, a, b;
 117
 118         init_once();
 119
 120         setup(x, bx0, bx1, rx0, rx1);
 121         setup(y, by0, by1, ry0, ry1);
 122
 123         i = perm[bx0];
 124         j = perm[bx1];
 125
 126         b00 = perm[i + by0];
 127         b10 = perm[j + by0];
 128         b01 = perm[i + by1];
 129         b11 = perm[j + by1];
 130
 131         /* calculate hermite inteprolating factors */
 132         sx = s_curve(rx0);
 133         sy = s_curve(ry0);
 134
 135         /* interpolate along the left edge */
 136         u = dotgrad2(grad2[b00], rx0, ry0);
 137         v = dotgrad2(grad2[b10], rx1, ry0);
 138         a = lerp(u, v, sx);
 139
 140         /* interpolate along the right edge */
 141         u = dotgrad2(grad2[b01], rx0, ry1);
 142         v = dotgrad2(grad2[b11], rx1, ry1);
 143         b = lerp(u, v, sx);
 144
 145         /* interpolate between them */
 146         return lerp(a, b, sy);
 147 }
 148
 149 float noise3(float x, float y, float z)
 150 {
 151         int i, j;
 152         int bx0, bx1, by0, by1, bz0, bz1;
 153         int b00, b10, b01, b11;
 154         float rx0, rx1, ry0, ry1, rz0, rz1;
 155         float sx, sy, sz;
 156         float u, v, a, b, c, d;
 157
 158         init_once();
 159
 160         setup(x, bx0, bx1, rx0, rx1);
 161         setup(y, by0, by1, ry0, ry1);
 162         setup(z, bz0, bz1, rz0, rz1);
 163
 164         i = perm[bx0];
 165         j = perm[bx1];
 166
 167         b00 = perm[i + by0];
 168         b10 = perm[j + by0];
 169         b01 = perm[i + by1];
 170         b11 = perm[j + by1];
 171
 172         /* calculate hermite interpolating factors */
 173         sx = s_curve(rx0);
 174         sy = s_curve(ry0);
 175         sz = s_curve(rz0);
 176
 177         /* interpolate along the top slice of the cell */
 178         u = dotgrad3(grad3[b00 + bz0], rx0, ry0, rz0);
 179         v = dotgrad3(grad3[b10 + bz0], rx1, ry0, rz0);
 180         a = lerp(u, v, sx);
 181
 182         u = dotgrad3(grad3[b01 + bz0], rx0, ry1, rz0);
 183         v = dotgrad3(grad3[b11 + bz0], rx1, ry1, rz0);
 184         b = lerp(u, v, sx);
 185
 186         c = lerp(a, b, sy);
 187
 188         /* interpolate along the bottom slice of the cell */
 189         u = dotgrad3(grad3[b00 + bz0], rx0, ry0, rz1);
 190         v = dotgrad3(grad3[b10 + bz0], rx1, ry0, rz1);
 191         a = lerp(u, v, sx);
 192
 193         u = dotgrad3(grad3[b01 + bz0], rx0, ry1, rz1);
 194         v = dotgrad3(grad3[b11 + bz0], rx1, ry1, rz1);
 195         b = lerp(u, v, sx);
 196
 197         d = lerp(a, b, sy);
 198
 199         /* interpolate between slices */
 200         return lerp(c, d, sz);
 201 }
 202
 203 float noise4(float x, float y, float z, float w)
 204 {
 205         return 0;       /* TODO */
 206 }
 207
 208
 209 float pnoise1(float x, int period)
 210 {
 211         int bx0, bx1;
 212         float rx0, rx1, sx, u, v;
 213
 214         init_once();
 215
 216         setup_p(x, bx0, bx1, rx0, rx1, period);
 217         sx = s_curve(rx0);
 218         u = rx0 * grad1[perm[bx0]];
 219         v = rx1 * grad1[perm[bx1]];
 220         return lerp(u, v, sx);
 221 }
 222
 223 float pnoise2(float x, float y, int per_x, int per_y)
 224 {
 225         int i, j, b00, b10, b01, b11;
 226         int bx0, bx1, by0, by1;
 227         float rx0, rx1, ry0, ry1;
 228         float sx, sy, u, v, a, b;
 229
 230         init_once();
 231
 232         setup_p(x, bx0, bx1, rx0, rx1, per_x);
 233         setup_p(y, by0, by1, ry0, ry1, per_y);
 234
 235         i = perm[bx0];
 236         j = perm[bx1];
 237
 238         b00 = perm[i + by0];
 239         b10 = perm[j + by0];
 240         b01 = perm[i + by1];
 241         b11 = perm[j + by1];
 242
 243         /* calculate hermite inteprolating factors */
 244         sx = s_curve(rx0);
 245         sy = s_curve(ry0);
 246
 247         /* interpolate along the left edge */
 248         u = dotgrad2(grad2[b00], rx0, ry0);
 249         v = dotgrad2(grad2[b10], rx1, ry0);
 250         a = lerp(u, v, sx);
 251
 252         /* interpolate along the right edge */
 253         u = dotgrad2(grad2[b01], rx0, ry1);
 254         v = dotgrad2(grad2[b11], rx1, ry1);
 255         b = lerp(u, v, sx);
 256
 257         /* interpolate between them */
 258         return lerp(a, b, sy);
 259 }
 260
 261 float pnoise3(float x, float y, float z, int per_x, int per_y, int per_z)
 262 {
 263         int i, j;
 264         int bx0, bx1, by0, by1, bz0, bz1;
 265         int b00, b10, b01, b11;
 266         float rx0, rx1, ry0, ry1, rz0, rz1;
 267         float sx, sy, sz;
 268         float u, v, a, b, c, d;
 269
 270         init_once();
 271
 272         setup_p(x, bx0, bx1, rx0, rx1, per_x);
 273         setup_p(y, by0, by1, ry0, ry1, per_y);
 274         setup_p(z, bz0, bz1, rz0, rz1, per_z);
 275
 276         i = perm[bx0];
 277         j = perm[bx1];
 278
 279         b00 = perm[i + by0];
 280         b10 = perm[j + by0];
 281         b01 = perm[i + by1];
 282         b11 = perm[j + by1];
 283
 284         /* calculate hermite interpolating factors */
 285         sx = s_curve(rx0);
 286         sy = s_curve(ry0);
 287         sz = s_curve(rz0);
 288
 289         /* interpolate along the top slice of the cell */
 290         u = dotgrad3(grad3[b00 + bz0], rx0, ry0, rz0);
 291         v = dotgrad3(grad3[b10 + bz0], rx1, ry0, rz0);
 292         a = lerp(u, v, sx);
 293
 294         u = dotgrad3(grad3[b01 + bz0], rx0, ry1, rz0);
 295         v = dotgrad3(grad3[b11 + bz0], rx1, ry1, rz0);
 296         b = lerp(u, v, sx);
 297
 298         c = lerp(a, b, sy);
 299
 300         /* interpolate along the bottom slice of the cell */
 301         u = dotgrad3(grad3[b00 + bz0], rx0, ry0, rz1);
 302         v = dotgrad3(grad3[b10 + bz0], rx1, ry0, rz1);
 303         a = lerp(u, v, sx);
 304
 305         u = dotgrad3(grad3[b01 + bz0], rx0, ry1, rz1);
 306         v = dotgrad3(grad3[b11 + bz0], rx1, ry1, rz1);
 307         b = lerp(u, v, sx);
 308
 309         d = lerp(a, b, sy);
 310
 311         /* interpolate between slices */
 312         return lerp(c, d, sz);
 313 }
 314
 315 float pnoise4(float x, float y, float z, float w, int per_x, int per_y, int per_z, int per_w)
 316 {
 317         return 0;
 318 }
 319
 320
 321 float fbm1(float x, int octaves)
 322 {
 323         float res = 0.0f, freq = 1.0f;
 324         for(int i=0; i<octaves; i++) {
 325                 res += noise1(x * freq) / freq;
 326                 freq *= 2.0f;
 327         }
 328         return res;
 329 }
 330
 331 float fbm2(float x, float y, int octaves)
 332 {
 333         float res = 0.0f, freq = 1.0f;
 334         for(int i=0; i<octaves; i++) {
 335                 res += noise2(x * freq, y * freq) / freq;
 336                 freq *= 2.0f;
 337         }
 338         return res;
 339 }
 340
 341 float fbm3(float x, float y, float z, int octaves)
 342 {
 343         float res = 0.0f, freq = 1.0f;
 344         for(int i=0; i<octaves; i++) {
 345                 res += noise3(x * freq, y * freq, z * freq) / freq;
 346                 freq *= 2.0f;
 347         }
 348         return res;
 349
 350 }
 351
 352 float fbm4(float x, float y, float z, float w, int octaves)
 353 {
 354         float res = 0.0f, freq = 1.0f;
 355         for(int i=0; i<octaves; i++) {
 356                 res += noise4(x * freq, y * freq, z * freq, w * freq) / freq;
 357                 freq *= 2.0f;
 358         }
 359         return res;
 360 }
 361
 362
 363 float pfbm1(float x, int per, int octaves)
 364 {
 365         float res = 0.0f, freq = 1.0f;
 366         for(int i=0; i<octaves; i++) {
 367                 res += pnoise1(x * freq, per) / freq;
 368                 freq *= 2.0f;
 369                 per *= 2;
 370         }
 371         return res;
 372 }
 373
 374 float pfbm2(float x, float y, int per_x, int per_y, int octaves)
 375 {
 376         float res = 0.0f, freq = 1.0f;
 377         for(int i=0; i<octaves; i++) {
 378                 res += pnoise2(x * freq, y * freq, per_x, per_y) / freq;
 379                 freq *= 2.0f;
 380                 per_x *= 2;
 381                 per_y *= 2;
 382         }
 383         return res;
 384 }
 385
 386 float pfbm3(float x, float y, float z, int per_x, int per_y, int per_z, int octaves)
 387 {
 388         float res = 0.0f, freq = 1.0f;
 389         for(int i=0; i<octaves; i++) {
 390                 res += pnoise3(x * freq, y * freq, z * freq, per_x, per_y, per_z) / freq;
 391                 freq *= 2.0f;
 392                 per_x *= 2;
 393                 per_y *= 2;
 394                 per_z *= 2;
 395         }
 396         return res;
 397 }
 398
 399 float pfbm4(float x, float y, float z, float w, int per_x, int per_y, int per_z, int per_w, int octaves)
 400 {
 401         float res = 0.0f, freq = 1.0f;
 402         for(int i=0; i<octaves; i++) {
 403                 res += pnoise4(x * freq, y * freq, z * freq, w * freq,
 404                                 per_x, per_y, per_z, per_w) / freq;
 405                 freq *= 2.0f;
 406                 per_x *= 2;
 407                 per_y *= 2;
 408                 per_z *= 2;
 409                 per_w *= 2;
 410         }
 411         return res;
 412 }
 413
 414
 415 float turbulence1(float x, int octaves)
 416 {
 417         float res = 0.0f, freq = 1.0f;
 418         for(int i=0; i<octaves; i++) {
 419                 res += fabs(noise1(x * freq) / freq);
 420                 freq *= 2.0f;
 421         }
 422         return res;
 423 }
 424
 425 float turbulence2(float x, float y, int octaves)
 426 {
 427         float res = 0.0f, freq = 1.0f;
 428         for(int i=0; i<octaves; i++) {
 429                 res += fabs(noise2(x * freq, y * freq) / freq);
 430                 freq *= 2.0f;
 431         }
 432         return res;
 433 }
 434
 435 float turbulence3(float x, float y, float z, int octaves)
 436 {
 437         float res = 0.0f, freq = 1.0f;
 438         for(int i=0; i<octaves; i++) {
 439                 res += fabs(noise3(x * freq, y * freq, z * freq) / freq);
 440                 freq *= 2.0f;
 441         }
 442         return res;
 443 }
 444
 445 float turbulence4(float x, float y, float z, float w, int octaves)
 446 {
 447         float res = 0.0f, freq = 1.0f;
 448         for(int i=0; i<octaves; i++) {
 449                 res += fabs(noise4(x * freq, y * freq, z * freq, w * freq) / freq);
 450                 freq *= 2.0f;
 451         }
 452         return res;
 453 }
 454
 455
 456 float pturbulence1(float x, int per, int octaves)
 457 {
 458         float res = 0.0f, freq = 1.0f;
 459         for(int i=0; i<octaves; i++) {
 460                 res += fabs(pnoise1(x * freq, per) / freq);
 461                 freq *= 2.0f;
 462                 per *= 2;
 463         }
 464         return res;
 465 }
 466
 467 float pturbulence2(float x, float y, int per_x, int per_y, int octaves)
 468 {
 469         float res = 0.0f, freq = 1.0f;
 470         for(int i=0; i<octaves; i++) {
 471                 res += fabs(pnoise2(x * freq, y * freq, per_x, per_y) / freq);
 472                 freq *= 2.0f;
 473                 per_x *= 2;
 474                 per_y *= 2;
 475         }
 476         return res;
 477 }
 478
 479 float pturbulence3(float x, float y, float z, int per_x, int per_y, int per_z, int octaves)
 480 {
 481         float res = 0.0f, freq = 1.0f;
 482         for(int i=0; i<octaves; i++) {
 483                 res += fabs(pnoise3(x * freq, y * freq, z * freq, per_x, per_y, per_z) / freq);
 484                 freq *= 2.0f;
 485                 per_x *= 2;
 486                 per_y *= 2;
 487                 per_z *= 2;
 488         }
 489         return res;
 490 }
 491
 492 float pturbulence4(float x, float y, float z, float w, int per_x, int per_y, int per_z, int per_w, int octaves)
 493 {
 494         float res = 0.0f, freq = 1.0f;
 495         for(int i=0; i<octaves; i++) {
 496                 res += fabs(pnoise4(x * freq, y * freq, z * freq, w * freq,
 497                                 per_x, per_y, per_z, per_w) / freq);
 498                 freq *= 2.0f;
 499                 per_x *= 2;
 500                 per_y *= 2;
 501                 per_z *= 2;
 502                 per_w *= 2;
 503         }
 504         return res;
 505 }