Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

In this Discussion

Updated Mandelbrot Set

edited January 2012 in Code Sharing Posts: 118

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.\n")
    print("ManCalcTime is the time it takes to calculate the set membership and escape values.\n")
    print("Parameter: RMode sets the type of color map the set is rendered with.\n")
    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
Tagged:

Comments

  • Posts: 118

    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\nPinch 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
    
    
Sign In or Register to comment.