commit 77982a0ace6365ff15dab1c51ed2f4e87d6c6e90
Author: Bruce Hill <bitbucket@bruce-hill.com>
Date:   Fri Mar 31 01:23:09 2017 -0700

    Initial commit.

diff --git a/conf.lua b/conf.lua
new file mode 100644
index 0000000..2e721d7
--- /dev/null
+++ b/conf.lua
@@ -0,0 +1,21 @@
+function love.conf(t)
+    t.version = "0.10.2"
+    t.window.title = "Noise Demo"
+    t.window.width = 1024
+    t.window.height = 768
+    t.window.msaa = 2
+    t.window.highdpi = true
+
+    t.accelerometerjoystick = false
+    t.externalstorage = false
+
+    t.modules.audio = false
+    t.modules.image = false
+    t.modules.joystick = false
+    t.modules.physics = false
+    t.modules.sound = false
+    t.modules.system = false
+    t.modules.touch = false
+    t.modules.video = false
+    t.modules.thread = false
+end
diff --git a/main.lua b/main.lua
new file mode 100644
index 0000000..bc93e65
--- /dev/null
+++ b/main.lua
@@ -0,0 +1,124 @@
+lg = love.graphics
+W,H = lg.getDimensions()
+local Noise = require "noise"
+local DRAW_RES = love.window.toPixels(4)
+local SCALE = 100
+local font = lg.newFont(love.window.toPixels(24))
+
+local decay = function(x) return .1^x end
+local function rng()
+    local r = love.math.newRandomGenerator(1)
+    return function() return r:random() end
+end
+local n1 = Noise.make1d(5,rng(),decay)
+local n2 = Noise.make2d(7,rng(),decay)
+local n3 = Noise.make3d(11,rng(),decay)
+local s1 = Noise.make1dShader(7,rng(),decay)
+local s2 = Noise.make2dShader(13,rng(),decay)
+local s3 = Noise.make3dShader(17,rng(),decay)
+
+local cw,ch = W/2,H/3
+local canv = lg.newCanvas(cw,ch)
+
+function love.load()
+    xOffset = 0
+    yOffset = 0
+    t0 = love.timer.getTime()
+    love.mouse.setCursor(love.mouse.getSystemCursor("hand"))
+end
+
+function love.mousepressed()
+    love.mouse.setCursor(love.mouse.getSystemCursor("sizeall"))
+end
+
+function love.mousereleased()
+    love.mouse.setCursor(love.mouse.getSystemCursor("hand"))
+end
+
+function love.draw()
+    lg.setColor(255,255,255)
+    lg.setLineWidth(love.window.toPixels(2))
+    lg.setLineJoin("none")
+    local p1 = {}
+    for x=0,W-H/3,DRAW_RES do
+        table.insert(p1,x)
+        table.insert(p1,n1(x/SCALE+xOffset)*H/3)
+    end
+    lg.line(p1)
+
+    local p2 = {}
+    local t = love.timer.getTime()-t0
+    for x=0,W-H/3,DRAW_RES do
+        table.insert(p2,x)
+        table.insert(p2,(1+n2(x/SCALE+xOffset,t))*H/3)
+    end
+    lg.line(p2)
+
+    for y=0,H/6*1.2,2*DRAW_RES do
+        local p3 = {}
+        for x=0,W-H/3,DRAW_RES do
+            table.insert(p3,x)
+            table.insert(p3,2/3*H+H/6*n3(x/SCALE+xOffset,y/SCALE+yOffset,t)+y)
+        end
+        lg.line(p3)
+    end
+    for x=-H/6,W-H/3,DRAW_RES*2 do
+        local p3 = {}
+        local slant = .5
+        for y=0,H/6*1.2,2*DRAW_RES do
+            local x2 = x + slant*y
+            if x2 < -2*DRAW_RES then goto continue end
+            table.insert(p3,x2)
+            table.insert(p3,2/3*H+H/6*n3(x2/SCALE+xOffset,y/SCALE+yOffset,t)+y)
+            if x2 > W/2 then break end
+            ::continue::
+        end
+        if #p3 > 2 then
+            lg.line(p3)
+        end
+    end
+
+    lg.setShader(s1)
+    s1:send("range_min", xOffset)
+    s1:send("range_max", xOffset+cw/SCALE)
+    lg.draw(canv,W-cw,0)
+    lg.setShader(s2)
+    s2:send("range_min", {xOffset,yOffset})
+    s2:send("range_max", {xOffset+cw/SCALE,yOffset+ch/SCALE})
+    lg.draw(canv,W-cw,ch)
+    lg.setShader(s3)
+    s3:send("range_min", {xOffset,yOffset})
+    s3:send("range_max", {xOffset+cw/SCALE,yOffset+ch/SCALE})
+    s3:send('z', t)
+    lg.draw(canv,W-cw,2*ch)
+    lg.setShader(nil)
+
+    lg.setColor(255,200,0)
+    lg.setFont(font)
+    lg.printf("1D_Noise(x)", 5,5, W)
+    lg.printf("2D_Noise(x,time)", 5,5+H/3, W)
+    lg.printf("3D_Noise(x,y,time)", 5,5+2*H/3, W)
+    lg.printf("1D_Noise_Shader(x)", 0,5, W-5, 'right')
+    lg.printf("2D_Noise_Shader(x,y)", 0,5+H/3, W-5, 'right')
+    lg.printf("3D_Noise_Shader(x,y,time)", 0,5+2*H/3, W-5, 'right')
+end
+
+function love.update(dt)
+    local key = love.keyboard.isDown
+    local dx = 400/SCALE*dt*((key('right') and 1 or 0) + (key('left') and -1 or 0))
+    local dy = 400/SCALE*dt*((key('down') and 1 or 0) + (key('up') and -1 or 0))
+    xOffset = xOffset + dx
+    yOffset = yOffset + dy
+end
+
+function love.keypressed(key)
+    if key == 'escape' then love.event.quit() end
+    if key == 'r' then love.load() end
+end
+
+function love.mousemoved(x,y,dx,dy)
+    if love.mouse.isDown(1) then
+        xOffset = xOffset - dx/SCALE
+        yOffset = yOffset - dy/SCALE
+    end
+end
diff --git a/noise.lua b/noise.lua
new file mode 100644
index 0000000..de528c6
--- /dev/null
+++ b/noise.lua
@@ -0,0 +1,150 @@
+local noise = {}
+local exp,sin,cos,floor,log,acos,sqrt = math.exp,math.sin,math.cos,math.floor,math.log,math.acos,math.sqrt
+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,sigma)
+    return .5 + .5*(x<0 and -1 or 1)*sqrt(1.-exp(-4./TAU * x*x))
+end
+
+local function _defaultArgs(resolution,random,decayFn)
+    if not resolution then resolution = 4 end
+    if not random then random = math.random end
+    if not decayFn then decayFn = function(x) return .1^x end end
+    return resolution,random,decayFn
+end
+
+local function _amplitudesAndOffsets(decayFn,numAmplitudes, numOffsets, random)
+    local sigma = 0
+    local amplitudes = {}
+    for i=1,numAmplitudes do
+        local a = decayFn((i-1+random())/numAmplitudes)
+        amplitudes[i] = a
+        sigma = sigma + a^2
+    end
+    sigma = math.sqrt(sigma/2)
+    local offsets = {}
+    for i=1,numOffsets do
+        offsets[i] = random()*TAU
+    end
+    return amplitudes,sigma,offsets
+end
+
+noise.make1d = function(resolution,random,decayFn)
+    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
+    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,resolution,random)
+    return function(x)
+        local noise = 0
+        for i,a in ipairs(amplitudes) do
+            noise = noise + a*cos(x/a + offsets[i])
+        end
+        return cdf(noise/sigma);
+    end
+end
+
+noise.make2d = function(resolution,random,decayFn)
+    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
+    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,2*resolution,random)
+    return function(x,y)
+        local noise = 0
+        for i,a in ipairs(amplitudes) do
+            local angle = ((i*GR) % 1)*TAU
+            local u = x*cos(angle) - y*sin(angle)
+            local v = x*cos(angle+TAU/4) - y*sin(angle+TAU/4)
+            noise = noise + a/2*(cos(u/a + offsets[2*i]) + cos(v/a + offsets[2*i-1]))
+        end
+        return cdf(noise/sigma);
+    end
+end
+
+noise.make3d = function(resolution,random,decayFn)
+    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
+    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,3*resolution,random)
+    return function(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 noise,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 n1 = {sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)}
+            local n2 = {sin(phi+TAU/4.)*cos(theta), sin(phi+TAU/4.)*sin(theta), cos(phi+TAU/4.)}
+            -- Cross product
+            local n3 = {n1[2]*n2[3] - n1[3]*n2[2],
+                        n1[3]*n2[1] - n1[1]*n2[3],
+                        n1[1]*n2[2] - n1[2]*n2[1]}
+            -- Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
+            local u = n1[1]*x + n1[2]*y + n1[3]*z
+            local v = n2[1]*x + n2[2]*y + n2[3]*z
+            local w = n3[1]*x + n3[2]*y + n3[3]*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.
+            noise = noise + a*(cos(u/a + offsets[3*i+1]) + cos(v/a + offsets[3*i+2]) + cos(w/a + offsets[3*i+3]))/3
+        end
+        return cdf(noise/sigma)
+    end
+end
+
+noise.make1dShader = function(resolution,random,decayFn)
+    local shader = lg.newShader("noise1d.glsl")
+    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
+    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,resolution,random)
+    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)
+    return shader
+end
+
+noise.make2dShader = function(resolution,random,decayFn)
+    local shader = lg.newShader("noise2d.glsl")
+    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
+    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,2*resolution,random)
+    local offsets2 = {}
+    for i=1,#offsets-1,2 do
+        table.insert(offsets2, {offsets[i],offsets[i+1]})
+    end
+    shader:send("sigma",sigma)
+    shader:send("resolution",resolution)
+    shader:send("offsets",unpack(offsets2))
+    shader:send("amplitudes",unpack(amplitudes))
+    shader:send("range_min", {0,0})
+    shader:send("range_max", {1,1})
+    return shader
+end
+
+noise.make3dShader = function(resolution,random,decayFn)
+    local shader = lg.newShader("noise3d.glsl")
+    resolution,random,decayFn = _defaultArgs(resolution,random,decayFn)
+    local amplitudes,sigma,offsets = _amplitudesAndOffsets(decayFn,resolution,3*resolution,random)
+    local offsets2 = {}
+    for i=1,#offsets-1,3 do
+        table.insert(offsets2, {offsets[i],offsets[i+1],offsets[i+2]})
+    end
+    shader:send("sigma",sigma)
+    shader:send("resolution",resolution)
+    shader:send("offsets",unpack(offsets2))
+    shader:send("amplitudes",unpack(amplitudes))
+    shader:send("range_min", {0,0})
+    shader:send("range_max", {1,1})
+    return shader
+end
+
+return noise
diff --git a/noise1d.glsl b/noise1d.glsl
new file mode 100644
index 0000000..1c1d40c
--- /dev/null
+++ b/noise1d.glsl
@@ -0,0 +1,28 @@
+#define TAU 6.283185307179586476925286766559005768394338798750211641949
+#define MAX_RESOLUTION 64
+extern int resolution;
+extern float sigma;
+extern float amplitudes[MAX_RESOLUTION];
+extern float offsets[MAX_RESOLUTION];
+extern float range_min, range_max;
+
+float cdf(float x) {
+    return .5 + .5*sign(x)*sqrt(1.-exp(-4./TAU * x*x));
+}
+
+float noise1d(float x) {
+    float noise = 0.;
+    for (int i=0; i < resolution; i++) {
+        float a = amplitudes[i];
+        noise += a*cos(x/a + offsets[i]);
+    }
+    return cdf(noise/sigma);
+}
+
+vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
+{
+    float x = mix(range_min,range_max,texture_coords.x);
+    float n = noise1d(x);
+    //return vec4(n,n,n,1.);
+    return texture_coords.y > n ? vec4(1.,1.,1.,1.) : vec4(0.,0.,0.,1.);
+}
diff --git a/noise2d.glsl b/noise2d.glsl
new file mode 100644
index 0000000..860c4b0
--- /dev/null
+++ b/noise2d.glsl
@@ -0,0 +1,32 @@
+#define TAU 6.283185307179586476925286766559005768394338798750211641949
+#define PHI 1.618033988749894848204586834365638117720309179805762862135
+#define MAX_RESOLUTION 64
+extern int resolution;
+extern float sigma;
+extern float amplitudes[MAX_RESOLUTION];
+extern vec2 offsets[MAX_RESOLUTION];
+extern vec2 range_min, range_max;
+
+float cdf(float x) {
+    return .5 + .5*sign(x)*sqrt(1.-exp(-4./TAU * x*x));
+}
+
+float noise2d(vec2 pos) {
+    float noise = 0.;
+    for (int i=0; i < resolution; i++) {
+        float angle = mod(float(i)*PHI, 1.)*TAU;
+        float u = pos.x*cos(angle) - pos.y*sin(angle);
+        float v = pos.x*cos(angle+TAU/4.) - pos.y*sin(angle+TAU/4.);
+        float a = amplitudes[i];
+        noise += a*mix(cos(u/a + offsets[i].x), cos(v/a + offsets[i].y), .5);
+    }
+    return cdf(noise/sigma);
+}
+
+
+vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 pixel_coords)
+{
+    vec2 coords = mix(range_min,range_max,texture_coords);
+    float n = noise2d(coords);
+    return vec4(n,n,n,1.);
+}
diff --git a/noise3d.glsl b/noise3d.glsl
new file mode 100644
index 0000000..900bf46
--- /dev/null
+++ b/noise3d.glsl
@@ -0,0 +1,60 @@
+#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.);
+}