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 + bz1], rx0, ry0, rz1);
 190         v = dotgrad3(grad3[b10 + bz1], rx1, ry0, rz1);
 191         a = lerp(u, v, sx);
 192
 193         u = dotgrad3(grad3[b01 + bz1], rx0, ry1, rz1);
 194         v = dotgrad3(grad3[b11 + bz1], 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 + bz1], rx0, ry0, rz1);
 302         v = dotgrad3(grad3[b10 + bz1], rx1, ry0, rz1);
 303         a = lerp(u, v, sx);
 304
 305         u = dotgrad3(grad3[b01 + bz1], rx0, ry1, rz1);
 306         v = dotgrad3(grad3[b11 + bz1], 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         int i;
 324         float res = 0.0f, freq = 1.0f;
 325         for(i=0; i<octaves; i++) {
 326                 res += noise1(x * freq) / freq;
 327                 freq *= 2.0f;
 328         }
 329         return res;
 330 }
 331
 332 float fbm2(float x, float y, int octaves)
 333 {
 334         int i;
 335         float res = 0.0f, freq = 1.0f;
 336         for(i=0; i<octaves; i++) {
 337                 res += noise2(x * freq, y * freq) / freq;
 338                 freq *= 2.0f;
 339         }
 340         return res;
 341 }
 342
 343 float fbm3(float x, float y, float z, int octaves)
 344 {
 345         int i;
 346         float res = 0.0f, freq = 1.0f;
 347         for(i=0; i<octaves; i++) {
 348                 res += noise3(x * freq, y * freq, z * freq) / freq;
 349                 freq *= 2.0f;
 350         }
 351         return res;
 352
 353 }
 354
 355 float fbm4(float x, float y, float z, float w, int octaves)
 356 {
 357         int i;
 358         float res = 0.0f, freq = 1.0f;
 359         for(i=0; i<octaves; i++) {
 360                 res += noise4(x * freq, y * freq, z * freq, w * freq) / freq;
 361                 freq *= 2.0f;
 362         }
 363         return res;
 364 }
 365
 366
 367 float pfbm1(float x, int per, int octaves)
 368 {
 369         int i;
 370         float res = 0.0f, freq = 1.0f;
 371         for(i=0; i<octaves; i++) {
 372                 res += pnoise1(x * freq, per) / freq;
 373                 freq *= 2.0f;
 374                 per *= 2;
 375         }
 376         return res;
 377 }
 378
 379 float pfbm2(float x, float y, int per_x, int per_y, int octaves)
 380 {
 381         int i;
 382         float res = 0.0f, freq = 1.0f;
 383         for(i=0; i<octaves; i++) {
 384                 res += pnoise2(x * freq, y * freq, per_x, per_y) / freq;
 385                 freq *= 2.0f;
 386                 per_x *= 2;
 387                 per_y *= 2;
 388         }
 389         return res;
 390 }
 391
 392 float pfbm3(float x, float y, float z, int per_x, int per_y, int per_z, int octaves)
 393 {
 394         int i;
 395         float res = 0.0f, freq = 1.0f;
 396         for(i=0; i<octaves; i++) {
 397                 res += pnoise3(x * freq, y * freq, z * freq, per_x, per_y, per_z) / freq;
 398                 freq *= 2.0f;
 399                 per_x *= 2;
 400                 per_y *= 2;
 401                 per_z *= 2;
 402         }
 403         return res;
 404 }
 405
 406 float pfbm4(float x, float y, float z, float w, int per_x, int per_y, int per_z, int per_w, int octaves)
 407 {
 408         int i;
 409         float res = 0.0f, freq = 1.0f;
 410         for(i=0; i<octaves; i++) {
 411                 res += pnoise4(x * freq, y * freq, z * freq, w * freq,
 412                                 per_x, per_y, per_z, per_w) / freq;
 413                 freq *= 2.0f;
 414                 per_x *= 2;
 415                 per_y *= 2;
 416                 per_z *= 2;
 417                 per_w *= 2;
 418         }
 419         return res;
 420 }
 421
 422
 423 float turbulence1(float x, int octaves)
 424 {
 425         int i;
 426         float res = 0.0f, freq = 1.0f;
 427         for(i=0; i<octaves; i++) {
 428                 res += fabs(noise1(x * freq) / freq);
 429                 freq *= 2.0f;
 430         }
 431         return res;
 432 }
 433
 434 float turbulence2(float x, float y, int octaves)
 435 {
 436         int i;
 437         float res = 0.0f, freq = 1.0f;
 438         for(i=0; i<octaves; i++) {
 439                 res += fabs(noise2(x * freq, y * freq) / freq);
 440                 freq *= 2.0f;
 441         }
 442         return res;
 443 }
 444
 445 float turbulence3(float x, float y, float z, int octaves)
 446 {
 447         int i;
 448         float res = 0.0f, freq = 1.0f;
 449         for(i=0; i<octaves; i++) {
 450                 res += fabs(noise3(x * freq, y * freq, z * freq) / freq);
 451                 freq *= 2.0f;
 452         }
 453         return res;
 454 }
 455
 456 float turbulence4(float x, float y, float z, float w, int octaves)
 457 {
 458         int i;
 459         float res = 0.0f, freq = 1.0f;
 460         for(i=0; i<octaves; i++) {
 461                 res += fabs(noise4(x * freq, y * freq, z * freq, w * freq) / freq);
 462                 freq *= 2.0f;
 463         }
 464         return res;
 465 }
 466
 467
 468 float pturbulence1(float x, int per, int octaves)
 469 {
 470         int i;
 471         float res = 0.0f, freq = 1.0f;
 472         for(i=0; i<octaves; i++) {
 473                 res += fabs(pnoise1(x * freq, per) / freq);
 474                 freq *= 2.0f;
 475                 per *= 2;
 476         }
 477         return res;
 478 }
 479
 480 float pturbulence2(float x, float y, int per_x, int per_y, int octaves)
 481 {
 482         int i;
 483         float res = 0.0f, freq = 1.0f;
 484         for(i=0; i<octaves; i++) {
 485                 res += fabs(pnoise2(x * freq, y * freq, per_x, per_y) / freq);
 486                 freq *= 2.0f;
 487                 per_x *= 2;
 488                 per_y *= 2;
 489         }
 490         return res;
 491 }
 492
 493 float pturbulence3(float x, float y, float z, int per_x, int per_y, int per_z, int octaves)
 494 {
 495         int i;
 496         float res = 0.0f, freq = 1.0f;
 497         for(i=0; i<octaves; i++) {
 498                 res += fabs(pnoise3(x * freq, y * freq, z * freq, per_x, per_y, per_z) / freq);
 499                 freq *= 2.0f;
 500                 per_x *= 2;
 501                 per_y *= 2;
 502                 per_z *= 2;
 503         }
 504         return res;
 505 }
 506
 507 float pturbulence4(float x, float y, float z, float w, int per_x, int per_y, int per_z, int per_w, int octaves)
 508 {
 509         int i;
 510         float res = 0.0f, freq = 1.0f;
 511         for(i=0; i<octaves; i++) {
 512                 res += fabs(pnoise4(x * freq, y * freq, z * freq, w * freq,
 513                                 per_x, per_y, per_z, per_w) / freq);
 514                 freq *= 2.0f;
 515                 per_x *= 2;
 516                 per_y *= 2;
 517                 per_z *= 2;
 518                 per_w *= 2;
 519         }
 520         return res;
 521 }