Orbiting the Hénon Attractor
Note: This was originally published as an Observable notebook.
The Hénon strange attractor is a system of two simple equations that gives rise to chaotic behavior.
I first encountered Jared Tarbell's beautiful visualizations of the attractor many years ago and recently explored making my own pictures using my favorite programming language at the time, Julia.
Here's a gallery of the most interesting pictures I created along the way. The code for these is described at the bottom of the post.
Click on any picture to see it in full resolution.
Behind the Scenes
Here's the Julia code that I wrote to explore this system.
1. Introducing the attractor
To explore a chaotic attractor we need to define it in code. The Hénon map has an angle parameter henon(α)
to return the attractor function for that angle.
# define the attractor
function henon(α)
c = cos(α)
s = sin(α)
attractor((x, y)) = (x*c  (y  x^2)*s, x*s + (y  x^2)*c)
end
We'll also need some utilities to help us along the way.
First, a function to generate random points
# randomly initialize an point within our bounds, then
# iterate it a few times to remove initial noisy jumps
function init(h, (xlo, xhi), (ylo, yhi))
init_x = xlo + (xhi  xlo)rand()
init_y = ylo + (yhi  ylo)rand()
p = (init_x, init_y)
for _ in 1:100
p = h(p)
end
p
end
Second, a linear scale function, similar to d3.scaleLinear, which interpolates values in the interval
# map from continuous (lo, hi) to integer 1:n
function lerp((lo, hi), n)
@assert lo < hi
x > 1 + round(Int, (x  lo) / (hi  lo) * (n  1))
end
Finally, a simple fuzz
function to randomly jitter points. This technique is responsible for the organic/analog feeling of the images, particularly in deep zooms.
function fuzz(x, amount)
x + amount * (rand()  0.5)
end
2. Computing orbits
Now we get to the heart of the system — it samples Hénon orbits and returns a 2D array of counts indicating how many times each pixel bin was visited by an orbit.
function orbit_counts(
α; # the henon angle parameter
sz = 1000, # size of the output canvas
n_samples = 100_000, # number of total points we want to record on our canvas
max_samples_per_orbit = 50_000, # the maximum number of samples for any individual orbit
x_bounds = (1.0, 1.0), # the x zoom bounds
y_bounds = (1.0, 1.0), # the y zoom bounds
fuzz_factor = 0.0 # a fuzz factor to jitter points before plotting
)
# create a henon attractor with our given angle
h = henon(α)
# linear scales, to integer coordinates in A
ix, iy = lerp(x_bounds, sz), lerp(y_bounds, sz)
# output array of visit counts for each bin
counts = zeros(sz, sz)
# loop until we collect the desired number of attractor samples
count = 0
while count < n_samples
# initialize a point
p = init(h, x_bounds, y_bounds)
for _ in 1:max_samples_per_orbit
# iterate the attractor
p = h(p)
# if the orbit diverges, break out of the loop for this orbit
if any(isinf, p)  any(isnan, p)
break
end
# if the point is within bounds, increment the count in the corresponding bin
if x_bounds[1] ≤ p[1] ≤ x_bounds[2] && y_bounds[1] ≤ p[2] ≤ y_bounds[2]
x = clamp(ix(fuzz(p[1], fuzz_factor)), 1, sz)
y = clamp(iy(fuzz(p[2], fuzz_factor)), 1, sz)
counts[y, x] += 1
count += 1
end
end
end
# return the array of counts
counts
end
3. Making pictures
Next, let's generate six 1,000 × 1,000 arrays of orbit counts to visualize:
orbits = [
orbit_counts(
1.5957,
sz = 1000,
n_samples = 1_000_000,
max_samples_per_orbit = 100_000,
x_bounds = (.45, .55),
y_bounds = (0.05, .15),
fuzz_factor = 0.0001
)
for _ in 1:6
]
To make our pictures more interesting we generate multiple "layers" of counts and composite them together into the final image. This lets us take advantage of the variation between layers to make our pictures more interesting.
To talk about colors and color spaces we're going to use the Julia Images package:
using Images
We define one last utility function, after which we're ready to visualize the Hénon attractor:
# normalize an array of positive numbers to (0, 1), being
# careful not to introduce NaNs via division by zero.
normalize(A) = A ./ max(1e6, maximum(A))
Finally, the pièce de résistance! The visualize
function takes a collection of 2D orbit arrays as input and produces an image as output:
function visualize(orbits)
# normalize each layer of orbit counts to (0, 1), then
# concatenate them to create a 3D array of weights
weights = cat(normalize.(orbits)..., dims=3)
# try visualizing each pixel based on the difference between its
# minimum value and maximum value across all of the orbit count layers.
# there are other approaches; this is a simple one to start with.
A = normalize(maximum(weights, dims=3) . minimum(weights, dims=3))
# we're going to use the LCH color space.
# lightness goes from 0100
L = 100(1 . A)
# contrast goes from 0300ish
C = 300A
# hue cycles with a period of 360 degrees
H = 360normalize(sum(weights .^ 3, dims=3))
# create an array of colors by combining our lightness, contrast, and hue
colors = RGBA.(LCHuv.(L, C, H), 1)
# we want to return a 2D array. By this point our array is
# width × height × 1, so we can just drop the third dimension
dropdims(colors, dims=3)
end
# visualize the six orbit count matrices we computed
visualize(orbits)
All of the pictures in this post were made with variations of this code.
Interesting directions

Use contrast normalization techniques, such as histogram equalization, which adjusts an image so that the full range of intensities is equally represented. Julia has a a great implementation in the ImageContrastAdjustment package.

Try applying tone mapping technques, such as Reinhard tone mapping, for which a simple implementation is
# reinhard tone mapping operator function tonemap(c, c_white) c * (1 + c / c_white^2) / (1 + c) end

Experiment with custom blend modes. Since the
visualize
function controls compositing, we can use standard blending operations in RGB or write our own.
That's it. If you make something cool, let me know!
Yuri Vishnevsky
July 2020