Mandelbulb Codea

--# Main
-- Mandlebulb
displayMode(FULLSCREEN)
function setup()
    --noSmooth()

    t = vec3(0.591177,-0.381549,1)

    nbTouches = 0
    touches = {}

    step = 0.0001
    NB_PASS = 20
    PASS = 0
    PASS_DIR = 0

    scene = Scene(WIDTH,HEIGHT)

    T_BEGAN = false
    T_ENDED = false

    zoomDir = -1

end

function draw()
    local sh = scene.surface.shader
    if PASS_DIR == 0 then PASS_DIR = 1 else PASS_DIR = 0 end
    sh.PASS_DIR = PASS_DIR

    if zooming then
        t.z = t.z + (1000/(100/t.z)/60*DeltaTime)*zoomDir
        sh.ZOOM = t.z
    end

    if T_BEGAN then
        NB_PASS = 1
        PASS = 0
        sh.resolution = vec2(128, 96)
        scene.buffer = image(128, 96)
        T_BEGAN = false
    end

    if T_ENDED then
        if nbTouches == 0 then
            scene:draw()
            NB_PASS = 20
            local w, h = spriteSize(scene.BG)
            setContext(scene.BG)
            sprite(scene.buffer, w/2, h/2, w)
            sh.resolution = vec2(scene.w, scene.h)
            scene.buffer = image(scene.w, scene.h)

            w, h = spriteSize(scene.buffer)
            setContext(scene.buffer)
            sprite(scene.BG, w/2, h/2, w)
            setContext()
        end
        T_ENDED = false
    else

        PASS = PASS + 1
        if PASS > NB_PASS then
            PASS = 1
        end

        sh.PASS = PASS
        sh.NB_PASS = NB_PASS

        scene:draw()
    end

end


--# Scene
Scene = class()

function Scene:init(w, h)

    --set to non retina
    self.w = w/2
    self.h = h/2

    self.buffer = image(self.w, self.h)
    self.BG = image(WIDTH, HEIGHT)

    local myShader = shader()
    myShader.vertexProgram = self:vertexShader()
    myShader.fragmentProgram = self:fragmentShader()

    self.surface = mesh()
    self.surface.shader = myShader

    local cos = math.cos
    local sin = math.sin

    local camPos = vec3(cos(t.y)*cos(t.x), sin(t.y), cos(t.y)*sin(t.x))*2
    local camUp = vec3(0., 1., 0.)
    local camTarget = vec3(0., 0., 0.)
    local camDir = (camTarget-camPos):normalize()
    local camSide = camDir:cross(camUp):normalize()
    local camCDS = camDir:cross(camSide)

    local sh = self.surface.shader
    sh.resolution = vec2(self.w, self.h)
    sh.NB_PASS = NB_PASS
    sh.PASS_DIR = PASS_DIR
    sh.PASS = PASS
    sh.sinTime = 1

    sh.camPos = camPos
    sh.camDir = camDir
    sh.camSide = camSide
    sh.camCDS = camCDS
    sh.ZOOM = t.z

    self.surface:addRect(0,0, WIDTH, HEIGHT)
end

function Scene:draw()
    --self.surface.shader.sinTime = .5+math.abs(math.sin(ElapsedTime))*2
    setContext(self.buffer)
    self.surface:draw()
    setContext()
    sprite(self.buffer, WIDTH/2, HEIGHT/2, WIDTH)
end

--# BuildShaders
function Scene:fragmentShader()
    local source =
