Rendering Implicit Surfaces and Distance Fields: Sphere Tracing

News (August, 31): We are working on Scratchapixel 3.0 at the moment (current version of 2). The idea is to make the project open source by storing the content of the website on GitHub as Markdown files. In practice, that means you and the rest of the community will be able to edit the content of the pages if you want to contribute (typos and bug fixes, rewording sentences). You will also be able to contribute by translating pages to different languages if you want to. Then when we publish the site we will translate the Markdown files to HTML. That means new design as well.

That's what we are busy with right now and why there won't be a lot of updates in the weeks to come. More news about SaP 3.0 soon.

We are looking for native Engxish (yes we know there's a typo here) speakers that will be willing to readproof a few lessons. If you are interested please get in touch on Discord, in the #scratchapixel3-0 channel. Also looking for at least one experienced full dev stack dev that would be willing to give us a hand with the next design.

Feel free to send us your requests, suggestions, etc. (on Discord) to help us improve the website.

And you can also donate). Donations go directly back into the development of the project. The more donation we get the more content you will get and the quicker we will be able to deliver it to you.

20 mns read.

The Maths of Sphere-Tracing

To understand the sphere-tracing algorithm you need to first remember the fact that the shapes that we will render with this method can be defined mathematically or implicitly. Recall that for instance the implicit form of the sphere is:

$$x^2 + y^2 + z^2 - r^2 = 0,$$

Figure 1: A is the set of points making up the surface of the implicit sphere. They all satisfy the equation $x^2 + y^2 + z^2 - r^2 = 0$.

where \(x\), \(y\) and \(z\) are the coordinates of a point in 3D space. All it says is that points that actually lie on the surface of the sphere, in other words all points making up the surface of the sphere or its isosurface, satisfy this equation. In other words, if you raise the point coordinates to the power of two sum them up, subtract the radius of the sphere raised to the power of two as well, then the result is 0 for all points lying on the surface of the sphere. All these points belong to a set that we will call \(A\) (figure 1).

Now once you understand certain shapes like sphere, cubes, cones, cylinders (and many more, we will show examples later on) can be represented with equations, the next step in developing a method similar to sphere-tracing, is to answer the following question:

"Is it possible to find a equation that returns the distance from a point (let's call it \(P\)) in space to the nearest point on the surface of an implicit shape?"

Figure 2: a DIF is a function that returns the distance from any point to the nearest point on the surface of the implicit object. A DUF is a function that returns a distance equal to or lower than the distance returned by the DIF.

In other words, we need to find the point in the set \(A\) that is the closest to \(P\). In figure 1, we denoted this equation \(d(P, A)\). This is not the mathematical notation used by Hart, but hopefully this is more meaningful if you are not a mathematician. Technically, if we can find an equation that returns the distance from \(P\) to the nearest point in \(A\), then we have what we called in the previous chapter a DIF or distance implicit function. Though remember that we also mentioned that while using a DIF is ideal, any function that returns a distance equal to or lower than the function that would return by a DIF would work as well. Hart called such functions DUFs or distance underestimate implicit function. We can write that:

$$DUF(P) \le DIF(P)$$
The distance \(d\) returned by a DUF can be considered the radius of a sphere guaranteed not to intersect the implicit surface. Such spheres are called unbounding spheres.

Figure 2 shows the example of two points \(P\) and \(P'\) from which we want to compute ideally the nearest point on the simplicity shape (a sphere in this particular case). However for \(P\) we use a DIF to compute that distance and it returns indeed the "optimum distance" that is the distance between a \(P\) and the nearest point on the surface of the sphere, though for \(P'\) we use a DIF that returns a distance lower than that optimum distance. In summary, the distance returned by a DUF function for a given point \(P\) is equal to or lower than the distance return by a DIF function. DIFs are better than DUFs because they return the optimum distance, so with a DIF we will be able to get to the intersection point in fewer steps that with a DUF. But as mentioned in the previous chapter, it is not always possible to find a DIF for a given implicit shape therefore it is always possible in such case to use a DUF instead (assuming we can also find such function for the shape in question). The following and of course most important question is then: "how do we find either a DIF or a DUF for any given shape". Is there a formula that would allow us to find these functions systematically?

