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

   1 #include <math.h>
   2 #include <stdlib.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 #define SQ(x)   ((x) * (x))
  42 #define LERP(a, b, t)   ((a) + ((b) - (a)) * (t))
  43
  44 static void init_noise()
  45 {
  46         int i;
  47         float len;
  48
  49         /* calculate random gradients */
  50         for(i=0; i<B; i++) {
  51                 perm[i] = i;    /* .. and initialize permutation mapping to identity */
  52
  53                 grad1[i] = (float)((rand() % (B + B)) - B) / B;
  54
  55                 grad2[i][0] = (float)((rand() % (B + B)) - B) / B;
  56                 grad2[i][1] = (float)((rand() % (B + B)) - B) / B;
  57                 len = sqrt(SQ(grad2[i][0]) + SQ(grad2[i][1]));
  58                 if(len != 0.0f) {
  59                         grad2[i][0] /= len;
  60                         grad2[i][1] /= len;
  61                 }
  62
  63                 grad3[i][0] = (float)((rand() % (B + B)) - B) / B;
  64                 grad3[i][1] = (float)((rand() % (B + B)) - B) / B;
  65                 grad3[i][2] = (float)((rand() % (B + B)) - B) / B;
  66                 len = sqrt(SQ(grad3[i][0]) + SQ(grad3[i][1]) + SQ(grad3[i][2]));
  67                 if(len != 0.0f) {
  68                         grad3[i][0] /= len;
  69                         grad3[i][1] /= len;
  70                         grad3[i][2] /= len;
  71                 }
  72         }
  73
  74         /* permute indices by swapping them randomly */
  75         for(i=0; i<B; i++) {
  76                 int rand_idx = rand() % B;
  77
  78                 int tmp = perm[i];
  79                 perm[i] = perm[rand_idx];
  80                 perm[rand_idx] = tmp;
  81         }
  82
  83         /* fill up the rest of the arrays by duplicating the existing gradients */
  84         /* and permutations */
  85         for(i=0; i<B+2; i++) {
  86                 perm[B + i] = perm[i];
  87                 grad1[B + i] = grad1[i];
  88                 grad2[B + i][0] = grad2[i][0];
  89                 grad2[B + i][1] = grad2[i][1];
  90                 grad3[B + i][0] = grad3[i][0];
  91                 grad3[B + i][1] = grad3[i][1];
  92                 grad3[B + i][2] = grad3[i][2];
  93         }
  94
  95         tables_valid = 1;
  96 }
  97
  98
  99 float noise1(float x)
 100 {
 101         int bx0, bx1;
 102         float rx0, rx1, sx, u, v;
 103
 104         init_once();
 105
 106         setup(x, bx0, bx1, rx0, rx1);
 107         sx = s_curve(rx0);
 108         u = rx0 * grad1[perm[bx0]];
 109         v = rx1 * grad1[perm[bx1]];
 110         return LERP(u, v, sx);
 111 }
 112
 113 float noise2(float x, float y)
 114 {
 115         int i, j, b00, b10, b01, b11;
 116         int bx0, bx1, by0, by1;
 117         float rx0, rx1, ry0, ry1;
 118         float sx, sy, u, v, a, b;
 119
 120         init_once();
 121
 122         setup(x, bx0, bx1, rx0, rx1);
 123         setup(y, by0, by1, ry0, ry1);
 124
 125         i = perm[bx0];
 126         j = perm[bx1];
 127
 128         b00 = perm[i + by0];
 129         b10 = perm[j + by0];
 130         b01 = perm[i + by1];
 131         b11 = perm[j + by1];
 132
 133         /* calculate hermite inteprolating factors */
 134         sx = s_curve(rx0);
 135         sy = s_curve(ry0);
 136
 137         /* interpolate along the left edge */
 138         u = grad2[b00][0] * rx0 + grad2[b00][1] * ry0;
 139         v = grad2[b10][0] * rx1 + grad2[b10][1] * ry0;
 140         a = LERP(u, v, sx);
 141
 142         /* interpolate along the right edge */
 143         u = grad2[b01][0] * rx0 + grad2[b01][1] * ry1;
 144         v = grad2[b11][0] * rx1 + grad2[b11][1] * ry1;
 145         b = LERP(u, v, sx);
 146
 147         /* interpolate between them */
 148         return LERP(a, b, sy);
 149 }
 150
 151 float noise3(float x, float y, float z)
 152 {
 153         int i, j;
 154         int bx0, bx1, by0, by1, bz0, bz1;
 155         int b00, b10, b01, b11;
 156         float rx0, rx1, ry0, ry1, rz0, rz1;
 157         float sx, sy, sz;
 158         float u, v, a, b, c, d;
 159
 160         init_once();
 161
 162         setup(x, bx0, bx1, rx0, rx1);
 163         setup(y, by0, by1, ry0, ry1);
 164         setup(z, bz0, bz1, rz0, rz1);
 165
 166         i = perm[bx0];
 167         j = perm[bx1];
 168
 169         b00 = perm[i + by0];
 170         b10 = perm[j + by0];
 171         b01 = perm[i + by1];
 172         b11 = perm[j + by1];
 173
 174         /* calculate hermite interpolating factors */
 175         sx = s_curve(rx0);
 176         sy = s_curve(ry0);
 177         sz = s_curve(rz0);
 178
 179         /* interpolate along the top slice of the cell */
 180         u = grad3[b00 + bz0][0] * rx0 + grad3[b00 + bz0][1] * ry0 + grad3[b00 + bz0][2] * rz0;
 181         v = grad3[b10 + bz0][0] * rx1 + grad3[b10 + bz0][1] * ry0 + grad3[b10 + bz0][2] * rz0;
 182         a = LERP(u, v, sx);
 183
 184         u = grad3[b01 + bz0][0] * rx0 + grad3[b01 + bz0][1] * ry1 + grad3[b01 + bz0][2] * rz0;
 185         v = grad3[b11 + bz0][0] * rx1 + grad3[b11 + bz0][1] * ry1 + grad3[b11 + bz0][2] * rz0;
 186         b = LERP(u, v, sx);
 187
 188         c = LERP(a, b, sy);
 189
 190         /* interpolate along the bottom slice of the cell */
 191         u = grad3[b00 + bz0][0] * rx0 + grad3[b00 + bz0][1] * ry0 + grad3[b00 + bz0][2] * rz1;
 192         v = grad3[b10 + bz0][0] * rx1 + grad3[b10 + bz0][1] * ry0 + grad3[b10 + bz0][2] * rz1;
 193         a = LERP(u, v, sx);
 194
 195         u = grad3[b01 + bz0][0] * rx0 + grad3[b01 + bz0][1] * ry1 + grad3[b01 + bz0][2] * rz1;
 196         v = grad3[b11 + bz0][0] * rx1 + grad3[b11 + bz0][1] * ry1 + grad3[b11 + bz0][2] * rz1;
 197         b = LERP(u, v, sx);
 198
 199         d = LERP(a, b, sy);
 200
 201         /* interpolate between slices */
 202         return LERP(c, d, sz);
 203 }
 204
 205
 206 float pnoise1(float x, int period)
 207 {
 208         int bx0, bx1;
 209         float rx0, rx1, sx, u, v;
 210
 211         init_once();
 212
 213         setup_p(x, bx0, bx1, rx0, rx1, period);
 214         sx = s_curve(rx0);
 215         u = rx0 * grad1[perm[bx0]];
 216         v = rx1 * grad1[perm[bx1]];
 217         return LERP(u, v, sx);
 218 }
 219
 220 float pnoise2(float x, float y, int per_x, int per_y)
 221 {
 222         int i, j, b00, b10, b01, b11;
 223         int bx0, bx1, by0, by1;
 224         float rx0, rx1, ry0, ry1;
 225         float sx, sy, u, v, a, b;
 226
 227         init_once();
 228
 229         setup_p(x, bx0, bx1, rx0, rx1, per_x);
 230         setup_p(y, by0, by1, ry0, ry1, per_y);
 231
 232         i = perm[bx0];
 233         j = perm[bx1];
 234
 235         b00 = perm[i + by0];
 236         b10 = perm[j + by0];
 237         b01 = perm[i + by1];
 238         b11 = perm[j + by1];
 239
 240         /* calculate hermite inteprolating factors */
 241         sx = s_curve(rx0);
 242         sy = s_curve(ry0);
 243
 244         /* interpolate along the left edge */
 245         u = grad2[b00][0] * rx0 + grad2[b00][1] * ry0;
 246         v = grad2[b10][0] * rx1 + grad2[b10][1] * ry0;
 247         a = LERP(u, v, sx);
 248
 249         /* interpolate along the right edge */
 250         u = grad2[b01][0] * rx0 + grad2[b01][1] * ry1;
 251         v = grad2[b11][0] * rx1 + grad2[b11][1] * ry1;
 252         b = LERP(u, v, sx);
 253
 254         /* interpolate between them */
 255         return LERP(a, b, sy);
 256 }
 257
 258 float pnoise3(float x, float y, float z, int per_x, int per_y, int per_z)
 259 {
 260         int i, j;
 261         int bx0, bx1, by0, by1, bz0, bz1;
 262         int b00, b10, b01, b11;
 263         float rx0, rx1, ry0, ry1, rz0, rz1;
 264         float sx, sy, sz;
 265         float u, v, a, b, c, d;
 266
 267         init_once();
 268
 269         setup_p(x, bx0, bx1, rx0, rx1, per_x);
 270         setup_p(y, by0, by1, ry0, ry1, per_y);
 271         setup_p(z, bz0, bz1, rz0, rz1, per_z);
 272
 273         i = perm[bx0];
 274         j = perm[bx1];
 275
 276         b00 = perm[i + by0];
 277         b10 = perm[j + by0];
 278         b01 = perm[i + by1];
 279         b11 = perm[j + by1];
 280
 281         /* calculate hermite interpolating factors */
 282         sx = s_curve(rx0);
 283         sy = s_curve(ry0);
 284         sz = s_curve(rz0);
 285
 286         /* interpolate along the top slice of the cell */
 287         u = grad3[b00 + bz0][0] * rx0 + grad3[b00 + bz0][1] * ry0 + grad3[b00 + bz0][2] * rz0;
 288         v = grad3[b10 + bz0][0] * rx1 + grad3[b10 + bz0][1] * ry0 + grad3[b10 + bz0][2] * rz0;
 289         a = LERP(u, v, sx);
 290
 291         u = grad3[b01 + bz0][0] * rx0 + grad3[b01 + bz0][1] * ry1 + grad3[b01 + bz0][2] * rz0;
 292         v = grad3[b11 + bz0][0] * rx1 + grad3[b11 + bz0][1] * ry1 + grad3[b11 + bz0][2] * rz0;
 293         b = LERP(u, v, sx);
 294
 295         c = LERP(a, b, sy);
 296
 297         /* interpolate along the bottom slice of the cell */
 298         u = grad3[b00 + bz0][0] * rx0 + grad3[b00 + bz0][1] * ry0 + grad3[b00 + bz0][2] * rz1;
 299         v = grad3[b10 + bz0][0] * rx1 + grad3[b10 + bz0][1] * ry0 + grad3[b10 + bz0][2] * rz1;
 300         a = LERP(u, v, sx);
 301
 302         u = grad3[b01 + bz0][0] * rx0 + grad3[b01 + bz0][1] * ry1 + grad3[b01 + bz0][2] * rz1;
 303         v = grad3[b11 + bz0][0] * rx1 + grad3[b11 + bz0][1] * ry1 + grad3[b11 + bz0][2] * rz1;
 304         b = LERP(u, v, sx);
 305
 306         d = LERP(a, b, sy);
 307
 308         /* interpolate between slices */
 309         return LERP(c, d, sz);
 310 }
 311
 312
 313 float fbm1(float x, int octaves)
 314 {
 315         int i;
 316         float res = 0.0f, freq = 1.0f;
 317         for(i=0; i<octaves; i++) {
 318                 res += noise1(x * freq) / freq;
 319                 freq *= 2.0f;
 320         }
 321         return res;
 322 }
 323
 324 float fbm2(float x, float y, int octaves)
 325 {
 326         int i;
 327         float res = 0.0f, freq = 1.0f;
 328         for(i=0; i<octaves; i++) {
 329                 res += noise2(x * freq, y * freq) / freq;
 330                 freq *= 2.0f;
 331         }
 332         return res;
 333 }
 334
 335 float fbm3(float x, float y, float z, int octaves)
 336 {
 337         int i;
 338         float res = 0.0f, freq = 1.0f;
 339         for(i=0; i<octaves; i++) {
 340                 res += noise3(x * freq, y * freq, z * freq) / freq;
 341                 freq *= 2.0f;
 342         }
 343         return res;
 344
 345 }
 346
 347
 348 float pfbm1(float x, int per, int octaves)
 349 {
 350         int i;
 351         float res = 0.0f, freq = 1.0f;
 352         for(i=0; i<octaves; i++) {
 353                 res += pnoise1(x * freq, per) / freq;
 354                 freq *= 2.0f;
 355                 per *= 2;
 356         }
 357         return res;
 358 }
 359
 360 float pfbm2(float x, float y, int per_x, int per_y, int octaves)
 361 {
 362         int i;
 363         float res = 0.0f, freq = 1.0f;
 364         for(i=0; i<octaves; i++) {
 365                 res += pnoise2(x * freq, y * freq, per_x, per_y) / freq;
 366                 freq *= 2.0f;
 367                 per_x *= 2;
 368                 per_y *= 2;
 369         }
 370         return res;
 371 }
 372
 373 float pfbm3(float x, float y, float z, int per_x, int per_y, int per_z, int octaves)
 374 {
 375         int i;
 376         float res = 0.0f, freq = 1.0f;
 377         for(int i=0; i<octaves; i++) {
 378                 res += pnoise3(x * freq, y * freq, z * freq, per_x, per_y, per_z) / freq;
 379                 freq *= 2.0f;
 380                 per_x *= 2;
 381                 per_y *= 2;
 382                 per_z *= 2;
 383         }
 384         return res;
 385 }
 386
 387
 388 float turbulence1(float x, int octaves)
 389 {
 390         int i;
 391         float res = 0.0f, freq = 1.0f;
 392         for(i=0; i<octaves; i++) {
 393                 res += fabs(noise1(x * freq) / freq);
 394                 freq *= 2.0f;
 395         }
 396         return res;
 397 }
 398
 399 float turbulence2(float x, float y, int octaves)
 400 {
 401         int i;
 402         float res = 0.0f, freq = 1.0f;
 403         for(i=0; i<octaves; i++) {
 404                 res += fabs(noise2(x * freq, y * freq) / freq);
 405                 freq *= 2.0f;
 406         }
 407         return res;
 408 }
 409
 410 float turbulence3(float x, float y, float z, int octaves)
 411 {
 412         int i;
 413         float res = 0.0f, freq = 1.0f;
 414         for(i=0; i<octaves; i++) {
 415                 res += fabs(noise3(x * freq, y * freq, z * freq) / freq);
 416                 freq *= 2.0f;
 417         }
 418         return res;
 419 }
 420
 421
 422 float pturbulence1(float x, int per, int octaves)
 423 {
 424         int i;
 425         float res = 0.0f, freq = 1.0f;
 426         for(i=0; i<octaves; i++) {
 427                 res += fabs(pnoise1(x * freq, per) / freq);
 428                 freq *= 2.0f;
 429                 per *= 2;
 430         }
 431         return res;
 432 }
 433
 434 float pturbulence2(float x, float y, int per_x, int per_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(pnoise2(x * freq, y * freq, per_x, per_y) / freq);
 440                 freq *= 2.0f;
 441                 per_x *= 2;
 442                 per_y *= 2;
 443         }
 444         return res;
 445 }
 446
 447 float pturbulence3(float x, float y, float z, int per_x, int per_y, int per_z, int octaves)
 448 {
 449         int i;
 450         float res = 0.0f, freq = 1.0f;
 451         for(i=0; i<octaves; i++) {
 452                 res += fabs(pnoise3(x * freq, y * freq, z * freq, per_x, per_y, per_z) / freq);
 453                 freq *= 2.0f;
 454                 per_x *= 2;
 455                 per_y *= 2;
 456                 per_z *= 2;
 457         }
 458         return res;
 459 }