Categories

## Astronomical Depth Buffer

Guess what. We’re going to talk about floating point arithmetic and what happens when using distances that range from the millimeter to Pluto’s distance to the Sun. Yes, yet again.

## Floating Point Arithmetic

OpenGL uses single precision floating point arithmetic. In any case, that’s what your GPU prefers to use. This means you can be as precise as 60 parts in a billion, but not more. That’s more than enough for drawing human-sized objects: the error is much smaller than the width of human hair. And it causes no problem when drawing the Solar-system: the error to Pluto’s orbit is 352 km, which is not visible when looking at the whole thing. But you are going to have issues with drawing a human if you do not know if it’s right here, or 352 km away.

We can avoid the main problem by doing most computations in double precision and centering around the camera before switching to simple precision for the rendering, as I briefly touched in a previous article. And that’s fine! We do not need more than single point precision to locate a pixel on the screen after all.

So far so good. But using a wide range of coordinates in the scene surfaces another issue: the naive use of the depth-buffer.

## The Depth-Buffer

3D rendering, at its core, is not very complex. You ask to draw some triangles, your GPU does some math to know where it should draw it on your screen and give colors to the pixels that end up in it. And that’s it! Well, almost. There is one tiny remaining problem.

Each pixel might be in the direction of various objects. But, most of the time, it should be of the color of the closest object. If we just draw all the objects overwriting the pixels each time, we might show something which is supposed to be hidden!

The most intuitive way to fix this is to just draw the objects in the appropriate order! This is z-sorting. But it is not a trivial task, and will not work in all cases (when A is partially behind B, B partially behind C and C partially behind A).

The more robust way is to remember for each pixel, how far it comes from. When drawing something over this pixel, we check whether the new value is closer. If it is, we use the new value, otherwise we keep the old one. This is the depth-buffer. For each pixel, we’ll track its current color (red, green, blue and alpha), and its current depth.

## Precision Issues

By default, OpenGL actually sets the depth-buffer (d) as the inverse of the distance from camera (z): d = 1 / z. This is a secondary effect of the perspective transform, which makes things look actually 3D.

There are various ways to fix this. The Outerra blog offers two approach for a logarithmic depth buffer:

1. explicitly set the depth-buffer from z (or w) in the fragment shader
2. override z when it will not affect perspective (which relies on w) and let OpenGL use it to set the depth-buffer

The first method works very well, but disables various optimizations, so it slows down rendering a bit. The second method lets OpenGL do its thing — so it is fast — but it does not work well with large triangles close to the camera.

For more information, you can also have a look at another article from the Outerra blog.

Writing this article has forced me to better understand the technique, and to fix a bug due to using the second method incorrectly (I had kept the multiplication by w from an earlier method).

With this out of the way, you can draw things very close-by or extremely far away in a single scene! No need for hacks such as splitting the scene in a close-up part and a far-away part!

Categories

## Solving Kepler’s Equation 5 Million Times a Second

I know it is a controversial opinion, but you might need to know where things are when running a simulation. As a bonus, it helps you in knowing what to draw on the screen.

Of course, you can just use a position vector (x, y, z) where x, y and z are the coordinates of the object. But things have the inconvenient tendency not to stay in the same place for all eternity. So we’ll want to know where that object is at the given time.

Now, the more generic approach is to run a physics simulation. We compute the acceleration of the object, update its speed depending on the acceleration and update its position depending on its speed. By the way, use Runge-Kutta 4 for this, not Euler’s method.

That would work for things in orbit too. But we can do much better in this case. This one guy, Johannes Kepler found simpler formulas to describe the movement of two orbiting objects (e.g. the Sun and the Earth, or the Earth and the Moon). Basically, the trajectories will be ellipses, and we can deduce the position of an object in its orbit at any given time.

We can describe the position of the object in its orbit as the angle between the direction of the periapsis (where the planet is closest to its sun) and the direction to the object, taken at the relevant focus. This is the true anomaly. Everything is an illusion, so we also like to pretend the ellipse is actually a circle and define the eccentric anomaly. But these angles are not linear with time. So we also define the mean anomaly, which works like an angle (goes from 0 to 2π) but is linear with time. Once we know the direction of the object, we will only need to know its distance, r.

