#define TAU 6.283185307179586476925286766559005768394338798750211641949
#define PHI 1.618033988749894848204586834365638117720309179805762862135
#define LOG_PHI 0.481211825059603447497758913424368423135184334385660519660
#define SQRT5 2.236067977499789696409173668731276235440618359611525724270
#define MAX_RESOLUTION 64
extern int resolution;
extern float sigma;
extern float amplitudes[MAX_RESOLUTION];
extern vec3 offsets[MAX_RESOLUTION];
extern float z;
extern vec2 range_min, range_max;

float cdf(float x) {
    return .5 + .5*sign(x)*sqrt(1.-exp(-4./TAU * x*x));
}

// https://www.graphics.rwth-aachen.de/media/papers/jgt.pdf
float noise3d(vec3 pos) {
    // Find the biggest fibonacci number F_n such that F_n < RESOLUTION
    int n = int(log((float(resolution)-1.)*SQRT5 + .5)/LOG_PHI);
    int dec = int(.5 + pow(PHI,n)/SQRT5); // F_n, using closed form Fibonacci
    int inc = int(.5 + dec/PHI); // F_(n-1)

    float noise = 0.;
    for (int i=0, j=0; i<resolution; ++i) {
        if (j >= dec) {
            j -= dec;
        } else {
            j += inc;
            if (j >= resolution)
                j -= dec;
        }
        // Convert golden ratio sequence into polar coordinate unit vector
        float phi = mod(float(i)*PHI,1.)*TAU;
        float theta = acos(mix(-1.,1.,mod(float(j)*PHI,1.)));
        // Make an orthonormal basis, where n1 is from polar phi/theta,
        // n2 is roated 90 degrees along phi, and n3 is the cross product of the two
        vec3 n1 = vec3(sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi));
        vec3 n2 = vec3(sin(phi+TAU/4.)*cos(theta), sin(phi+TAU/4.)*sin(theta), cos(phi+TAU/4.));
        vec3 n3 = cross(n1,n2);
        // Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
        float u = dot(n1, pos);
        float v = dot(n2, pos);
        float w = dot(n3, pos);
        // Pull the amplitude from the shuffled array index ("j"), not "i",
        // otherwise neighboring unit vectors will have similar amplitudes!
        float a = amplitudes[j];
        //float a = pow(mod(float(i+1)*(PHI-1.), 1.), .3);
        // Noise is the average of cosine of distance along each axis, shifted by offsets and scaled by amplitude.
        noise += a*(cos(u/a + offsets[i].x) + cos(v/a + offsets[i].y) + cos(w/a + offsets[i].z))/3.;
    }
    return cdf(noise/sigma);
}

vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
{
    vec3 coords = vec3(mix(range_min,range_max,texture_coords), z);
    float n = noise3d(coords);
    return vec4(n,n,n,1.);
}