The foundation of Hart's method relies precisely on proving that indeed, it is possible using maths to find a generic solution to this problem. His method relies on a mathematical concept called the Lipschitz constant denoted \(M\) in this lesson:

$$f(a) - f(b) \le M ||a - b||.$$

Here \(a\) and \(b\) are just two values that are acceptable for the function \(f(x)\). What this equation says is that there is a constant called \(M\) such that the product of \(M\) with the absolute value of \(a-b\) is greater or equal to \(f(a) - f(b)\). \(M\) here as we already mentioned is called the Lipschitz constant and is also a simple real value like \(a\) and \(b\). It can be interpreted in this particular case as some sort of constant that scales the term \(||a-b||\) up or down. There is nothing really complicated to this idea. This theorem is only true when \(f(x)\) is continuous and differentiable (to simplify), though we won't get into these technical details nor demonstrate the theorem itself as it can be quite daunting. Suffice to say that this theorem is true for implicit functions.

Note that there is no upper limit to \(M\) but there is lower bound or minimum. In other words, there is a minimum for \(M\) such that the \(M||a-b|| \ge f(a) - f(b)\) holds true for any \(a\) and \(b\). But any value for \(M\) that is greater than that minimum will work too (for example if the minimum possible value for \(M\) is 1, then any value greater than 1 will work for \(M\) as well).

Figure 3: the Lipschitz constant \(M\) can be interpreted as the slope of the function \(f(x)\).

We can re-arrange the terms for the Lipschitz constant equation as follows:

$$ \begin{array}{l} f(a) - f(b) \le M ||a - b||,\\ M \ge \dfrac {f(a) - f(b)} {||a - b||}. \end{array} $$

This is interesting because the term on the right hand side of the second equation (second line) has a name in mathematics and also has some very important properties. It is called a difference quotient is in calculus the method by which we can compute a function derivative. The difference quotient is written as:

$$\dfrac {f(b) - f(a)} {b - a} .$$

