Updated Mandelbrot Set

Updated the Mandlebrot example with Zoom support, multi-pass, and improved colour scheme. The code seems to be too long for this forum, so I’ll post it in 2 (?) separate pieces.

-- Mandelbrot set explorer, by Dr. Phillip Alvelda
-- it uses a few tricks to speed up the set calcs

-- Updated with Zoom support, multi-pass, and improved colour scheme
-- Herwig Van Marck

supportedOrientations(LANDSCAPE_ANY)

-- Setup all the variables we'll be using
function variableSetup()
    Color=1
    Binary=2
    Grayscale=3

    RMode = Color -- Default rendering mode. 
                  -- Can also be "Grayscale" or "Binary"

    tileNumber = 100 -- Image resolution. 
                     -- Try varying this between 50 to 700

    tilesize = WIDTH/tileNumber

    ManCalcStart = 0 -- Variables track time it takes to calc one line
    ManCalcEnd = 0
    ManCalcTime = 0
    gx,gy = 0
    MinRe = -2.0 -- Window borders in real coordinates
    MaxRe = 0.7
    MinIm = -1.35
    MaxIm = 1.35
    Re_factor = (MaxRe-MinRe)/(WIDTH -1) -- Scale factors 
                                         -- From real to window
                                        
    Im_factor = (MaxIm-MinIm)/(HEIGHT -1)

    MaxIterations = 255 -- Max iteration count for membership check
    MaxIterations_ = 255
    UCounter = 0
    count = 1

    red={} -- Tables to store the color maps for the set rendering
    green={}
    blue={}

    myImage = image(WIDTH,HEIGHT) -- The raster display image
    tmpImage = nil
    zooming=false
    Resolution=50
end

-- Print the initial instructions, setup the color map
function setup() 
    setInstructionLimit(0)
    
    variableSetup()
    
    watch("ManCalcTime") 
    watch("DeltaTime")
    
    iparameter("RMode",1, 3,1) 
    --iparameter("Resolution",50,800,50)
    iparameter("MaxIterations_",10,1024,255)
    
    print("This program plots an image of the Mandelbrot set using the image() class.\
")
    print("ManCalcTime is the time it takes to calculate the set membership and escape values.\
")
    print("Parameter: RMode sets the type of color map the set is rendered with.\
")
    print("Parameter: MaxIterations_ sets the level of image detail.")
    print("Higher iterations make prettier pictures at the expense of slower update.")
    
    InitColorMap()
    --ManSetCalc()
    zoom=Zoom()
end


function InitColorMap()
    for i=0,MaxIterations do
        if RMode == Grayscale then -- Grayscale
        
            red[i]=255-i*255/MaxIterations
            green[i]=255-i*255/MaxIterations
            blue[i]=255-i*255/MaxIterations
            
        elseif RMode == Color then 
            
            -- Color gradients designed to use their prime
            -- number mulltiples to create a range of varied 
            -- and uniqe color bands, even when MaxIterations
            -- gets very high.
            if (i==MaxIterations) then
                red[i]=0
                green[i]=0
                blue[i]=0
            else
                red[i] = math.abs((i*11) % 512-255)
                green[i] = math.abs((i*5) % 512-255)
                blue[i] = math.abs((i*7) % 512-255)
            end
        elseif RMode == Binary then 
            
            -- Alternating black and white
            if (i==MaxIterations) then
                red[i]=0
                green[i]=0
                blue[i]=0
            else
                if (i % 2) == 0 then
                    red[i]=255
                    green[i]=255
                    blue[i]=255  
                else 
                    red[i]=0
                    green[i]=0
                    blue[i]=0
                end
            end
        end         
    end    
end

function touched(touch)
    zoom:touched(touch)
end