[[
#define PI 3.141592653589793
#define BAILOUT sqrt(PI*.5)
#define MAX_ITER 10
#define MAX_DIST 2.2
#define MAX_MARCH 30

#define DETAIL 0.0025

precision highp float;

uniform float sinTime;
uniform float ZOOM;

uniform float NB_PASS;
uniform float PASS;
uniform int PASS_DIR;

uniform vec2 resolution;
uniform vec3 camPos;
uniform vec3 camDir;
uniform vec3 camSide;
uniform vec3 camCDS;



float mandelbulb(vec3 p, inout float color)
{
    vec3 v = p;
    vec3 c = v;

    float r=0.0;
    float d=1.0;
    for(int n=0; n<=MAX_ITER; ++n)
    {
        r = length(v);
        if(r>BAILOUT) break;

        d *= pow(r, 7.)*8.+1.;
        float zr2 = dot(v, v);
        float zq2 = dot(v.xz,v.xz);
        float izr2 = 1.0/zr2;
        float izq2 = 1.0/zq2;

        vec4 k2 = vec4( izr2*v.y*v.y, izr2*zq2, izq2*v.z*v.z, izq2*v.x*v.x );
        vec4 k4 = k2*k2;
        vec4 k6 = k4*k2;
        vec4 k8 = k4*k4;

        float ko = izr2*v.y*sqrt(zq2);
        float ki = izq2*v.z*v.x;

        float cozo8 = k8.x + k8.y - 28.0*k6.x*k2.y - 28.0*k2.x*k6.y + 68.0*k4.y*k4.x;
        float cozi8 = k8.z + k8.w - 28.0*k6.z*k2.w - 28.0*k2.z*k6.w + 68.0*k4.w*k4.z;
        float sizo8 = 8.0*ko*(k6.x - k6.y - 7.0*k2.y*k4.x + 7.0*k2.x*k4.y);
        float sizi8 = 8.0*ki*(k6.z - k6.w - 7.0*k2.w*k4.z + 7.0*k2.z*k4.w);

        float zr8 = zr2*zr2; zr8 = zr8*zr8;
        //v = zr8*vec3( sizo8*sizi8, cozo8, sizo8*cozi8 );
        v = zr8*vec3( sizo8, cozi8, cozo8 );

        //faster (no sqrt) but not really interesting
        //v = zr8*vec3( sizi8, cozi8, cozo8 );
        //v = zr8*vec3( sizi8, cozo8, cozi8 );

        v += c*sinTime;
        color = abs(v.x)+abs(v.y)+abs(v.z);
    }
    return .5*log(r)*r/d;
}


void main( void )
{
    vec2 pos = (gl_FragCoord.xy - resolution.xy) / resolution.x;

    if ( pos.x > 1. || pos.y > 1.)
    {
        gl_FragColor = vec4(0, 0, 0, 1);
        return;
    }


    float n;
    if (PASS_DIR == 1)
        n = mod(1. + PASS + gl_FragCoord.y + gl_FragCoord.x, NB_PASS);
    else
        n = mod(PASS + gl_FragCoord.y - gl_FragCoord.x, NB_PASS);

    if ( n != 0.)
    {
        gl_FragColor = vec4(0, 0, 0, 0); return;
    }

    vec3 rayDir = normalize((ZOOM*sinTime)*camDir + pos.x*camSide + pos.y*camCDS);

    float detail = (DETAIL/ZOOM)*2.;
    if (detail>.0025) detail = 0.0025;

    float m = 1.0;
    float dist = 0.0;
    float total_dist = 0.9;

    float rdt = dot(camPos, rayDir);
    float rdr = dot(camPos, camPos) - (1.3);
    float Disc = (rdt*rdt)-rdr;
    float col2 = 0.;
    float col = 0.;
    float color = 0.;
    if (Disc>0.)
    {
        for(int i=0; i<MAX_MARCH; ++i)
        {
            total_dist+=dist;
            dist = mandelbulb(camPos + rayDir * total_dist, color);
            m -= 1./float(MAX_MARCH);

            if(dist<detail || total_dist>MAX_DIST) break;
        }
        col2 = m*(2.-total_dist)*.5;
        col = m*color*.4;

    }
    if (total_dist>MAX_DIST) col = 0.;

    gl_FragColor = vec4(col, col2, col2, 1.);
}

]]
    return source
end

function Scene:vertexShader()
    local source =
[[
attribute vec4 position;
void main()
{
    gl_Position = position;
}
]]
    return source
end


--# Touches
--                   Touch Handling                 --
------------------------------------------------------

