Skip to content

Instantly share code, notes, and snippets.

@Mononofu
Created October 10, 2020 21:39
Show Gist options
  • Save Mononofu/60263aa6caf07d9228c4342f37edbf93 to your computer and use it in GitHub Desktop.
Save Mononofu/60263aa6caf07d9228c4342f37edbf93 to your computer and use it in GitHub Desktop.
Rendering a mandelbrot zoom sequence with Julia
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
@Mononofu
Copy link
Author

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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment