mathr / blog / #

Cardioid and bulb checking

Implicit Curves

The Mandelbrot set has two prominent solid regions. There is a cardioid, which is associated with fixed point (period 1) attractors, and a circle to the left, which is associated with period 2 attractors. The rest of the cardioid- and circle-like components in the Mandelbrot set are distorted.

These shapes can be described as implicit functions. For example, the circle is centered on \(-1+0i\) and has radius \(\frac{1}{4}\), and the function \[C_2(x, y) = (x - (-1))^2 + (y - 0)^2 - \left(\frac{1}{4}\right)^2\] is negative inside, zero on the boundary, and positive outside, and the same applies to the more complicated function for the cardioid: \[C_1(x, y) = \left( \left(x - \frac{1}{4}\right)^2 + y^2 \right)^2 + \left(x - \frac{1}{4}\right) \left(\left(x - \frac{1}{4}\right)^2 + y^2\right) - \frac{1}{4} y^2 \]

These implicit functions can be used to accelerate Mandelbrot set rendering. You can test if each \(c=x+iy\) is in the cardioid or circle quickly and easily, saving iterating the pixel all the way to the maximum iteration count (being interior to the Mandelbrot set means iterations of \(z \to z^2 + c\) will never escape to infinity).

Bounding Boxes

But we can accelerate further. If the whole viewing rectangle of a zoomed in view is far away from the circle and cardioid, then these per-pixel cardioid and circle tests are a waste of time, as they will never say they are inside. By analysing the coordinates of an axis-aligned bounding box (AABB) it's possible to decide when it's worth doing per pixel tests (that is, when the boundary of the shapes passes through the box - otherwise it's 100% interior or (more likely) 100% exterior to the shapes). This can save \(O(W H)\) work.

For example for the circle, if the lower edge of the box is above \(\frac{1}{4}\), or the upper edge of the box is below \(-\frac{1}{4}\), or the right edge of the box is left of \(-\frac{5}{4}\), or the left edge of the box is right of \(-\frac{3}{4}\), clearly it cannot overlap the circle. And if all corners are inside the circle, the whole box must be inside the circle. If some corners are inside and some are outside, then the boundary passes through.

But if all corners are outside it gets complicated: the box could be surrounding the whole circle, or a bulge of the circle could pop into an edge of the box. So the next step is to consider the vertices of the circle (the points most left/right/top/bottom): axis alignment means that if a bulge pops into the box, the vertex must lead the way. All in all there are many cases to consider, but it's not insurmountable.

Similarly for the cardioid, with the added complication that the vertices are not rational. However, squaring the coordinates does give dyadic rationals, so comparing \(y^2\) with \(\frac{3}{64}\) and \(\frac{27}{64}\) can do the trick.

Perturbation of Implicit Curves

For deep zooms, coordinates need high precision (lots of digits, most of which are the same for nearby pixels). Perturbation techniques mean using a high precision reference, with low precision differences to nearby points. This can also be applied to the implicit functions for the cardioid and circle: symbolically expand and cancel the large terms \(X, Y\) leaving only small terms of the scale of \(x, y\) in: \[c(X, Y, x, y) = C(X + x, Y + y) - C(X, Y)\] then evaluate \(C(X, Y)\) in high precision, round it to low precision, and add \(c(X, Y, x, y)\) evaluated in low precision. For accuracy, some coefficients in \(c\) will need to be calculated at high precision before rounding to low precision.

For example, the cardioid: \[c_1(X, Y, x, y) = a_x x + a_y y + a_{x^2} x^2 + a_{xy} xy + a_{y^2} y^2 + a_{x^3} x^3 + a_{x^2y} x^2y + a_{xy^2} xy^2 + a_{y^3} y^3 + x^4 + 2x^2y^2 + y^4 \] where \[a_x = (32XY^2+32X^3-6X+1)/8; a_y = (32Y^3+(32X^2-6)Y)/8; a_{x^2} = (16Y^2+48X^2-3)/8; \ldots \] I used wxMaxima to find all of the \(a\) coefficients, and calculate them one time per view using the reference, along with \(C_1(X, Y)\). For accuracy, with fixed point calculations you need about 4 times the number of fractional bits for intermediate calculations, and the values at low precision need to be relatively high accuracy (in my tests 24 bits was not enough to achieve good images, 53 seemed ok, and I use 64 bits just to be safe).

Perturbation of Bounding Boxes

The previous discussion about rejecting interior checks for the whole view can be applied with perturbation too, but some magic numbers need to be calculated at high precision before rounding to low precision, namely the special points (vertices and cusps): \[ X + \frac{5}{4}, X + 1, X + \frac{3}{4}, X + \frac{1}{8}, X - \frac{1}{4}, X - \frac{3}{8} \] and \[ Y + \frac{1}{4}, Y - \frac{1}{4}, Y^2 - \frac{3}{64}, Y^2 - \frac{27}{64} \] the addends are all dyadic rationals so can be represented exactly in binary fixed point or floating point.

Parametric Curves and Distance Estimates

The circle and cardioid also have parametric forms. Here's the cardioid: \[C_1(t) = \left(\frac{\sin(t)^2-(\cos(t)-1)^2+1}{4}, \frac{(\cos(t)-1)\sin(t)}{2}\right) \] If you could work out the distance to the nearest point of the curve, then all views with a smaller circumradius and same center would be 100% exterior (or 100% interior). For the cardioid, considering that the dot product of the tangent of the curve at the nearest point and the vector from the point to the nearest point of the curve must be zero (perpendicular), it means solving \[(4\cos(t)-4\cos(2t))y+(4\sin(2t)-4\sin(t))x-\sin(t)=0\] which can be rearrange to a 9 or 10 degree polynomial in \(\tan(t)\) using trigonometric identities. This is altogether a hard problem to solve, most practical is bisection of the trigonometric form on segments between by \(t = k\frac{\pi}{3}\). Linear distance estimate using Taylor expansion gives a closed form \(d(x,y)\), but it's not accurate, especially near the cusp. Quadratic Taylor expansion gives a high degree polynomial to solve. On the other hand, exact distance to a circle is easy.

Bibliography