Distance estimation.
authorClaude Heiland-Allen <claude@mathr.co.uk>
Sat, 21 Dec 2013 15:56:28 +0000 (15:56 +0000)
committerClaude Heiland-Allen <claude@mathr.co.uk>
Sat, 21 Dec 2013 15:56:28 +0000 (15:56 +0000)
commite7ff7b7158c999b8fcc24e627c39226233181e38
treed886f24eadb9037d5eea5131bcd8cb8846a8b4f6
parent6255e2b5e432f65e1318d0704b9dbe3b1b0041a8
Distance estimation.

Our images look noisy and grainy near the boundary of the Mandelbrot set.
The escape time bands get closer and closer, while the pixel spacing is
fixed.  The pixel grid samples isolated points of a mathematically abstract
image defined on the continuous plane.  The Nyquist-Shannon sampling theorem
shows that sampling isolated points from a continuum is a valid approximation
only so long as the values don't change too quickly between the points.
Aliasing occurs when the values do change too quickly compared to the
sampling rate, with the grainy noisy visual effects as we have seen.
Because the escape time bands increase in number without bound as we
approach the boundary of the Mandelbrot set, no sampling rate can be
high enough.

We can measure the rate of change using differential calculus.  If we have
a function f(z), we can measure how the output of f changes when we change its
input.  The ratio (f(z+h)-f(z))/h with small h tells us if f expands or contracts
the region around z.  This is made mathematically precise by taking the limit
of this ratio as h tends to 0.  For example, the differential of z^2 with
respect to z can be calculated like this:

    ((z + h)^2-z^2)/h = (z^2 + 2 h z + h^2 - z^2)/h = 2 z + h -> 2 z

We iterate a function of two variables F(z,c) = z^2 + c, so there are two
possible derivatives: one with respect to z, and one with respect to c.  We
want to differentiate with respect to c, to see how F expands or contract
near each point in our image.  Our iteration is defined by

    F^{n+1} = F(F^n(z, c), c) = F^n(z,c)^2 + c

We can use mathematical induction: suppose we know the derivative dc_n for F^n.
By the chain rule, the derivative of a composition of functions is the product
of the derivatives at each step.  We're composing F^n with squaring, and we
know the derivatives for both, moreover the derivative of c with respect to c
is 1, so we end up with

    dc_{n+1} = 2 F^n(z, c) dc_n + 1

F^n(z,c) is just the iterate z_n, so we can calculate the derivative as we
go along.  We now just need a base case for our induction, and we start with
z_0 = 0 which doesn't involve c at all, so dc_0 = 0 too.  We store the final
dc value in our augmented pixel struct.

When it comes to saving the image, we have the same problem as before: how to
scale the value.  We know the escape time bands get closer together, and we
know this causes problems related to the pixel spacing.  So, multiplying the
derivative by the pixel spacing is a first step (the pixel spacing is in fact
the derivative with respect to our image device coordinates).  We know the
problems occur when the bands are closer, meaning the variation is more
rapid, so we divide something by this product.  What we divide by will remain
somewhat magical for now, involving the final z iterate, but our final result
is a distance estimate for the Mandelbrot set scaled to image device
coordinates.

This still is far from the [0,1] output range we need for PPM,
but we can choose how many bits of precision we want to keep for small values
(12 in our case), which gives us a fixed-point number (essentially an integer
type with a scaling factor that remains constant).  We clamp the value to
prevent wrapping around from putchar() truncation which would now give us
meaningless values for overflow.

The distance estimate has the useful property (proved in the Koebe 1/4
theorem) that it tells us roughly how far our point is from the boundary of
the Mandelbrot set.  We can use it to mask the grainy noise from aliasing,
with the pleasing side-effect of highlighting the intricate shape of the
boundary of the Mandelbrot set with all its filaments and shrunken copies.
To contrast with the dark boundary, we fill the interior with white.

    $ for exp in $(seq -w 0 16)
      do
        ./render \
          -1.2864956446267794623797546990264761111 \
           0.4339437613619272464268001250115015311 \
          2e-${exp} 512 ${exp} &&
        ./colour ${exp} > ${exp}.ppm
      done
colour.c
mandelbrot.h
mandelbrot_imp.c
render.c