function ManSetCalc() 
    ManCalcStart= os.clock() 
    local y, x, n, isInside
    local imf=Im_factor*tilesize
    local ref=Re_factor*tilesize
    local iterations
    
    -- note each call to ManSetCalc() only updates 
    -- a single line of the image at y=count
    y = myImage.height - count
    c_im = MinIm + y*imf      
      
    for x = 1, tileNumber-1 do
        c_re = MinRe + x*ref
        --check to see if x,y are inside the main 
        --cardioid or period2 bulb
          
        --in which case they are in the set and 
        --no iterative check is necessary
          
        --results in a 5x speedup when zoomed out
        q=(c_re-0.25)*(c_re-0.25) + c_im*c_im
        
        if ( ((c_re+1)*(c_re+1) + c_im * c_im < 0.0625) or 
              (q*(q+(c_re-0.25))<0.25*c_im*c_im) ) then         
            myImage:set(x,y,red[MaxIterations],
                            green[MaxIterations],
                            blue[MaxIterations],255)
        else
            -- Iterate the Mandelbrot set function from the starting
            -- position on the complex plane to see if it stays 
            -- inside a radius less than 2 units from the origin.
            Z_re = c_re
            Z_im = c_im
            isInside = true
            
            for n=0, MaxIterations, 1 do
                Z_re2 = Z_re*Z_re
                Z_im2 = Z_im*Z_im
                iterations=n
                -- Note the radius limit is 4 here after a rewrite
                -- to avoid having a sqrt() in the inner loop.
                if Z_re2 + Z_im2 > 4 then
                    isInside = false   
                    break
                end
                Z_im = 2*Z_re*Z_im + c_im
                Z_re = Z_re2 - Z_im2 + c_re
            end
             
            -- Store the iteration result in the image buffer           
            myImage:set(x,y,red[iterations],
                            green[iterations],
                            blue[iterations],255) 
        end
    end 
    
    ManCalcEnd = os.clock()
    ManCalcTime = ManCalcEnd-ManCalcStart
end

function restart()
    tileNumber=Resolution
    tilesize = WIDTH/tileNumber

    --tmpImage = myImage
    -- Redefine image buffer to new resolution
    myImage = image(tileNumber+1,tileNumber+1) 
        
    -- Reset rendering to top of image 
    -- to be sure count is within buffer bounds  
    count=1 
    LastResolution = Resolution
    LastRMode = RMode
end

Part two:

function draw()
    zoom:draw()
    smooth() -- noSmooth() if you like a blockier image
    background(0, 0, 0, 255)
    -- If RMode iparamter changes reinitialize the colormap
    if RMode ~= LastRMode then 
        InitColorMap()
        restart()
    end
    
    if MaxIterations_ ~= MaxIterations then
        MaxIterations=MaxIterations_
        InitColorMap()
        restart()
    end
    -- If Resolution iparameter changes then reinitialize image buffer
    --if Resolution ~= LastResolution then
    --    restart()
    --end
    
    if (tmpImage) then
        sprite(tmpImage,WIDTH*0.5,HEIGHT*0.5, WIDTH, HEIGHT)
    end
    -- Render the image buffer
    if zoom.started or zoom.started2 then
        sprite(myImage,WIDTH*0.5,HEIGHT*0.5, WIDTH, HEIGHT)
        zooming=true
    else
        if (zooming) then
            if (zoom.center==vec2(zoom.initx,zoom.inity)) and
                (zoom.offset == vec2(0,0)) and
                (zoom.zoom==1) then
                    print("Reset")
                    MinRe,MaxRe,MinIm,MaxIm = -2.0,0.7,-1.35,1.35
            else
            -- recalc coordinates using current zoomed window
              MinRe,MaxRe,MinIm,MaxIm=
                MinRe+(zoom.offset.x-(zoom.center.x)/zoom.zoom)*(MaxRe-MinRe)/WIDTH,
                MinRe+(zoom.offset.x-(zoom.center.x- WIDTH)/zoom.zoom)*(MaxRe-MinRe)/WIDTH,
                MinIm+(zoom.offset.y-(zoom.center.y)/zoom.zoom)*(MaxIm-MinIm)/HEIGHT,
                MinIm+(zoom.offset.y-(zoom.center.y- HEIGHT)/zoom.zoom)*(MaxIm-MinIm)/HEIGHT
            end
            Re_factor = (MaxRe-MinRe)/(WIDTH - 1)
            Im_factor = (MaxIm-MinIm)/(HEIGHT - 1)
            Resolution=50
            restart()
            zooming=false
        else
            pushMatrix()
            resetMatrix()
            sprite(myImage,WIDTH*0.5,HEIGHT*0.5, WIDTH, HEIGHT)
            popMatrix()
        end
    end
    
    -- Calculate and store current line of image data
    ManSetCalc() 
    count = 1 + count % tileNumber
    if (count==1) then
        if not(zooming) then
            zoom:clear()
        end
        if (Resolution<800) then
            Resolution = Resolution * 2
            tmpImage = myImage
            restart()
        end
    end
