Rendering a mandelbrot zoom sequence with Julia

mandelbrot.jl

using ColorSchemes
using Dates
using Images
using Printf
using SIMD

const palette = RGB{Float32}.(colorschemes[:inferno].colors)
const n_colors = length(palette)
const black = RGB{Float32}(0, 0, 0)

function interpolate_color(n)
    if n == 0
        return black
    end

    c1 = palette[1 + Int32(floor(n * 10)) % n_colors]
    c2 = palette[1 + Int32(floor(n * 10 + 1)) % n_colors]
    frac = n % 1
    RGB{Float32}(c1 * frac + c2 * (1 - frac))
end

function iterate_mandelbrot(steps::Int16, c::Array{Complex{Float64}, 2})
    # Compute when each pixel diverges.
    set_zero_subnormals(true)

    flat_c = vec(c)
    z = flat_c .* 0
    escaped_at = zeros(Int16, size(z))

    n = length(z)

    for i = 1:n
        z_val = z[i]
        c_val = flat_c[i]
        @inbounds for step = 1:steps
            z_val = z_val ^ 2 + c_val
            if (abs(z_val) > 256)
                escaped_at[i] = Int16(step)
                break
            end
        end
        z[i] = z_val
    end
    reshape(z, size(c)), reshape(escaped_at, size(c))
end

const mid = complex(-0.743643887035763, 0.13182590421259918)

png_lock = ReentrantLock()

Threads.@threads for zoom = 0:1000
    start = now()

    scale = 4.5 * 1.025^(1-zoom)
    # scale = 4.5 * 1.5^(1-zoom)
    step = scale / 350.0
    im_width = 1920
    im_height = 1080

    x0 = ([0:im_width-1;] ./ (im_width - 1) .- 0.5) .* complex(scale, 0)
    y0 = ([0:im_height-1;] ./ (im_height - 1) .- 0.5) .* complex(0, scale / 16 * 9)

    na = [CartesianIndex()]

    c = mid .+ (x0[na,:] .+ y0[:,na])

    z, escaped_at = iterate_mandelbrot(Int16(2000), c)

    plot_start = now()

    # Interpolate colors according to https://en.wikipedia.org/wiki/Plotting_algorithms_for_the_Mandelbrot_set#Continuous_(smooth)_coloring
    filtered_z = ifelse.(escaped_at .== 0, 10, abs.(z))
    log_z = log.(filtered_z .* filtered_z) ./ 2
    nu = log.(log_z ./ (log(2))) ./ log(2)

    smooth_escaped_at = ifelse.(escaped_at .== 0, 0, escaped_at .+ 1 .- nu)

    mandelbrot = map(interpolate_color, smooth_escaped_at)
    mandelbrot = colorview(RGB, mandelbrot)

    filename = @sprintf "img/mandelbrot_%03d.png" zoom
    lock(png_lock) do
        save(filename, mandelbrot)
    end

    not_diverged = sum(escaped_at .== 0)

    println("zoom $zoom done: $(plot_start - start) to compute, $(now() - plot_start) to plot; $(not_diverged / length(escaped_at)) never escaped")
end

Second attempt, with an explicit for loop - somehow, this is 10x faster than broadcasting operations over the arrays:

z .= ifelse.(escaped_at .> 0, z, z.^2 .+ flat_c)
escaped_at .= ifelse.((escaped_at .== 0) .& (abs.(z) .> 256),
                      Int16(step), escaped_at)