Download
541 lines 22.3 KiBLuaModified 2017-10-28 raw
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541
-- Hill Noise library
-- Copyright 2017 Bruce Hill
--
-- This library provides functions for continuous pseudorandom 1D, 2D, and 3D noise functions
-- with approximately uniform distribution on the interval (0,1) and customizable resolution,
-- levels of detail, and shape characteristics.
--
-- Functions: make1d([opts]), make2d([opts]), make3d([opts])
--   each of these returns a noise function, and a gradient function of the noise function
-- Shader functions: make1dShader([opts]), make2dShader([opts]), make3dShader([opts])
--   each of these returns a noise shader, and a gradient shader of the noise function
--
-- Usage:
--     local Noise = require 'hill_noise'
--     local noise,noise_gradient = Noise.make3d{resolution=9}
--     local x,y,z = 1,2,3
--     local n = noise(x,y,z)
--     local dndx,dndy,dndz = noise_gradient(x,y,z)
--
-- The noise construction functions all take the following options:
--     * random: a random number function (default: math.random)
--     * seed: if "random" is not provided, and the love.math module is loaded, a new
--       pseudorandom number generator will be used with the specified seed
--     * amplitudes: a list of sine wave amplitudes to use
--     * resolution: if "amplitudes" is not provided, the number of sine waves to use
--       (default: 5,7,9 for 1d,2d,3d functions, and 10,15,20 for 1d,2d,3d shaders)
--     * distribution: if "amplitudes" is not provided, a function to evenly sample
--       amplitudes from on the interval (0,1) (default: 0.1^x)
--
-- The noise and noise gradient functions all run in O(resolution), and more resolution
-- may be needed for higher dimensional noise and additional detail.
--
-- For the LOVE game engine, there are versions that produce 1D, 2D, and 3D noise shaders
-- as well, which are much faster and run on the GPU.
--     local noise_shader, noise_gradient_shader = Noise.make3dShader{resolution=11}
--     local canvas = love.graphics.newCanvas()
--     -- At each location on the texture, the noise function is sampled at a linear
--     --   interpolation between range_min and range_max, using the texture coordinates.
--     --   For 3D noise, the "z" value is passed as a parameter.
--     noise_shader:send("range_min", {0,20})
--     noise_shader:send("range_max", {0,20})
--     noise_shader:send('z', love.timer.getTime())
--     love.graphics.setShader(noise_shader)
--     love.graphics.draw(canvas)
--     love.graphics.setShader()
-- The noise gradient shaders populate the RGB channels with:
--   * for 1D: sigmoid(2*dn/dx), sigmoid(2*dn/dx), sigmoid(2*dn/dx)
--   * for 2D: sigmoid(1.5*dn/dx), sigmoid(1.5*dn/dy), 0
--   * for 3D: sigmoid(4*dn/dx), sigmoid(4*dn/dy), sigmoid(4*dn/dz)
-- (coefficients were arbitrarily chosen to best squash typical gradients for the default
--  parameters evenly into (0,1))
-- Additionally, the 1D noise shader has an extra option "draw_outline" that, if set to
-- true, makes the 1D shader render output as white or black if the y-coordinate is
-- above/below the  noise value for a given x-coordinate.
-- 
-- Permission is hereby granted, free of charge, to any person obtaining a copy of
-- this software and associated documentation files (the "Software"), to deal in
-- the Software without restriction, including without limitation the rights to
-- use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
-- of the Software, and to permit persons to whom the Software is furnished to do
-- so, subject to the following conditions:
-- 
-- The above copyright notice and this permission notice shall be included in all
-- copies or substantial portions of the Software.
-- 
-- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-- IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-- FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-- AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-- LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-- OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-- SOFTWARE.

local noise = {}
local exp,sin,cos,floor,log,acos,sqrt,abs = math.exp,math.sin,math.cos,math.floor,math.log,math.acos,math.sqrt,math.abs
local GR, PI, TAU, SQRT5, LOG_GR = (sqrt(5)+1)/2, math.pi, 2*math.pi, sqrt(5), log((sqrt(5)+1)/2)

local function cdf(x)
    return .5 + .5*(x<0 and -1 or 1)*sqrt(1.-exp(-4./TAU * x*x))
end

local function cdf_prime(x, dx)
    return (0.31831 * exp(-2/PI * x*x) * abs(x)*dx)/sqrt(1-exp(-2/PI*x*x))
end

