dc-terrain / libdc / include / dc_noise.h

#ifndef DC_NOISE_H
#define DC_NOISE_H

// Open Simplex 2S noise (3D), public domain.
// Based on K.jpg's OpenSimplex2S implementation, simplified for our use case.
// Provides seeded 3D simplex-style noise with fractal octaves.

#include <cmath>
#include <cstdint>

class SimplexNoise {
public:
    SimplexNoise() { set_seed(0); }
    explicit SimplexNoise(int seed) { set_seed(seed); }

    void set_seed(int seed) {
        // Build permutation table from seed using LCG
        uint32_t s = static_cast<uint32_t>(seed);
        for (int i = 0; i < 256; i++) perm_[i] = i;
        for (int i = 255; i > 0; i--) {
            s = s * 1664525u + 1013904223u;
            int j = static_cast<int>((s >> 16) % (i + 1));
            int tmp = perm_[i];
            perm_[i] = perm_[j];
            perm_[j] = tmp;
        }
        for (int i = 0; i < 256; i++) {
            perm_[256 + i] = perm_[i];
        }
    }

    // Single-evaluation 3D noise, returns roughly [-1, 1]
    float noise3(float x, float y, float z) const {
        // Classic simplex noise 3D
        // Skew input space
        const float F3 = 1.0f / 3.0f;
        const float G3 = 1.0f / 6.0f;

        float s = (x + y + z) * F3;
        int i = fastfloor(x + s);
        int j = fastfloor(y + s);
        int k = fastfloor(z + s);

        float t = (i + j + k) * G3;
        float X0 = i - t;
        float Y0 = j - t;
        float Z0 = k - t;
        float x0 = x - X0;
        float y0 = y - Y0;
        float z0 = z - Z0;

        int i1, j1, k1;
        int i2, j2, k2;

        if (x0 >= y0) {
            if (y0 >= z0) {
                i1=1; j1=0; k1=0; i2=1; j2=1; k2=0;
            } else if (x0 >= z0) {
                i1=1; j1=0; k1=0; i2=1; j2=0; k2=1;
            } else {
                i1=0; j1=0; k1=1; i2=1; j2=0; k2=1;
            }
        } else {
            if (y0 < z0) {
                i1=0; j1=0; k1=1; i2=0; j2=1; k2=1;
            } else if (x0 < z0) {
                i1=0; j1=1; k1=0; i2=0; j2=1; k2=1;
            } else {
                i1=0; j1=1; k1=0; i2=1; j2=1; k2=0;
            }
        }

        float x1 = x0 - i1 + G3;
        float y1 = y0 - j1 + G3;
        float z1 = z0 - k1 + G3;
        float x2 = x0 - i2 + 2.0f * G3;
        float y2 = y0 - j2 + 2.0f * G3;
        float z2 = z0 - k2 + 2.0f * G3;
        float x3 = x0 - 1.0f + 3.0f * G3;
        float y3 = y0 - 1.0f + 3.0f * G3;
        float z3 = z0 - 1.0f + 3.0f * G3;

        int ii = i & 255;
        int jj = j & 255;
        int kk = k & 255;

        float n0 = contrib(ii, jj, kk, x0, y0, z0);
        float n1 = contrib(ii + i1, jj + j1, kk + k1, x1, y1, z1);
        float n2 = contrib(ii + i2, jj + j2, kk + k2, x2, y2, z2);
        float n3 = contrib(ii + 1, jj + 1, kk + 1, x3, y3, z3);

        return 32.0f * (n0 + n1 + n2 + n3);
    }

    // Fractal Brownian Motion: multiple octaves
    float fbm3(float x, float y, float z, int octaves,
               float lacunarity = 2.0f, float persistence = 0.5f) const {
        float sum = 0.0f;
        float amp = 1.0f;
        float freq = 1.0f;
        float max_amp = 0.0f;
        for (int o = 0; o < octaves; o++) {
            sum += noise3(x * freq, y * freq, z * freq) * amp;
            max_amp += amp;
            amp *= persistence;
            freq *= lacunarity;
        }
        return sum / max_amp; // normalize to ~[-1, 1]
    }

    // Public accessor for GPU upload
    const int* perm_data() const { return perm_; }

private:
    int perm_[512];

    static int fastfloor(float x) {
        int xi = static_cast<int>(x);
        return (x < xi) ? xi - 1 : xi;
    }

    // 12 gradient directions for 3D simplex
    static float grad3(int hash, float x, float y, float z) {
        int h = hash & 15;
        float u = h < 8 ? x : y;
        float v = h < 4 ? y : (h == 12 || h == 14 ? x : z);
        return ((h & 1) ? -u : u) + ((h & 2) ? -v : v);
    }

    float contrib(int ii, int jj, int kk, float x, float y, float z) const {
        float t = 0.6f - x * x - y * y - z * z;
        if (t < 0.0f) return 0.0f;
        t *= t;
        int gi = perm_[(ii + perm_[(jj + perm_[kk & 255]) & 255]) & 255];
        return t * t * grad3(gi, x, y, z);
    }
};

#endif // DC_NOISE_H