Craters

Git repo

The git repository for my code has been made public. Anyone is free to view and learn from my code. I will try to keep the code in sync with what is posted on the blog.

Geometry

With that said, I'd like to show you how to procedurally generate craters on a planetary scale.
To find a mathematical description of what a crater should be, let's take a look at this spherical coordinate system:

As you might know, a circle C in the XY plane can be described by the equation: x² + y² = r² with r being the radius. Let's assume C is placed at a z coordinate that puts it on the surface of the sphere. If I wanted to create a crater there, I would have to loop over all the vertices, checking if they're inside C and decreasing their height.

The distance between a point and the center of the circle is defined by:

So a vertex is inside C if:

What I've described so far only works for circles centered around the z axis. If it has to work for any circle centered around an axis (φ,θ), then a change of basis is required. Because this change is only defined by two rotations, it can be easily done using quaternions. Now, explaining how quaternions work and how are they defined isn't the object of this post and there are many great articles about them on the internet. Basically, quaternions are mathematical objects that represent the transformation induced by a rotation. For example, a rotation of α around an axis u can be expressed with the following quaternion:

The rotation can be applied to a vector v by conjugating it with the quaternion:

In our case, the transformation is defined by a compound quaternion:

To obtain the coordinates of a vertex (x,y,z) in the new basis, we need to apply the inverse transformation, that is, conjugate it with the inverse quaternion.

Probability distributions

Next, we need to choose probability distributions for the radius and the angular positions of a crater. As you might imagine, these distributions differ depending on the history of the celestial body. The Moon, for example, is tidally locked to Earth, which means that it acts like a shield and most of its craters are located on its far side. Other bodies might have craters that are more or less uniformly distributed across its surface.

As for the radius, several papers have been written describing the distribution of crater sizes on the Moon and on Mars. Essentially, they showed that, the bigger a crater is, the less likely it is to exist. Instead of coding the distribution that's described in the aforementioned paper, I decided to use the gamma distribution, which comes with Python's random library:

You might think that using uniform distributions for the angles φ∈[0,π] and θ∈[0,2π] would yield a uniform distribution of craters on the planet but that's wrong. The (φ,θ) values tend to concentrate around the poles. Here's an image of a height map I was able to generate and a sphere's "wireframe":

Because of the nature of the sphere's parametrization, points near the poles are simply more likely to get randomly chosen. This problem is better explained on Wolfram's page, which also offers a solution:

where v ∈ [0,1] is a random variable.

Blending with the terrain

Finally, we need to define how to decrease the height of the vertices inside a crater. The relevant variable to use is the relative distance of a vertex to the center of the crater ρ/r so the decrease in height should be a function of this variable. I decided to simply use a polynomial (x=ρ/r):

I decided to give it an elevated rim by changing the sign of the function for high ρ/r values. In terms of code I start by finding the closest vertex to the crater's center and its height becomes a reference for the entire crater (gray horizontal line on the drawing). It means that terrain heights outside the crater isn't taken into account, which can create some very unnatural geometry:

 To fix this, the crater needs to be blended into the terrain. At the edge of the crater (ρ=r) the natural terrain height is the boundary condition, then it blends into the elevated part of the crater (ρ=0.9r). Also, the inner part of the crater takes into account a fraction of the original terrain height because real craters are not perfectly smooth.

This is a planet I generated which is supposed to look at the Moon:

I feel like the result is a little disappointing but that's a consequence of choosing colors based on heights instead of using textures. Actually, what makes a crater look like a crater, and you can see it well on the moon photo, are the lighting conditions and, in particular, the shadows inside the crater.

Improving performance

Earlier in this post I mentioned that all of the vertices were filtered to see which ones are inside of a crater. This is very inefficient and makes generating craters awfully slow. Since we know the crater center's angular position (φ,θ) and its radius r, we can calculate angular bounds within which the crater is confined. Here's a 3D illustration I made with GeoGebra:

Let γ be the circle's angular radius and β the angle between a meridian tangent to the circle and its center. Here's a beautiful figure made with sketchometry:

We can already notice that all the vertices (φ,θ) inside the circle verify:

An analogous inequality could be written for θ.
Now, γ is only a function of r but β is a function of r and φc.
β is defined for φ∈[γ,π-γ] and r < R. To find its expression it's a good idea to fix one variable and vary the other.

With a little imagination you can guess that:

Knowing γ and β, vertices can be filtered twice: first with angles and then with distance to the crater's center. It makes crater generation much faster, especially because most craters are small so angle filtering eliminates most of the circle's vertices.

Closing words

I didn't expect to spend so much time working on craters, especially because they don't bring much value to my project right now. Meanwhile, I managed define a lot of the theory needed to simulate erosion and deposition and I'm going to start writing some code soon.

Hope you enjoyed this post, stay tuned!