function touched(touch)
    local sh = scene.surface.shader
    local state = touch.state
    local touches = touches
    local deltaX = touch.deltaX
    local deltaY = touch.deltaY



    if state == ENDED then
        nbTouches = nbTouches - 1
        touches[touch.id] = nil
        if not zooming then T_ENDED = true end
    else
        if state == BEGAN then
            nbTouches = nbTouches + 1
            if not zooming then T_BEGAN = true end
        end

        touches[touch.id] = touch
    end

    if touch.tapCount == 2 and state == ENDED then
        zoomDir = -zoomDir
    end

    if touch.tapCount == 3 and state == ENDED then
        if zooming then
            zooming = false
            T_ENDED = true
        else
            zooming = true
            T_BEGAN = true
            T_ENDED = false
        end
    end

    local cos = math.cos
    local sin = math.sin
    local pi = math.pi


    if nbTouches == 1 then
        t.x = t.x + ((deltaX/WIDTH)*pi*2)/(t.z)
        t.y = t.y + ((deltaY/HEIGHT)*pi*2)/(t.z)

        if t.y>=pi/2 then t.y = t.y - (deltaY/HEIGHT)*pi*2 end
        if t.y<=-pi/2 then t.y = t.y - (deltaY/HEIGHT)*pi*2 end

        local camPos = vec3(cos(t.y)*cos(t.x), sin(t.y), cos(t.y)*sin(t.x))*2
        local camUp = vec3(0., 1., 0.)
        local camTarget = vec3(0., 0., 0.)
        local camDir = (camTarget-camPos):normalize()
        local camSide = camDir:cross(camUp):normalize()
        local camCDS = camDir:cross(camSide)

        sh.camPos = camPos
        sh.camSide = camSide
        sh.camDir = camDir
        sh.camCDS = camCDS
    end

-- Handle pinch to zoom
    if nbTouches == 2 then
        local newPt = {}
        local lastPt = {}
        local ins = table.insert
        for _, p in pairs(touches) do
            local px = p.x
            local py = p.y
            ins(newPt, vec2(px, py))
            if state == BEGAN then
                ins(lastPt, vec2(px, py))
            else
                ins(lastPt, vec2(px - p.deltaX, py - p.deltaY))
            end
        end

        local delta = lastPt[1]:dist(lastPt[2]) - newPt[1]:dist(newPt[2])
        t.z = t.z - delta/(200/t.z)
        if t.z < .1 then t.z = .1 end
    end


    sh.ZOOM = t.z
end
--# Backup
--[[
-- http://blog.hvidtfeldts.net/index.php/2011/11/distance-estimated-3d-fractals-vi-the-mandelbox/
#define MR2 .01
#define SCALE -2.
#define C1 abs(SCALE-1.)
#define C2 pow(abs(SCALE), float(1-MAX_ITER))
#define SCALE_VEC vec4(SCALE/MR2, SCALE/MR2, SCALE/MR2, abs(SCALE)/MR2)

float mandelbox(vec3 position, inout float color)
{
    vec4 p = vec4(position.xyz, 1.0);
    vec4 p0 = p;
    float n;
    for (int i = 0; i < MAX_ITER; i++)
    {
        p.xyz = clamp(p.xyz, -1.0, 1.0) * 2.0 - p.xyz;
        n = dot(p.xyz, p.xyz);
        p *= max(MR2/n, MR2);
        p = p*SCALE_VEC + p0;
    }
    return sqrt(n) / p.w;
}

--  slightly faster but also inacurate
--  http://www.iquilezles.org/www/articles/mandelbulb/mandelbulb.htm
float fastMandelbulb(vec3 p)
{
    vec3 v = p;
    vec3 c = v;

    float r=0.0;
    float d=1.0;
    for(int n=0; n<=MAX_ITER; ++n)
    {
        r = length(v);
        if(r>BAILOUT) break;

        d *= pow(r, 7.)*8.+1.;

        float x = v.x; float x2 = x*x; float x4 = x2*x2;
        float y = v.y; float y2 = y*y; float y4 = y2*y2;
        float z = v.z; float z2 = z*z; float z4 = z2*z2;

        float k3 = x2 + z2;
        float k2 = inversesqrt( pow(k3, 7.) );
        float k1 = x4 + y4 + z4 - 6.0*y2*z2 - 6.0*x2*y2 + 2.0*z2*x2;
        float k4 = x2 - y2 + z2;

        v.x = 64.0*x*y*z*(x2-z2)*k4*(x4-6.0*x2*z2+z4)*k1*k2;
        v.y = -16.0*y2*k3*k4*k4 + k1*k1;
        v.z = -8.0*y*k4*(x4*x4 - 28.0*x4*x2*z2 + 70.0*x4*z4 - 28.0*x2*z2*z4 + z4*z4)*k1*k2;

        v += c;

    }
    return .5*log(r)*r/d;
}
]]