And here are the formulas to get from time to the true anomaly, and how we then obtain the distance from focus (e is the eccentricity of the ellipse):

• mean anomaly (M): M = M₀ + n × (t – t₀)
• eccentric anomaly (E): M = E – e × sin E
• true anomaly (f): $latex \displaystyle f=2\arctan\sqrt{\frac{1+e}{1-e}} \tan\frac E2$
• distance from focus (r): $latex \displaystyle r=\frac p{1 + e \cos f}$

You won’t be calculating this in your head, but the expressions are pretty easy to translate to working code. Well, except for one. Notice how the eccentric anomaly is not expressed as a function of the mean anomaly? Well, that’s the hard part. And it’s called Kepler’s equation.

We do not know a simple way to express E as a formula depending only on e and M. Instead, we must proceed by successive approximations until we decide we are close enough. For this, we will need two things:

• a starting value E₀ for E
• a way to improve the current approximation of E

Most introductory books and online articles I have found will tell you to take E₀ = M. This has been obsolete since 1978. It works well enough, but fails in corner cases (near-parabolic orbits around periapsis). And modern starting values are close enough to limit the number of improvements you then need to make. I reviewed 21 papers, selected the 3 most relevant ones, and implemented 21 methods for picking starting values (3 papers for the hyperbolic case with 5 methods for the hyperbolic case).

To improve the starting value, old references will make you do E := M + e sin E, but this is pretty naive. Most modern references will use Newton’s method, which helps you solve an equation faster as long as you know its first derivative. But you can go further by using the second derivative using Halley’s method. You can go arbitrarily far using Householder’s method.

Great, now I have 21 × 4 = 84 approaches to try (5 × 3 = 15 for the hyperbolic case). Which one is the most accurate? Which one is the fastest? And, more importantly, can I find a fast and accurate approach?

You can spend years studying the theoretical properties of these functions. Or, we just try them for real. It is the most efficient way to know how they actually behave on real hardware with real implementations of numbers.

I wrote a benchmark tool to test both the accuracy and the speed of each approach. It turns out that only a handful of approaches are accurate enough for my taste:

# Elliptical case
## Naive method
## Newton's method
## Halley's method
580 c  3.580e-16  gooding_10
764 c  3.926e-16  gooding_11
577 c  3.604e-16  mikkola_1
742 c  3.653e-16  mikkola_2
## Householder's third order method
606 c  3.645e-16  gooding_10
809 c  3.623e-16  gooding_11
641 c  3.935e-16  mikkola_1
764 c  3.904e-16  mikkola_2

# Hyperbolic case
## Newton's method
## Halley's method
1454 c  3.731e-15  mikkola_1
1578 c  3.725e-15  mikkola_2
## Householder's third order method
1704 c  3.618e-15  mikkola_1
1757 c  3.618e-15  mikkola_2

The first column is the number of CPU cycles per solve of Kepler’s equation. The second column is the worst relative error over the test set. The third column is a short name for the method for picking a starting value for E.

As we can see above, no approach seems to work well enough with the naive improvement, or with Newton’s method. And using third derivatives does not do much. gooding_11 and mikkola_2 are relatively slow. This leaves gooding_10 and mikkola_1. I picked mikkola_1 because the paper also treated the hyperbolic case. It comes from A cubic approximation for Kepler’s equation by Seppo Mikkola.

And after all this, I have a very robust way to know where stuff is at any time! Of course, you should not forget to take care of floating point precision issues when implementing all of this…

## Appendix: Full benchmark

