High quality rendering of the Apollonian gasket
This is the third post in a series where I explore the mathematics of how it's possible to render an Apollonian gasket using modern fractal software. In the first one I worked out the first steps of the implementation, and in the second one how one can control the geometry and do basic animations.
The results of the second post consisted of an IFS to render gaskets, whose convergence wasn't very good, and was derived in a very formal manner, without geometrical considerations. It turns out these issues were related to each other, so that understanding better what the IFS does and improving it automatically made the quality better. The updates are all in the GitHub repository.
Improvements to the gasket creation parameters
Before moving to the quality improvements, I'd like to first point out some small improvements with respect to how the gasket is created. In the previous post, one needed to specify the values of the largest circles. That was a bit complicated, because it was pretty difficult to know a priori which combination of radii was valid. Because of this, I decided to explore the valid parameter space to see if I could parametrize it better. The valid parameters are the ones inside the deformed triangle seen in the graph below:
Why is this so? Well, the smallest possible value for the radius of the largest circle inscribed in the gasket is 2*sqrt(3)-3, which is around 0.464, and the largest is of course 1, the size of the unit disk where everything is inside. You can find this 2*sqrt(3)-3 value by applying Descartes' theorem where one curvature is -1 and the others are all equal. Now, if y denotes the second largest circle, it's bounded above by either x (which, by definition is the largest) or by 1-x, which is the largest circle that can fit in radius 1, together with x, whichever of this values is the smaller bounds x, therefore the y-values are bounded above by min(x,1-x).
Getting the lower bound is pretty simple once one has a geometrical intuition of what's happening. If you keep reducing y, eventually another circle will pop up which has larger radius, so the smallest value occurs when you have two circles of the same size. Applying Descartes' theorem for curvatures -1,1/x,1/y and 1/y, and solving for y we get the formula 4(x-1)x/(x+1)², which bounds the radii from below.
Because of this, I decide to parametrize this space in a natural way, with one parameter being 0 to represent the regular gasket, which has the smallest largest radius, and 1 the completely deformed one. Once this first parameter is chosen, we can find the radius value with linear interpolation, and I chose 0 to represent the top of the y-space, and 1 to represent the bottom.
Also, another point to notice is that this actually doesn't generate all gaskets. Even though we have a third parameter to rotate it, not all gaskets can be generated by specifying the circle sizes and rotating, because of a phenomenon that physicists like to call chirality. Intuitively, no matter how you rotate your right hand, it never is equal to your left hand. To make your right hand equal to the left, you need a "mirror" or "flip" transformation. Because of this, I included in the code the option to flip the gasket as well. It's pretty easy, one only needs to take the complex conjugate of the parameters.
Improving the quality
The way I was able to get better quality renders was by drastically improving the construction. If we take the Apollonian window, for example, which is the gasket generated by the ideal triangle at -1, 0 and 1. The set of ideal triangles that create the same gasket are the Kleininan group generated by the Möbius transformations A(z) = z/(-2i*z+1) and B(z)= (z+1)/(-3z+1).
The important thing about the B(z) transformation is that we have B(0)=1, B(1)=-1, and B(-1)=0. Applying B(z) three times gives us the identity transformation. We can also verify the following geometrical facts. When applied to the ideal triangle T=(-1,0,1), we have
- A(z) takes T to (-0.2+0.4i,0,0.2+0.4i)
- B⁻¹AB(z) takes T to (-1, -0.2+0.4i,i)
- BAB⁻¹(z) takes T to (i, 0.2+0.4i, 1)
We can understand this as: we have the first ideal triangle T, and then each of these three transformations will "shrink" it to a different ideal triangle. Just like the linear transformations in the Sierpinski triangle shrink a larger triangle into three smaller ones, the previous three Möbius transformations map the upper half plane into three distinct circles. Because of this, the upper part of the Apollonian window can be thought of as a distorted Sierpinski triangle.
To render the lower half, we simply replace A(z) by A⁻¹(z), and it can be verified, that the analogous transforms will map the initial ideal triangle T to the equivalent regions in the lower half plane. Then, why did the IFS in the previous post worked? Well, the IFS in fact consists of the functions f(z)=A(z) and g(z)=BA⁻¹(z). It turns out that the ideal triangles reached by all compositions of these functions (without their inverses) are the same as in the other set of transformations, that's why running the chaos game through these functions draws the gasket.
But why was the quality so poor? This problem was already observed in the flame algorithm paper. It's essentially because of symmetry. Good IFS's should consist of contracting transformations. While our IFS is minimal, we have (gf)³ = z, the identity transformation. Because of this, there are some random paths that do not bring the point closer to the limit set, and the convergence is poor, and even leads to washed out colors, if you're going to color the fractal.
In our case, the villain is the A(z) transformation. To render the whole gasket, we need both A(z) and A⁻¹(z) functions. But if they are applied in succession, they spoil the convergence of the IFS. The way to solve this is to use what the fractal art community often calls "xaos" or "relative weights". Fractal software like qosmic gives you the possibility to control the probability that a point in the chaos game goes to another transformation depending on the last one applied, essentially converting the chaos game into a Markov chain.
Using this feature, instead of the minimal two transforms, we can instead use six, three for the upper gasket and three for the lower one, and use relative weights so that every point that goes into the upper half never goes into the lower one. That way, convergence is very fast, and choosing color speed parameters carefully, we can make the color depend only on the last transformation applied, that's how the first figure in this post was plotted.
In fact, using relative weights, we can do even better. It can be noticed that intersection points of larger circles have some "holes". That's because to get there we need one transformation to be applied many times in sequence, and that's very unlikely to happen. We can tweak relative weights so that the last transformation applied is more likely to be applied again instead of the other two, and get the holes filled more.
Zooming into the gasket
The formalism developed here is very powerful. With the previous post, we could only control that the gasket would be rendered inside the complex unit disk, but with what we developed here, we have even more control over the regions of the complex plane where the points will end up. This allows us to apply some tricks in the style of Indra's pearls.
Applying successive transformations, we get smaller and smaller regions of the complex plane where we can be sure that parts of the gasket can be rendered. An example for the Apollonian window is given below:
This suggests a strategy to zoom into the gasket. We first iterate several transformations at random to find a point where we will not zoom into the black. Then, once we have that point, we can find the smallest circle that contains the entire screen, and then we set up the gasket starting from the ideal triangle defined by these points. By all our considerations, it's definitely the same gasket, and we get better quality because we only render what the screen is seeing.
To detect whether the screen is inside a given region, which can be either a half plane, the inside of a circle, or the outside of a circle, I used a technique similar to the signed distance field. It's implemented in the code below:
Basically, a circle or a line is always given by an equation of the form A(x²+y²)+Bx+Cy+D=0. When you are given the points, these coefficients can be calculated by solving a linear system, which means computing determinants. You can check if a point is "inside" or "outside" by substituting it into the formula and checking if the sign is positive or negative.
To check if the screen is inside a circle or a region of the plane, we can simply check if the vertices are all inside. To check if the screen is fully outside a circle requires more work, since the outside of a circle is not a convex region. This case is used only at the beginning of the zoom, though, and I just used standard circle rectangle intersection techniques for this case.
With all this done, I could successfully render a 4K video zooming into the gasket, showing its details in full glory. This is conclusive proof that the techniques here can be used for high quality rendering:
Conclusion
There's something very interesting about the video. To render it, I went as far as I could without having trouble with numerical precision. Although double floating point precision of around 12 decimal places is impressive for most applications, it's really nothing for fractals.
That's because since fractals are nearly scale invariant, to get the zoom working we need to increase the magnification scale exponentially, so we get to several digits really fast. In my case, the program started having issues finding where to place the transformations at around 100.000 scale. That makes sense, because the collision function computes determinants of the order of the area. At 10⁻⁵ distances, we are talking about areas of around 10⁻¹⁰! From this, the numerical errors start to accumulate and a lot of quality is lost.
I suppose this could be improved by using better collision algorithms, but I did not use the signed distance field technique randomly. The thing is, because of the algebraic closure of Möbius transformations, it is in principle possible to do everything with rational numbers, and using arbitrary precision libraries, everything done here could be replicated using integer arithmetic, allowing the zoom to go much much farther.
In any case, the actual bottleneck for this was to render all the frames. It took several hours to render 4K frames for the gasket, and the most immediate improvement would actually be using the GPU to render faster, before we could work on improving the precision.