technically this represents the mean (or average) value of the derivative of \(f\) over the interval \([a, b]\) (if you want to learn more about this topic such for the mean value theorem on the web). When \(a\) and \(b\) are far apart, we only get an approximation of the derivative of \(f\) at \(a\) (whether we compute the derivative at \(a\) or \(b\) or in the middle depends on whether we use forward, backward or central difference but we won't get into these details here). To go from the approximation to the exact solution requires to make the difference between \(a\) and \(b\) very small which we can write as:

$$f'(a) = \lim\limits_{||a-b|| \to 0} \dfrac{f(b)-f(a)}{b-a} = \lim\limits_{h \to 0} \dfrac {f(a+h) - f(a)} {h},$$

with \(h = a - b\). In words it says that the difference quotient gives the derivative (denoted \(f'\)) of the function \(f\) when the difference between \(a\) and \(b\) (which we replaced with the term \(h\)) approaches 0 (the limit).

This is great for us because in short, it also means that the Lipshitz constant \(M\) is equal to the function's \(f\) derivative:

$$M \ge \lim\limits_{||a-b|| \to 0} \dfrac{f(b)-f(a)}{b-a} \rightarrow f'(x).$$

Practically this observation is going to be useful later on to compute \(M\). Though before we get to that point we will first explain why is the Lipschitz constant useful.

Let's assume that \(b\) is a point on the surface of an implicit shape and \(a\) some random point in space (figure 4). If that is true then we can write that:

$$||b -a|| = d(a, A).$$ where \(d(a, A)\) here is our distance estimation function, the function that returns the distance from \(a\) to the nearest distance on the surface of the implicit object, which in our case is \(b\). Simple. Remember that \(d(a, A)\) is the function we are looking for. Now let's look at this other equation:

$$f(a) - f(b).$$

where \(f\) is our implicit function (such as \(x^2+y^2+z^2-r^2 = 0\) if we render spheres). If \(b\) lies exactly on the surface of the implicit object then we know that f(x) = 0. The implicit function returns the signed distance of a point to the isosurface, so if the point lies on the isosurface that distance is 0. In which case our previous equation simplifies to (assuming \(b\) is on the shape's surface):

$$f(a) - f(b) = f(a).$$

Let's put all the pieces of the puzzles together. We have:

$$f(a) - f(b) \le M||a-b||.$$

Since \(||a-b|| = d(a, A)\) we can write:

$$f(a) - f(b) \le M d(a, A).$$

and since \(f(a) - f(b)\) simplifies to \(f(a)\) when \(b\) lies on the shape's surface, which is the case for which we try to find a solution, we get:

$$f(a) \le M d(a, A).$$

We finally get a solution for our distance estimation function which is:

$$d(a, A) \ge \dfrac {f(a)} {M}.$$

The key part here is to realize that \(d(a, A)\) returns the "ideal" or optimum distance between point \(a\) and the nearest point on the surface of the object. Another way of saying this is that \(d(a,A)\) is a DIF. And recall that:

$$DIF \ge DUF.$$

Thus if the term on the right side of the equation (\(f(a)/M\)) is lower than that "optimum" distance \(d(a,A)\), then we can use that term as a DUF (if you don't understand this idea clearly read again):

$$ \begin{array}{l} DIF \ge DUF.\\ d(a, A) \ge \dfrac {f(a)} {M}, \end{array} $$

where \(DIF = d(a, A)\) and \(DUF = f(a)/M\).

In conclusion, we can find a DUF for any implicit function, if we divide the function itself by the function's Lipschitz constant, hence the need of computing a value for that constant \(M\).

So how do we compute a value for \(M\)? Remember that:

$$M \ge f'(x).$$

Now \(f'(x)\) is a function so its value is likely to vary as \(x\) varies. Though note that if \(M\) is equal to the function derivative's (\(f'(x)\)) maximum value then we will satisfy the above condition:

\(M\) is at least equal to \(f'(x)\)'s maximum value.

Figure 5: a graph of a function \(f(x)\) and its derivative f'(x). The derivates reaches its maximum value when the tangent on \(f(x)\) is the steepest.

So how do we find \(f'(x)\) maximum value? With mathematics of course. Remember from your lessons on calculus that the function derivative can be interpreted as the slope of the function \(f\). This idea is illustrated in figure 3 and 5. The slope is basically a line tangent to a point on the function where the derivative of that function is evaluated. Note that as you move from one point to the next along the x-axis, the slope (aka steepness) of the tangent varies. Sometimes it's very steep (the tangent is almost vertical), sometimes it's almost horizontal. The tangent is very steep when the function changes very rapidly from one point to the other and it's quite flat when this change is small. This variation of the slope can also be displayed as a graph. This idea is illustrated in figure 5. The red curve in figure 5 is a graph of the function \(f(x) = -3x^3 - 3x^2 + 5x + 4\). For the sake of the exercise, we will just consider the section of the curve for which \(x\) varies between roughly -1.6 and 1.2. The green curve is a graph of that function first order derivative \(f'(x)=-9x^2-6x+5\) (in the lesson we assume you know how to compute the derivative of exponential functions. You can find a lot of tutorials on this topic on the web). It shows the variation of the tangent's slope as a function of \(x\). Note that when the tangent is the steepest, the function derivative reaches a maximum. That maximum, is the value we are looking for. Alternatively, note that \(f'(x) = 0\) when the tangent is entirely flat. We now understand when is maximum value for f'(x) reached but how do we compute its actual value?

If you look at figure 5, note that the tangent at that point where \(f'(x)\) reaches its maximum value is also flat. And what did we learn about flat tangents? That we have flat tangents wherever the function's derivative is equal to 0. Thus if we wish to find where \(f'(x)\) has a flat tangent, all we need to do is compute the derivative of \(f'(x)\) i.e. \(f''(x)\) (\(f(x)\) second derivative) and then find the solution for which \(f''(x)=0\) (find the value of \(x\) for which \(f(x) = 0\)). This is illustrated in the following figure (figure 6).

Figure 6: solving \(f''(x)=0\) tells us where \(f'(x)\) has a flat tangent.

We how that \(f'(x)=-9x^2-6x+5\) so \(f''(x)=-18x-6\). This equation has a simple solution in this case which is \(x = -6/18=-1/3\). Now all we need to do in order to find \(f'(x)\)'s maximum value, is to plug our solution in the first order derivative function. We get: \(f'(-1/3) = -9(-1/3)^2-6(-1/3) + 5=-1+2+5=6\). Et voila!

This demonstration showed you a solution to find a value for the function's Lipschitz constant \(M\). We know \(M\) needs to be at least equal to the function derivative maximum value,

$$M \ge f'(x),$$

and we know that in order to compute that value, we need to find the solution to the function second order derivative (\(f''(x)\)). We then plug this solution in \(f'(x)\) and that gives us a solution for \(M\).

Hopefully you can now understand the mathematical principle upon which sphere-tracing is built. Don't worry if it still looks abstract to you, we will look at a few real examples in the rest of this lesson.

Note that we have only studied a 1D function in this part of the lesson but we deal with 3D shapes. The concept of tangent in 1D translates to the concept of vector in 3D space and gradient.

Practical Examples

Figure 7: nearest distance to point on implicit sphere.

Interestingly enough the method we just described doesn't need to be used all the times. There are some shapes for which it is in fact pretty straightforward to find the perfect DIF without making much of an effort. For example, it's pretty simple to see that the nearest distance to a sphere is just the distance from the point \(P\) to the center of the sphere (the origin for now) minus the sphere radius:

$$d(P, Sphere) = ||P|| - r =0.$$

Where \(||P||\) in code is computed with sqrtf(x*x+y*y+z*z). One important note. For now we assume that all shapes are centred around the origin (thus \(C\) the center of the sphere has coordinates <0,0,0>). Later in the lesson we will learn how to deal with transformations such as scale, rotation and translation. But for now, we will assume a canonical shape centred around the origin.

Figure 8: nearest distance to point on implicit plane.

Figure 9: nearest distance to point on implicit cube.

Finding the nearest point on a plane is simple too. We need a point on the plane (any point will do) and compute the vector that goes from that point to our position. Then we compute the dot product of this resulting vector with the plane normal. In the lesson on geometry, we explained that the dot of a normalized vector \(V_N\) (the plane's normal) with an arbitrary vector \(V_A\) that is not necessarily normalized gives the length of \(V_A\) projected onto the axis defined by \(V_N\). If you look at figure 8, you will notice that this length is just the distance \(d\) from our positon in space to the nearest point on the plane.

$$d(P, Plane) = (P - P_{plane}) \cdot N$$

where \(N\) is the plane normal and \(P_{plane}\) is a point on the plane.

Computing the nearest point to a cube can also be found using geometry. If you look at figure 9, you will see that we can encounter two different scenarios: when the point is inside the cube (\(P_3\)) when the point is outside of the cube (\(P_1\), and \(P_2\)).

Figure 10: nearest distance to point on implicit cylinder.

The solution for a torus is also quite simple. The main idea is to reduce the 3D point \(P\) to a 2D point lying in a plane passing through the torus: the intersection of the plane with the torus defines a circle (as illustrated in figure 10). By reducing the problem to a 2D point and a circle lying in the same plane, we can apply the same technique than the one we used for spheres: subtract the radius of the circle from the distance from the point to the center of the circle. This will give us the distance of the point to the circle and thus the torus.

To compute the coordinates of \(P\) as a 2D point, first compute the distance of the point from the world origin in the xz-plane and subtract to this result the torus larger radius (\(R_0\). This gives the x-coordinate of the point \(P\) on the plane with regard to the center of the circle. Then set the y-coordinate of this 2D point with \(P_y\). The length of the point gives the distance from the 2D point to the center of the circle. Subtract the inner radius of the torus to this length and you will get the distance of \(P\) to the nearest point on the torus surface. The equation is: $$d(P, Torus) = \sqrt{ (\sqrt{P.x^2 + P.z^2} - R_0^2) + P.y^2 } - R_1^2.$$

When and Where Will I Use the Lipschitz Constant Thing?

We didn't use the techniques based on the Lipschitz constant for these simple shapes, so you may think, "what was the point of studying all that if we don't use it?". In the last chapter of this lesson, we will show that when geometry can't be used to develop a DIF or a DUF, then we can always revert to the more complex but generic solution presented in this chapter.