end

-- Zoom library
-- Herwig Van Marck
-- usage:
--[[
function setup()
    zoom=Zoom(WIDTH/2,HEIGHT/2)
end
function touched(touch)
    zoom:touched(touch)
end
function draw()
    zoom:draw()
end
]]--

Zoom = class()

function Zoom:init(x,y)
    -- you can accept and set parameters here
    self.touches = {}
    self.initx=x or 0;
    self.inity=y or 0;
    self:clear()
    print("Tap and drag to move\
Pinch to zoom")
end

function Zoom:clear()
    self.lastPinchDist = 0
    self.pinchDelta = 1.0
    self.center = vec2(self.initx,self.inity)
    self.offset = vec2(0,0)
    self.zoom = 1
    self.started = false
    self.started2 = false
end

function Zoom:touched(touch)
    -- Codea does not automatically call this method
    if touch.state == ENDED then
        self.touches[touch.id] = nil
    else
        self.touches[touch.id] = touch
        if (touch.tapCount==2) then
            self:clear()
        end
    end
end

function Zoom:processTouches()
    local touchArr = {}
    for k,touch in pairs(self.touches) do
        -- push touches into array
        table.insert(touchArr,touch)
    end

    if #touchArr == 2 then
        self.started = false
        local t1 = vec2(touchArr[1].x,touchArr[1].y)
        local t2 = vec2(touchArr[2].x,touchArr[2].y)

        local dist = t1:dist(t2)
        if self.started2 then
        --if self.lastPinchDist > 0 then 
            self.pinchDelta = dist/self.lastPinchDist          
        else
            self.offset= self.offset + ((t1 + t2)/2-self.center)/self.zoom
            self.started2 = true
        end
        self.center = (t1 + t2)/2
        self.lastPinchDist = dist
    elseif (#touchArr == 1) then
        self.started2 = false
        local t1 = vec2(touchArr[1].x,touchArr[1].y)
        self.pinchDelta = 1.0
        self.lastPinchDist = 0
        if not(self.started) then
            self.offset = self.offset + (t1-self.center)/self.zoom
            self.started = true
        end
        self.center=t1
    else
        self.pinchDelta = 1.0
        self.lastPinchDist = 0
        self.started = false
        self.started2 = false
    end
end

function Zoom:clip(x,y,w,h)
    clip(x*self.zoom+self.center.x- self.offset.x*self.zoom,
        y*self.zoom+self.center.y- self.offset.y*self.zoom,
        w*self.zoom+1,h*self.zoom+1)
end

function Zoom:text(str,x,y)
    local fSz = fontSize()
    local xt=x*self.zoom+self.center.x- self.offset.x*self.zoom
    local yt=y*self.zoom+self.center.y- self.offset.y*self.zoom
    fontSize(fSz*self.zoom)
    local xtsz,ytsz=textSize(str)
    tsz=xtsz
    if tsz<ytsz then tsz=ytsz end
    if (tsz>2048) then
        local eZoom= tsz/2048.0
        fontSize(fSz*self.zoom/eZoom)
        pushMatrix()
        resetMatrix()
        translate(xt,yt)
        scale(eZoom)
        text(str,0,0)
        popMatrix()
        fontSize(fSz)
    else
        pushMatrix()
        resetMatrix()
        fontSize(fSz*self.zoom)
        text(str,xt,yt)
        popMatrix()
        fontSize(fSz)
    end
end

function Zoom:draw()
    -- compute pinch delta
    self:processTouches()
    -- scale by pinch delta
    self.zoom = math.max( self.zoom*self.pinchDelta, 0.2 )

    translate(self.center.x- self.offset.x*self.zoom,
        self.center.y- self.offset.y*self.zoom)
    
    scale(self.zoom,self.zoom)

    self.pinchDelta = 1.0
end