local default_resolutions = {["1d"]=5, ["2d"]=7, ["3d"]=9, ["1dShader"]=10, ["2dShader"]=15, ["3dShader"]=20}
local function _defaultArgs(opts, kind)
    opts = opts or {}
    local amplitudes, random = opts.amplitudes, opts.random
    if random then
        assert(type(random) == 'function', "Random number function must be a function.")
    elseif not random then
        if opts.seed and love and love.math then
            local rng = love.math.newRandomGenerator(opts.seed)
            random = function(...) return rng:random(...) end
        else
            random = math.random
        end
    end
    if not amplitudes then
        local resolution = opts.resolution or default_resolutions[kind]
        local distribution = opts.distribution or (function(x) return .1^x end)
        amplitudes = {}
        for i=1,resolution do
            amplitudes[i] = distribution((i-.5)/resolution)
        end
    end
    return amplitudes, random
end

local function calculateSigma(amplitudes)
    local sigma = 0
    for _,a in ipairs(amplitudes) do
        sigma = sigma + a*a
    end
    return math.sqrt(sigma/2)
end

local function makeOffsets(count, random)
    local offsets = {}
    for i=1,count do
        offsets[i] = random()*TAU
    end
    return offsets
end

noise.make1d = function(opts)
    local amplitudes, random = _defaultArgs(opts, "1d")
    local sigma = calculateSigma(amplitudes)
    local offsets = makeOffsets(#amplitudes, random)
    local function noise(x)
        local n = 0
        for i,a in ipairs(amplitudes) do
            n = n + a*cos(x/a + offsets[i])
        end
        n = n/sigma
        return cdf(n)
    end

    local function gradient(x)
        local n, dndx = 0, 0
        for i,a in ipairs(amplitudes) do
            n = n + a*cos(x/a + offsets[i])
            dndx = dndx - sin(x/a + offsets[i])
        end
        dndx = dndx/sigma
        n = n/sigma
        return cdf_prime(n,dndx)
    end
    return noise, gradient
end

noise.make2d = function(opts)
    local amplitudes, random = _defaultArgs(opts, "2d")
    local sigma = calculateSigma(amplitudes)
    local offsets = makeOffsets(2*#amplitudes, random)
    sigma = sigma/sqrt(2)
    local function noise(x,y)
        local n = 0
        for i,a in ipairs(amplitudes) do
            local angle = ((i*GR) % 1)*TAU
            local sina, cosa = sin(angle), cos(angle)
            local u = x*cosa - y*sina
            local v = -x*sina - y*cosa
            local k, w = offsets[2*i], offsets[2*i-1]
            n = n + a*(cos(u/a + k) + cos(v/a + w))
        end
        n = n/(2*sigma)
        return cdf(n)
    end
    local function gradient(x,y)
        local n, dx, dy = 0,0,0
        for i,a in ipairs(amplitudes) do
            local angle = ((i*GR) % 1)*TAU
            local sina, cosa = sin(angle), cos(angle)
            local u = x*cosa - y*sina
            local v = -x*sina - y*cosa
            local k, w = offsets[2*i], offsets[2*i-1]
            n = n + a*(cos(u/a + k) + cos(v/a + w))
            dx = dx + (sina*sin(v/a + w) - cosa*sin(u/a + k))
            dy = dy + (sina*sin(u/a + k) + cosa*sin(v/a + w))
        end
        n = n/(2*sigma)
        dx = dx/(2*sigma)
        dy = dy/(2*sigma)
        return cdf_prime(n,dx), cdf_prime(n,dy)
    end
    return noise, gradient
end

noise.make3d = function(opts)
    local amplitudes, random = _defaultArgs(opts, "3d")
    local resolution = #amplitudes
    local sigma = calculateSigma(amplitudes)
    local offsets = makeOffsets(3*#amplitudes, random)
    sigma = sigma/sqrt(3)
    local function noise(x,y,z)
        -- Find the biggest fibonacci number F_n such that F_n < resolution
        local n = floor(log((resolution-1)*SQRT5 + .5)/LOG_GR)
        local dec = floor(.5 + (GR^n)/SQRT5) -- F_n, using closed form Fibonacci
        local inc = floor(.5 + dec/GR) -- F_(n-1)

        local n,i,j = 0,0,0
        for i=0,resolution-1 do
            if j >= dec then
                j = j - dec
            else
                j = j + inc
                if j >= resolution then
                    j = j - dec
                end
            end
            -- Convert golden ratio sequence into polar coordinate unit vector
            local phi = ((i*GR) % 1)*TAU
            local theta = acos(-1+2*((j*GR) % 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
            local n1x,n1y,n1z = sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)
            local n2x,n2y,n2z = sin(phi+TAU/4.)*cos(theta), sin(phi+TAU/4.)*sin(theta), cos(phi+TAU/4.)
            -- Cross product
            local n3x,n3y,n3z = n1y*n2z - n1z*n2y,
                        n1z*n2x - n1x*n2z,
                        n1x*n2y - n1y*n2x
            -- Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
            local u = n1x*x + n1y*y + n1z*z
            local v = n2x*x + n2y*y + n2z*z
            local w = n3x*x + n3y*y + n3z*z
            -- Pull the amplitude from the shuffled array index ("j"), not "i",
            -- otherwise neighboring unit vectors will have similar amplitudes!
            local a = amplitudes[j+1]
            -- Noise is the average of cosine of distance along each axis, shifted by offsets and scaled by amplitude.
            n = n + a*(cos(u/a + offsets[3*i+1]) + cos(v/a + offsets[3*i+2]) + cos(w/a + offsets[3*i+3]))
        end
        return cdf(n/(3*sigma))
    end
    local function gradient(x,y,z)
        -- Find the biggest fibonacci number F_n such that F_n < resolution
        local n = floor(log((resolution-1)*SQRT5 + .5)/LOG_GR)
        local dec = floor(.5 + (GR^n)/SQRT5) -- F_n, using closed form Fibonacci
        local inc = floor(.5 + dec/GR) -- F_(n-1)

        local n,i,j = 0,0,0
        local dndx,dndy,dndz = 0,0,0
        for i=0,resolution-1 do
            if j >= dec then
                j = j - dec
            else
                j = j + inc
                if j >= resolution then
                    j = j - dec
                end
            end
            -- Convert golden ratio sequence into polar coordinate unit vector
            local phi = ((i*GR) % 1)*TAU
            local theta = acos(-1+2*((j*GR) % 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
            local n1x,n1y,n1z = sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)
            local n2x,n2y,n2z = sin(phi+TAU/4.)*cos(theta), sin(phi+TAU/4.)*sin(theta), cos(phi+TAU/4.)
            -- Cross product
            local n3x,n3y,n3z = n1y*n2z - n1z*n2y,
                        n1z*n2x - n1x*n2z,
                        n1x*n2y - n1y*n2x
            -- Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
            local u = n1x*x + n1y*y + n1z*z
            local v = n2x*x + n2y*y + n2z*z
            local w = n3x*x + n3y*y + n3z*z
            -- Pull the amplitude from the shuffled array index ("j"), not "i",
            -- otherwise neighboring unit vectors will have similar amplitudes!
            local a = amplitudes[j+1]
            -- Noise is the average of cosine of distance along each axis, shifted by offsets and scaled by amplitude.
            n = n + a*(cos(u/a + offsets[3*i+1]) + cos(v/a + offsets[3*i+2]) + cos(w/a + offsets[3*i+3]))

            local kx,ky,kz = -sin(u/a+offsets[3*i+1]),-sin(v/a+offsets[3*i+2]),-sin(w/a + offsets[3*i+3])
            dndx = dndx + (n1x*kx + n2x*ky + n3x*kz)
            dndy = dndy + (n1y*kx + n2y*ky + n3y*kz)
            dndz = dndz + (n1z*kx + n2z*ky + n3z*kz)
        end
        n = n / (3*sigma)
        return cdf_prime(n, dndx/(3*sigma)), cdf_prime(n, dndy/(3*sigma)), cdf_prime(n, dndz/(3*sigma))
    end
    return noise, gradient
end

if love and love.graphics then
    local shader1d = [[
    #define TAU 6.283185307179586476925286766559005768394338798750211641949
    #define MAX_RESOLUTION 64
    #define SIGMOID(x) (1./(1.+exp(-(x))))
    #define CDF(x) (.5 + .5*sign(x)*sqrt(1.-exp(-4./TAU * (x)*(x))))
    #define CDF_PRIME(x, dx) ((0.31831 * exp(-4./TAU * (x)*(x)) * abs(x)*(dx))/sqrt(1.0-exp(-4./TAU * (x)*(x))))
    extern int resolution;
    extern float sigma, range_min, range_max;
    extern float amplitudes[MAX_RESOLUTION];
    extern float offsets[MAX_RESOLUTION];
    extern bool draw_outline = false;

    vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
    {
        float x = mix(range_min,range_max,texture_coords.x);
        float noise = 0.;
    #ifdef GRADIENT
        float dndx = 0.;
    #endif
        for (int i=0; i < resolution; i++) {
            float a = amplitudes[i];
            noise += a*cos(x/a + offsets[i]);
    #ifdef GRADIENT
            dndx -= sin(x/a + offsets[i]);
    #endif
        }
        noise /= sigma;
    #ifdef GRADIENT
        dndx = CDF_PRIME(noise, dndx/sigma);
        dndx = SIGMOID(dndx*2.);
        if (draw_outline) {
            return texture_coords.y > (1.-dndx) ? vec4(1.,1.,1.,1.) : vec4(0.,0.,0.,1.);
        } else {
            return vec4(dndx,dndx,dndx,1.);
        }
    #else
        noise = CDF(noise);
        if (draw_outline) {
            return texture_coords.y > (1.-noise) ? vec4(1.,1.,1.,1.) : vec4(0.,0.,0.,1.);
        } else {
            return vec4(noise,noise,noise,1.);
        }
    #endif
    }
    ]]

    noise.make1dShader = function(opts)
        local amplitudes, random = _defaultArgs(opts, "1dShader")
        local resolution = #amplitudes
        local sigma = calculateSigma(amplitudes)
        local offsets = makeOffsets(#amplitudes, random)
        local shaderCode = "#define RESOLUTION "..tostring(resolution).."\n"..shader1d
        local noiseShader = love.graphics.newShader(shaderCode)
        local gradShader = love.graphics.newShader("#define GRADIENT\n"..shaderCode)
        do -- Dumb hack to work around a bug in Love 0.10.2 not sending the last value
            table.insert(amplitudes, 1.)
            table.insert(offsets, 1.)
        end
        for _,shader in ipairs{noiseShader, gradShader} do
            shader:send("sigma",sigma)
            shader:send("resolution",resolution)
            shader:send("offsets",unpack(offsets))
            shader:send("amplitudes",unpack(amplitudes))
            shader:send("range_min", 0)
            shader:send("range_max", 0)
            if opts and opts.draw_outline then
                shader:send("draw_outline", true)
            end
        end
        return noiseShader,gradShader
    end

    local shader2d = [[
    #define TAU 6.283185307179586476925286766559005768394338798750211641949
    #define PHI 1.618033988749894848204586834365638117720309179805762862135
    #define CDF(x) (.5 + .5*sign(x)*sqrt(1.-exp(-4./TAU * (x)*(x))))
    #define CDF_PRIME(x, dx) ((0.31831 * exp(-4./TAU * (x)*(x)) * abs(x)*(dx))/sqrt(1.0-exp(-4./TAU * (x)*(x))))
    #define SIGMOID(x) (1./(1.+exp(-2.*(x))))
    extern float sigma;
    extern float amplitudes[RESOLUTION];
    extern vec2 offsets[RESOLUTION];
    extern vec2 range_min, range_max;

    vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
    {
        vec2 pos = mix(range_min,range_max,texture_coords);
        float noise = 0.;
    #ifdef GRADIENT
        float dndx = 0., dndy = 0.;
    #endif
        for (int i=0; i < RESOLUTION; i++) {
            float angle = mod(float(i)*PHI, 1.)*TAU;
            float cosa = cos(angle), sina = sin(angle);
            float u = pos.x*cosa - pos.y*sina;
            float v = -pos.x*sina - pos.y*cosa;
            float a = amplitudes[i];
            float k = offsets[i].x, w = offsets[i].y;
            noise += a*(cos(u/a + k) + cos(v/a + w));
    #ifdef GRADIENT
            dndx += -cosa*sin(u/a + k) + sina*sin(v/a + w);
            dndy += sina*sin(u/a + k) + cosa*sin(v/a + w);
    #endif
        }

        float _sigma = 2.*sigma;
        noise /= _sigma;
    #ifdef GRADIENT
        dndx = CDF_PRIME(noise, dndx/_sigma);
        dndx = SIGMOID(dndx*1.5);

        dndy = CDF_PRIME(noise, dndy/_sigma);
        dndy = SIGMOID(dndy*1.5);

        return vec4(dndx,dndy,0.,1.);
    #else
        noise = CDF(noise);
        return vec4(noise,noise,noise,1.);
    #endif
    }
    ]]

    noise.make2dShader = function(opts)
        local amplitudes, random = _defaultArgs(opts, "2dShader")
        local resolution = #amplitudes
        local sigma = calculateSigma(amplitudes)
        local offsets = makeOffsets(2*#amplitudes, random)
        local shaderCode = "#define RESOLUTION "..tostring(resolution).."\n"..shader2d
        local noiseShader = love.graphics.newShader(shaderCode)
        local gradShader = love.graphics.newShader("#define GRADIENT\n"..shaderCode)
        sigma = sigma/sqrt(2)
        local offsets2 = {}
        for i=1,#offsets-1,2 do
            table.insert(offsets2, {offsets[i],offsets[i+1]})
        end
        do -- Dumb hack to work around a bug in Love 0.10.2 not sending the last value
            table.insert(amplitudes, 1.)
            table.insert(offsets2, {1,1})
        end
        for _,shader in ipairs{noiseShader, gradShader} do
            shader:send("sigma",sigma)
            shader:send("offsets",unpack(offsets2))
            shader:send("amplitudes",unpack(amplitudes))
            shader:send("range_min", {0,0})
            shader:send("range_max", {1,1})
        end
        return noiseShader,gradShader
    end

    local shader3d = [[
    #define TAU 6.283185307179586476925286766559005768394338798750211641949
    #define PHI 1.618033988749894848204586834365638117720309179805762862135
    #define LOG_PHI 0.481211825059603447497758913424368423135184334385660519660
    #define SQRT5 2.236067977499789696409173668731276235440618359611525724270
    #define CDF(x) (.5 + .5*sign(x)*sqrt(1.-exp(-4./TAU * (x)*(x))))
    #define CDF_PRIME(x, dx) ((0.31831 * exp(-4./TAU * (x)*(x)) * abs(x)*(dx))/sqrt(1.0-exp(-4./TAU * (x)*(x))))
    #define SIGMOID(x) (1./(1.+exp(-2.*(x))))
    extern float sigma, z;
    extern float amplitudes[RESOLUTION];
    extern vec3 offsets[RESOLUTION];
    extern vec2 range_min, range_max;

    // https://www.graphics.rwth-aachen.de/media/papers/jgt.pdf
    vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
    {
        vec3 pos = vec3(mix(range_min,range_max,texture_coords), z);
        // Find the biggest fibonacci number F_n such that F_n < RESOLUTION
        int fib_n = int(log((float(RESOLUTION)-1.)*SQRT5 + .5)/LOG_PHI);
        int dec = int(.5 + pow(PHI,fib_n)/SQRT5); // F_n, using closed form Fibonacci
        int inc = int(.5 + dec/PHI); // F_(fib_n-1)

        float n = 0.;
    #ifdef GRADIENT
        float dndx = 0., dndy = 0., dndz = 0.;
    #endif
        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.
            n += a*(cos(u/a + offsets[i].x) + cos(v/a + offsets[i].y) + cos(w/a + offsets[i].z));
    #ifdef GRADIENT
            vec3 k = vec3(-sin(u/a+offsets[i].x),-sin(v/a+offsets[i].y),-sin(w/a + offsets[i].z));
            dndx += (n1.x*k.x + n2.x*k.y + n3.x*k.z)/3.;
            dndy += (n1.y*k.x + n2.y*k.y + n3.y*k.z)/3.;
            dndz += (n1.z*k.x + n2.z*k.y + n3.z*k.z)/3.;
    #endif
        }
        float _sigma = 3.*sigma;
        n /= _sigma;
    #ifdef GRADIENT
        dndx = CDF_PRIME(n, dndx/_sigma);
        dndx = SIGMOID(dndx*4.);

        dndy = CDF_PRIME(n, dndy/_sigma);
        dndy = SIGMOID(dndy*4.);

        dndz = CDF_PRIME(n, dndz);
        dndz = SIGMOID(dndz*4.);

        return vec4(dndx,dndy,dndz, 1.);
    #else
        n = CDF(n);
        return vec4(n,n,n,1.);
    #endif
    }
    ]]

    noise.make3dShader = function(opts)
        local amplitudes, random = _defaultArgs(opts, "3dShader")
        local resolution = #amplitudes
        local sigma = calculateSigma(amplitudes)
        local offsets = makeOffsets(3*#amplitudes, random)
        local shaderCode = "#define RESOLUTION "..tostring(resolution).."\n"..shader3d
        local noiseShader = love.graphics.newShader(shaderCode)
        local gradShader = love.graphics.newShader("#define GRADIENT\n"..shaderCode)
        sigma = sigma/sqrt(3)
        local offsets2 = {}
        for i=1,#offsets-1,3 do
            table.insert(offsets2, {offsets[i],offsets[i+1],offsets[i+2]})
        end
        do -- Dumb hack to work around a bug in Love 0.10.2 not sending the last value
            table.insert(amplitudes, 1.)
            table.insert(offsets, {1,1,1})
        end
        for _,shader in ipairs{noiseShader, gradShader} do
            shader:send("sigma",sigma)
            shader:send("offsets",unpack(offsets2))
            shader:send("amplitudes",unpack(amplitudes))
            shader:send("range_min", {0,0})
            shader:send("range_max", {1,1})
        end
        return noiseShader, gradShader
    end
end

return noise