# Elliptical case
## Naive method
264 c  2.500e-01  smith_1
276 c  5.303e-01  smith_2
360 c  2.000e-01  smith_3
444 c  9.701e-02  smith_4
385 c  1.667e-01  smith_5
480 c  6.495e+01  smith_6
246 c  2.500e-01  gooding_1
299 c  2.000e-01  gooding_2
356 c  1.667e-01  gooding_3
246 c  5.303e-01  gooding_4
460 c  9.701e-02  gooding_5
253 c  2.417e-01  gooding_6
325 c  7.497e-02  gooding_7
334 c  7.497e-02  gooding_7b
378 c  1.277e-01  gooding_8
423 c  7.499e-02  gooding_9
455 c  2.376e-01  gooding_10
608 c  8.020e-02  gooding_11
286 c  1.439e-01  gooding_12
448 c  4.342e-01  mikkola_1
629 c  4.426e-01  mikkola_2
## Newton's method
522 c  8.517e-03  smith_1
539 c  1.379e-02  smith_2
578 c  8.111e-03  smith_3
697 c  8.111e-03  smith_4
706 c  8.111e-03  smith_5
900 c  1.028e+02  smith_6
481 c  8.517e-03  gooding_1
583 c  8.111e-03  gooding_2
561 c  8.111e-03  gooding_3
476 c  1.379e-02  gooding_4
686 c  8.111e-03  gooding_5
537 c  3.528e-02  gooding_6
536 c  1.378e-02  gooding_7
585 c  1.378e-02  gooding_7b
702 c  2.524e-04  gooding_8
680 c  1.379e-02  gooding_9
731 c  1.117e-10  gooding_10
940 c  2.072e-08  gooding_11
528 c  5.673e-03  gooding_12
697 c  1.825e-06  mikkola_1
815 c  1.825e-06  mikkola_2
## Halley's method
395 c  1.250e-01  smith_1
419 c  2.553e-03  smith_2
488 c  6.250e-02  smith_3
596 c  1.710e-02  smith_4
557 c  4.167e-02  smith_5
989 c  5.130e+01  smith_6
429 c  1.250e-01  gooding_1
503 c  6.250e-02  gooding_2
524 c  4.167e-02  gooding_3
426 c  2.553e-03  gooding_4
587 c  1.710e-02  gooding_5
410 c  6.127e-03  gooding_6
467 c  2.548e-03  gooding_7
440 c  2.548e-03  gooding_7b
508 c  9.153e-05  gooding_8
552 c  2.550e-03  gooding_9
598 c  3.580e-16  gooding_10
821 c  3.926e-16  gooding_11
430 c  1.028e-03  gooding_12
594 c  3.604e-16  mikkola_1
742 c  3.653e-16  mikkola_2
## Householder's third order method
438 c  1.562e-02  smith_1
439 c  6.779e-04  smith_2
501 c  7.812e-03  smith_3
667 c  2.138e-03  smith_4
533 c  5.208e-03  smith_5
1053 c  6.535e+01  smith_6
453 c  1.562e-02  gooding_1
501 c  7.812e-03  gooding_2
594 c  5.208e-03  gooding_3
429 c  6.779e-04  gooding_4
595 c  2.138e-03  gooding_5
454 c  1.638e-03  gooding_6
526 c  6.736e-04  gooding_7
585 c  6.736e-04  gooding_7b
583 c  9.493e-09  gooding_8
652 c  6.741e-04  gooding_9
681 c  3.645e-16  gooding_10
846 c  3.623e-16  gooding_11
486 c  2.724e-04  gooding_12
651 c  3.935e-16  mikkola_1
752 c  3.904e-16  mikkola_2

# Hyperbolic case
## Newton's method
780 c  7.245e+25  basic
1064 c  1.421e-03  mikkola_1
1201 c  1.970e-04  mikkola_2
1014 c  8.559e-01  gooding
892 c  1.515e+19  pfelger
## Halley's method
1134 c  1.729e+17  basic
1396 c  3.731e-15  mikkola_1
1526 c  3.725e-15  mikkola_2
1265 c  9.160e+10  gooding
1237 c  1.791e+13  pfelger
## Householder's third order method
1260 c  1.323e+07  basic
1631 c  3.618e-15  mikkola_1
1595 c  3.618e-15  mikkola_2
1467 c  2.276e+17  gooding
1365 c  9.997e+01  pfelger