© ESO 2022

1. Introduction

As of 2021, international space missions have visited a total of 16 minor planets and comets. They revealed complex worlds, shaped by evolutionary processes at different scales: catastrophic collisions, regolith transport, micro-impacts gardening, and sublimation. These objects have been characterized in high resolution and form the basis of our knowledge of asteroid and comet topography. However, they represent only a tiny fraction of the population known from remote observations (i.e., more than one million asteroids; Mainzer et al. 2015).

Consequently, all scientific investigations that require knowledge of a small body shape must rely on hypothetical objects: at best, a low-resolution shape model can be derived from light curves or adaptive optics images, but often that is not the case. In the next few sections, we describe why this missing shape information is a fundamental problem in the science of small bodies, for a variety of investigations, but also for mission planning and tasks such as close proximity operations. We then present a numerical approach to constructing virtual objects at a high spatial resolution that can be used to test the influence of morphological parameters for any scientific or operational investigation that requires a shape model.

The shape of a small body, for instance, an asteroid or comet nucleus, is the combined outcome of processes that created the object and drove its subsequent evolution. These processes transform the object locally and globally in ways that are not well constrained. For instance, impact cratering is fairly well described by experiments and has been tested in situ (Deep Impact, Hayabusa 2), and it is possible to predict the resulting topography if we know the material properties of the target well enough (see Sect. 3). Other processes, such as mass wasting or spin-up deformation, are hypothesized and likely to occur, but its mechanism is not well known at present. However, by transforming the surface morphology, these processes induce effects that can be detected even if the asteroid surface features are too small to be observed from the Earth:

The Yarkovsky-O’Keefe-Radzievskii-Paddack effect (YORP; Rubincam 2000) changes the spin rate of an object due to the anisotropic re-emission of radiation, particularly in the thermal range. This effect is completely controlled by the shape of the object. If the net spin acceleration is positive, it can lead to a deformation of the body up to complete disruption. This is one of the potential mechanisms leading to the formation of binary asteroids.

The Yarkovsky effect alone uses similar arguments to explain the change in orbital parameter for small asteroids (Bottke et al. 2006). This is particularly relevant for explaining the “V-shape” dispersion pattern in collisional families (Bolin et al. 2017).

The disk-integrated photometry of an object is generally the only information available, as most asteroids or comets cannot be resolved on astronomical imaging data. However, the shape plays a role here as well. At a high phase angle, a rough topography will cast more shadows than a smooth one, introducing an excess phase-darkening that could not be explained solely with material properties (Vincent 2019).

In the coming decade, most space agencies are planning exciting missions to a series of asteroids or comets. For instance, the asteroid deflection project AIDA will impact Dimorphos, the secondary element of the Didymos binary system, using NASA’s Double Asteroid Redirection Test (DART) (Cheng et al. 2018), while ESA’s Hera will observe the aftermath (Michel et al. 2018). Together they will measure the change in the binary system orbit as well as the volume and morphology of the impact crater. Many other missions to asteroids and comets have been selected: NASA’s Lucy (Levison & Lucy Science Team 2016) and Psyche (Elkins-Tanton & Bell, 2017), JAXA’s Destiny+ (Fujimoto & Tasker 2019), CNSA’s ZhengHe (Zhang et al. 2019), and ESA’s Comet Interceptor (Snodgrass & Jones 2019). Those missions are either aimed at carrying out close orbits, landings, or sample return of the corresponding target objects. Each of these missions face the challenge that their target morphology and topography is unknown at the time of planning the first operations, and requires very flexible and powerful simulation tools for the project planning and deployment.

Having a realistic virtual asteroid or comet allows for testing mission scenarios using an object that looks real and whose distribution of surface features as well as the material properties matches the type of object we plan the mission for (at least to the best of the current state of knowledge about that object). This is fundamental for many applications such as:

Navigation (Martin et al. 2014): missions using autonomous navigation may need to test the accuracy of their targeting algorithms depending on the terrain type (many craters, smooth regolith, boulder field, etc.). At close distance, autonomous navigation requires smart algorithms that can recognize typical surface features and react in the most fuel-efficient manner.

Classification: morphological analysis requires a detailed mapping of all surface features. A virtual model of a realistic surface allows for thorough testing of feature analysis and classification algorithms, which can then help with autonomous landing site selection. This is useful for testing algorithms that may also be used in different contexts (i.e., asteroid mining, ...)

Shape reconstruction (Lamy et al. 2008): stereo mapping and shape-from-shading algorithms require multiple images from different viewing geometries and specific illumination conditions. Having a realistic model of the topography allows for us to test whether operations scenarios meet all the requirements. For instance, at a high phase angle, a smooth terrain may be illuminated, but cratered regions can have too many shadows to be mapped accurately.

In this paper, we present a new approach to synthesize realistic asteroid shapes with geospatial features such as craters. Our system models the coarse shape of an asteroid by a global implicit surface. Using a discrete global grid system (DGGS), we then create the features over the coarse shape according to a power law. We also propose novel kernels to create the morphology of the features so that it follows the observation of real images.

More precisely, the main contributions of our work are as follows: 1. We extend 3D spot noise (Li et al. 2020) by introducing specific spatial kernels to generate typical types of craters; 2. We introduce several intuitive parameters in the spatial kernels and form a degradation model that matches the observable morphology (changes) to accurately describe their degradation state; 3. We define an inversion sampling method to replicate the power-law distribution of craters on an asteroid’s surface and therefore constrain the age and evolution of a particular terrain or a surface; 4. Finally, our method provides multiple layers of control, allowing the user to author craters and tune their spatial distribution rules within each layer.

Altogether, we provide a unified framework¹ to reproduce arbitrary sized craters and allow for the reconstruction of complex surfaces and shapes of asteroids at different ages. Another advantage of our approach is that it is possible to utilize the GPU, which allows for very fast generation of the 3D shapes. Moreover, the implicit surface representation allows for arbitrarily high resolution of the synthesized meshes representing the shape.

2. Related works

Research into algorithms for generating realistic scenic terrain has a long history. The vast majority of existing techniques addresses only 2.5D heightfield scenic terrain, so that heightmaps/digital elevation maps can be applied easily on a plane. An alternative method is to use a mathematical implicit function to model the coarse terrain, and overlapping terrain features (i.e., lines, curves) or textures on the predefined coarse terrain. However, adding multiscale structures or features on the coarse terrain is a challenging task as implicit surfaces do not provide an explicit 2D parameterization. Existing solutions either rely on interactive authoring or require a parameter-free texturing method to add features through displacement mapping.

Hyper-textures (Ebert et al. 2003) as a typical parameter-free texturing method is based on procedural noise, which allows for the synthesis of a theoretically infinite amount of fractal multiscale volumetric details. This method works by perturbing an initial, smooth volumetric model or implicit surface by a fractal noise function. Although this method generates fractal surfaces with self-similar features, such as bedrock, it fails at reproducing large structures such as craters or rocks shapes. Cellular textures (Worley 1996) can be easily adapted to generate Voronoi diagrams and synthesize polygon-like patterns but are hardly controllable for the user. Another type of procedural noise is based on the convolution of a spatial filter function (kernel) with uniform random impulses (points) to generate sparse convolution noise patterns (Lewis 1987). The advantage is that spectral control can be achieved by a spectral definition of the kernel function (spatial filter). Some micro-structural features, therefore, can be produced by the sparse convolution noise. van Wijk (1991) and de Leeuw & van Liere (1997) proposed spot noise to generate structural features by replacing the kernel from spectral-oriented functions into spatial-oriented functions. More recent procedural noise functions (Galerne et al. 2012; Gilet et al. 2014; Tricard et al. 2019) rely on the power spectrum to synthesis micro-details of the example textures, but they hardly capture large structures from the power spectrum. While such procedural textures can, in theory, synthesize an infinite amount of details by accumulating them, the self-similar geometries resulting from such a fractal process cannot capture the larger structures and, more importantly, there is no guarantee to obtain geologically correct patterns observed on real terrains.

More recent works propose a hierarchical, implicit construction tree combining 3D implicit primitives (e.g., arches Becher et al. 2017, caves Paris et al. 2019, sea or rock cliffs Paris et al. 2020) to model a variety of local landforms, which can then layed over the 2.5D heightfield or procedural terrains to author scenic terrain details. While effective for modeling general large-scale landforms with mesoscale or small-scale details, most existing terrain modeling methods suffer from a lack of physical precision and fail at reproducing their effects (3D features) on more complex irregular shapes, such as asteroids.

Generating physically correct (or realistic) virtual planets or small bodies has received little attention in the computer graphics community, but it is a current topic of attention in the astronomical sciences. One possible reason to explain this gap is that the synthetization of virtual planets is a narrow field for the computer graphics community and the lack of focus on methods that allow for artistic control. However, astronomy science looks to synthesizing realistic or physically correct objects, which calls for drawing patterns or physical laws from real observations. For instance, the rough shape of a celestial body is obtained from the light curve inversion (Jorda et al. 2010). This was used to plan the flyby of asteroid 21 Lutetia by ESA’s Rosetta in 2010. Ideally, each of these missions starts with a scouting phase during which the target is mapped at high resolution and a detailed shape model is reconstructed. This forms the basis of all future operations and scientific planning. However, that is not always possible, particularly for fly-by missions that spend only a limited time in the vicinity of their target.

Therefore, it is important to plan the mission and test data analysis pipelines with different “plausible” models of the target that look as realistic as possible to plan for potential contingency. For instance, the bilobate shape of comet 67P created large shadows that severely affected the sequencing of observations, originally planned using the best guess convex shape available before the mission (Lamy et al. 2006). Specific software had to be developed in the early phase of the mission to account for this situation and synthesize “virtual” 3D models (Vincent 2018). More recently, Vincent (2019) used an multi-ellipsoids model to approximate the cometary pits to generate virtual comet nuclei at different stages of evolution. This leads to a derivation of the effects of topography on the global phase darkening reported by ground-based observers.

When generating realistic virtual small bodies, the challenges stem from the fact that asteroid and comet surfaces are produced from a host of physical processes (including impact cratering, seismic shaking, and day and night cycling erosion), which depends on the materials involved (such as silicate-rich chondrules, metal-rich matrix) and produces a variety of morphologies (from bowl-shaped depression produced by an impact to slump deposits, debris chutes or gullies) and scales (from a few meters to hundreds of kilometers).

3. Our approach

Real asteroid surfaces often feature complex disruptions of the topography, such as craters or boulder fields. These landscapes are the result of past and present geological processes, including impact melt or debris, crater parts buried by ejecta, mass wasting due to instability or erosion, all of which yield valuable information about the formation of the solar system. Usually, before a space mission, the targeted small body, S, is approximated by simple analytical expressions such as ellipsoids (see Eq. (1)), in cases where the length in three perpendicular axes can be observed from the telescope or light curve or spheres. In principle, an implicit sphere or ellipsoid defines a 2D manifold embedded in the three-dimensional space ℝ³ as:

(1)

where r represents the radius, 0 is the iso-value of the implicit surface, and a, b, c represent the axial lengths of the ellipsoid. For instance, we depict the ellipsoidal shape of the DART mission’s target asteroid Dimorphos B (Sarli et al. 2017) in Fig. 1.

Fig. 1. Rough estimate of the shape of Dimorphos (Didymos B), which is the target of the Double Asteroid Redirection Test (DART) mission (Sarli et al. 2017). The estimate of the physical size of Dimorphos is: long axis 164 m, second axis 126 m, and third axis 105 m, corresponding to the parameters a, b, c in Eq. (1). These are all the details that are known until a few minutes before the impact.

However, Tanga et al. (2009) proposed numerical simulations to explain the formation of asteroid shapes and pointed out that a significant fraction of the observed asteroids is expected to consist of fragmented bodies joined together by gravity. Therefore, real asteroids usually exhibit a much more complex structure, so that simple representations like ellipsoids are not sufficient. We chose a different approach based on our two-step AstroGen algorithm (Li et al. 2018) to synthesize the complex irregular shape of real asteroids directly, thus avoiding computationally expensive simulations such as Tanga et al. (2009). The main advantages are: 1. We use non-uniform sphere packings (Weller & Zachmann 2010) and define an implicit shape using the metaballs method (Blinn 1982) to represent the coarse shape (see Fig. 2). One benefit is that an implicit function can represent any given water-tight constraint shape. Hence, it is possible to adapt this implicit representation to any existing polygonal shape like the low-poly models from the asteroid database (McMahon 1996) (see Fig. 3). 2. Our method allows for easy creation of variants by modifying the constraint shape (see Fig. 3, Didymos A) or by varying parameters in the implicit function. This enables us to generate diverse look-alike shapes for various applications, that is, different scenarios for virtual testbed simulations. 3. Our method allows for making a direct representation of the small body’s internal mass distribution. This can then be used, for instance, to approximate the gravitational field of such irregularly shaped small bodies (Srinivas et al. 2017).

Fig. 2. Metaballs (Blinn 1982) modeling pipeline. The upper left image shows a given constraint shape, which can be an extreme low polygon mesh. The upper right image shows the result of packing (non-uniform) spheres into the constraint shape. The bottom left shows the isosurface of the implicitly defined surface, and the bottom right shows the final surface after the blending operators. This can then be turned into a high-resolution, high-detail mesh.

Fig. 3. Results of our shape synthetization system. From the constraint shape, namely, the real low-poly coarse shape of the asteroid Didymos A (Naidu et al. 2020) in the top left image, our system (AstroGen) is able to generate a series of look-alike shapes, here still without craters.

Since the initial formation of asteroids, impact cratering is identified as the most significant geomorphic process for altering their topography (Fassett et al. 2017). Our work comes from the observation that high velocity impactors collide onto the small body and form craters exhibiting a variety of morphologies (i.e., bowl-shaped craters, polygonal craters, floor-fractured craters, secondary craters, and crater chains) and diverse size-frequency distributions (SFD). To simulate the impact cratering process, we first introduce multiple kernels of a locally controlled spot noise model (Pavie et al. 2016) to synthesize different types of craters. Second, we use an inversion sampling method to reconstruct the craters’ spatial distribution on the asteroid’s surface. In addition, we propose a discrete global grid system (DGGS) to sample and polygonize the implicit function using marching cubes (Lorensen & Cline 1987). The DGGS provides an efficient method to evaluate the locally controlled spot noise model, in order to synthesize large structures such as craters.

3.1. DGGS

In order to visualize the implicit surface from our approach, we built a multilayered 3D grid in space. Each layer is a grid with a different cell size, ranging from very small up to very large cells. Each layer/grid will be used to create craters within a specific range of diameters. For convenience, we chose the cell sizes so that each cell of the coarse layer contains 2³ cells of the next more fine-grained layer. Such a partitioning of space allows us to partition the 3D shape into a multiresolution hierarchy of indexed regular cells, which we will call a “discrete global grid system” in the following. A cell here is kind of a generalization of the notion of pixel in the traditional 2D spot noise (de Leeuw & van Liere 1997). The direct benefit of this data structure is that it simplifies the hierarchical and neighbor indexing substantially. We also use it to evaluate the implicit function on each grid point and then extract the surface of the small body by using, for instance, the marching cubes method to approximate the isosurface corresponding to the surface. Furthermore, the surface can be partitioned naturally into a set of distinct patches.

Finally, we can synthesize particular patterns that follow the implicit surface using the DGGS. In order to superimpose the implicit surface with the kernel functions, we construct a local reference frame in each and every cell that is aligned with the tangent vectors at the center of the cell (Li et al. 2020) (see Fig. 4). This serves as a reference frame to perform the blending of terrain features (e.g., craters, boulders). In addition, the DGGS provides unique indices that indicate each cell’s size and geospatial location. This allows us to easily associate the particular patch of the surface with relevant properties (i.e., material types) by indexing or retrieving from a geospatial database.

Fig. 4. Visualization of our 3D spot noise pattern on an asteroid surface. We created a 3D grid that partitions the virtual asteroid (Eros, in this example) into cells, which allows us to distribute 3D noise patterns across the shape and align them to the surface.

3.2. Crater morphology

Impact cratering is an important surface process on most of the small bodies as well as terrestrial planets in the solar system. The morphology and formation mechanism of impact craters is the key process in the field of planetary science. The research of impact craters mainly focuses on two aspects: crater measurement and numerical simulations. The former studies investigate the observed real shape and spatial size-frequency distribution of craters on the surface of asteroids. They are able to provide us with vital clues about the early collisional history and physical properties of the target asteroid (see Fig. 5). The latter use impact experiments and numerical simulations to quantify the response of rocks to particular impacts. For instance, the iSALE shock physics code (Wünnemann et al. 2006) can numerically model an impact under different material properties or structures (Raducan et al. 2019).

Fig. 5. Asteroid Eros as observed by NASA’s NEAR/Shoemaker mission in 2000 (Thomas et al. 2002). The surface is covered with craters of different sizes and stages of degradation.

Using the DGGS, the surface of asteroids is partitioned by equally sized cells, and within each cell, the locally controlled spot noise (LCSN) is applied to assign craters. The LCSN model is a kind of sparse convolution noise introduced by Pavie et al. (2016). It relies on the summation of spatial kernels at uniform positions in texture space and the large structure captured in the spatial kernel can be used to synthesize large spatial structures. Therefore, the LCSN is able to synthesize bulk spatial structures that cannot be reproduced from the power spectrum alone. Subsequently, this property was extended by Li et al. (2020) into 3D space to synthesize large 3D terrain primitives, especially craters. Here, we further improve upon this method by supporting realistic and, more importantly, physically correct crater synthesis. In order to generate craters with different morphologies (i.e., normal, fresh, or complex craters), we introduce four spatial kernels to the LCSN method:

(2)

with ω ∈ [0, 1] a user-controllable weight parameter, and kernel_l is the spatial kernel that defines the shape of craters. In total, we define four different kinds of kernels, l = 1, ...4, which will represent the four different types of craters known in the literature (Mainzer et al. 2015). Also, p_i is an offset impulse following a Poisson distribution which randomly shifts the position of the kernel within each cell (Pavie et al. 2016) and i is the index of each cell. The image from the Multispectral Image Camera of the NEAR spacecraft views the surface of Eros at a resolution of 1m per pixel (see Fig. 5), which can be used as a reference for the possible morphologies of real craters (Thomas et al. 2002). It can be observed that the real asteroid surface is marked by abundant evidence of regolith, a dearth of small impact craters, but many large bowl-shaped craters with clear and distinctive flat floors. However, larger asteroids such as Vesta (Marchi et al. 2012) and Ceres (Hiesinger et al. 2016), observed by the Framing Camera onboard NASA’s Dawn spacecraft, reveal many diverse surfaces with a large number and variety of crater morphologies, including bowl-shaped craters, polygonal craters, floor-fractured craters, crater terraces, central peaks, smooth crater floors, craters with lobate flow features, and secondary craters and crater chains (Hiesinger et al. 2016). Simple craters range in size from centimeters to tens of kilometers; the maximum size is inversely related to the planetary mass (Dombard et al. 2010). Simple craters are characterized by a consistently concave-upward shape with a nearly parabolic or conical interior profile. Craters larger than 40km usually exhibit special morphologies not observed at smaller craters, such as central pits, possible pit floors, and floor fractures (Hiesinger et al. 2016).

Because of a rather large simple-to-complex transition diameter, most small asteroids display simple craters, namely, bowl-shaped depressions produced by an impact. Simple craters can exhibit a bench-like step in the walls or a central mound if formed on a surface with layers of varying strength (Hiesinger et al. 2016). Otherwise, they show few other internal topographic features, with the exception of slump deposits, debris chutes, gullies, or sporadic boulder trails (Hiesinger et al. 2016). Therefore, we first introduce the following two kernels to represent the simple bowl-shape normal craters:

(3)

where kernel₁ is the standard Gaussian function (q₁ = 1, see Fig. 6, blue line) and σ represents the standard deviation (, ). We can define two parameters H_t, D_a to represent the height (true height) and the apparent crater diameter of the simple bowl-shape crater. The parameter q can be used to tune the smoothness of the blending of the boundary of the crater with the terrain and σ controls the length of the flat floor in the simple crater. Another type of simple crater may have a central flat floor of crater-fill deposit that originates from mass wasting and eolian or fluvial infilling or ponding of impact melt (Dombard et al. 2010). To create those, we set q₂ = 2 and obtain kernel₂ to represent this type of simple crater (see Fig. 6, red line). We therefore evaluate the height and diameter of two kernels (see Eqs. (4) and (5)) and this evaluation is motivated by the full width at half maximum (FWHM) of Gaussians where (Bovik et al. 1990).

(4)

(5)

Fig. 6. Profile lines of all kernels used to generate the different morphologies of different types of craters. The simple crater 1 (blue line), simple crater 2 (red line), fresh crater (green line), and complex crater (yellow line) correspond to kernel1, 2, 3, 4, resp. The kernel1 is a standard Gaussian function with the true height Ht = 1, the apparent diameter . In kernel2, Ht = 1 and (see Eqs. (4) and (5)). The fresh crater is created using kernel3 with Ht = 1,  Da = 1.31 and extra parameters Drim ≈ 2,  Hrim ≈ 0.309, which represent the rim-to-rim diameter and the rim height, respectively (see Eqs. (7)–(11)). The new degradation length, the rim-to-rim diameter and the rim height of the complex crater are defined as Drim = 2,  Hrim = 0.471,  Hc = 0.2 (see Eqs. (13)–(16)).

Another interesting characteristic of craters on asteroids, distinct from other terrestrial planets, is that they appear deprived of impact melts due to the lower average impact speeds. Therefore, they better show the morphology of the floor (kernel₂ contains the central floor) and walls (Marchi et al. 2015). Consequently, another type of crater we are aiming at is typically bowl-shaped but with a small rim (bulge effects on the edge, named walls) called “fresh crater”. The height of the rim (walls) depends strongly on the material properties and internal structure of the subsurface (Marchi et al. 2015). Fresh craters are polygonal to circular bowl-shaped depressions with smooth to ragged raised rims and surrounded by a layer of continuous ejecta that overlies the uplifted crater rim and flanks, usually extending about one crater diameter from the rim (Hiesinger et al. 2016). Hence, we specially designed kernel₃ to generate younger fresh craters with rim and flanks:

(6)

where κ and b control the height above the 0-plane; β and γ define the radii of the rim above the 0-plane; and a controls the depth of the crater under the 0-plane.

Similarly, we can define five variables H_rim,  H_t,  D_rim, D_a, L_d to represent the accurate morphology of a fresh bowl-shape crater (see Fig. 6, green line). Since the parameters of kernel₃ are controlled by our method, we let kernel₃ = 0 and in that way, we obtain the crater morphology (W is the Lambert W function; Corless et al. 1996):

(7)

(8)

(9)

(10)

(11)

Often we can simplify the evaluation by setting κ = a.

Simple craters mostly lack the central uplift, peak rings, or terraced wall of larger, complex craters. The transition diameter from simple craters to complex craters, termed the “simple-to-complex transition”, varies according to different planetary bodies. In order to synthesize such large complex craters with central uplift, we propose kernel₄ to be:

(12)

where κ and b control the height above the 0-plane, β and γ define the radii of bulge above the 0-plane, a controls the depth of the crater under the 0-plane. From kernel₄, we can derive four variables, H_rim,  H_c,  D_rim,  L_d, to accurately define the morphology of the complex crater with central uplift (see Fig. 6, yellow line). In the following, we set a threshold kernel₄ = 0.01 as the border of the complex crater to evaluate the parameter L_d:

(13)

(14)

(15)

(16)

The definition of our four kernels enables us to accurately control the morphology of different crater types. The easy-to-use parameters may help to study the impact mechanisms, impact effects, or even the evolutionary history of craters. For instance, the numerical simulation of impact cratering on different materials (e.g., granite, quartzite) shows a similar curve in peak pressure, but different densities will result in subtle differences in the crater morphology (Raducan et al. 2019).

We procedurally generate different types of craters as implicit primitives (kernels) on the implicit surface (see Fig. 7). The algorithm can easily generate more shapes; the user has control over what is most relevant to his or her scenario (e.g., only simple craters on small asteroids).

Fig. 7. Examples of different crater types that can be generated on a small body. The uppermost image shows the simple crater (type 1 and 2) using kernel1, 2. The middle image represents fresh craters generated from the kernel3, and the bottom image depicts the morphology of the complex crater from kernel4.

3.3. Crater degradation model

The degradation of impact craters on the surface of airless bodies under the cumulative effects has been a subject of interest for a long time. It can be used to estimate the age of the surface and is very useful for assessing the evolution of small body surfaces. In general, the research of craters is focused on dealing with simple craters and regardless of erosion or collapse of the border. Pelletier (2008) proposed a diffusive and theoretical model of the degradation process to predict the interior or exterior slopes of craters at a given size. The basic assumption is: the older the surface, the more advanced the degradation of the craters. This model has been frequently applied to simulate the retention age of terrestrial hillslopes, namely, the time elapsed since its formation or global surface reset.

More recently, Fassett & Thomson (2014) developed a more accurate diffusive crater degradation model which can be compared directly to the real profile of craters on Mercury. For this purpose, first, they extracted the real crater profile from the digital terrain models obtained from the space mission. Then, they built a numerical model to fit the real crater profiles, which enables the estimation of various types of a crater’s degradation state (see Fig. 11). Craters from the diffusive model display a wide range of shapes, which can be translated into a wide range of crater depth-to-diameter ratios. Based on Marchi et al. (2015), we usually expect the variation of depth with respect to diameter to follow a near linear power law, depth = ratio ⋅ Dⁿ, where n ≈ 1 and ratio is a constant for a given terrain.

To demonstrate that our model can describe morphologies such as those expected from erosion physics, we present a degradation model derived from Lunar studies, for example. It can be applied to other terrains covered in loose regolith, for instance, on large asteroids. However, this would not apply to comets, which experience very different types of erosion, mostly driven by sublimation – and this is an ongoing topic of research.

Degradation models, such as the one proposed by Fassett & Thomson (2014), can be easily integrated into our crater degradation model: by mapping the variations of all over the surface (in our case, mapping to each cell of the surface), we can bring the constraint on the degradation status of simple craters, which reflect the relative age and surface properties of different regions. So we can easily define the depth-to-diameter ratio for simple craters 1  &  2 as (see Eqs. (4) and (5)):

An advantage of our method is that we can assign random range values to the variable a to constrain the degradation state all over the surface. For instance, in Fig. 8 the solid and dashed simple craters 1 and 2 were obtained by different depth-to-diameter ratios to simulate the degradation state (also see the degradation model in Fig. 11 on the simple crater with diameter D = 300 m). Accordingly, our degradation model for the simple crater can be defined as: , where Δa represents a small change of the variable a. Therefore, we can tune the variable a to a specific value that corresponds to the crater’s degradation state.

Fig. 8. Profile lines showing crater degradation. We can model the degradation status of craters by changing the parameters in our kernel functions. In simple craters 1 & 2, the solid line represents the original crater, while the dashed line represents the same crater but with a lower (degraded) depth-to-diameter ratio (ratio of solid/dashed line = 0.4/0.1). In fresh craters, we set the degradation vector α = (Δk = −0.3, Δγ = −0.4, Δa = −0.25) (defined in Sect. 3.3). For complex craters, we set the degradation vector α = (Δk = −0.05, Δγ = −0.2, Δa = −0.1).

The degradation of fresh craters and complex craters can be divided into different parts (also see the degradation model in Fig. 11 with craters of diameter D = 1 km and 3 km, respectively). For the fresh crater (modeled by kernel₃), the rim height, H_rim, will degrade as , and the interior part will also degrade. In other words, the diameter of the rim, D_rim, will increase as , the crater depth, H_t, will become shallower as , and the wing of the crater, L_d, will blur with the decrease of the rim height. Therefore, we can define the degradation model of fresh craters as: . In case of complex craters (kernel₄), we can similarly define the rim height, H_rim. Accordingly, the central uplift area, H_c, will decrease and the rim diameter, D_rim, will increase. Finally, to model the degradation of fresh and complex craters, we can build a degradation vector α = (Δa, Δκ, Δγ) to change parameters in the corresponding kernel functions (kernel₃, kernel₄). In Fig. 8, we show the effect of our degradation models applied to each kernel (simple, fresh, or complex craters).

4. Evolution of cratered asteroid surface

Cratered asteroid surfaces provide time-integrated snapshots of the cumulation of impacts over a certain time frame. Therefore, the cumulative number of craters superposed on a given terrain should monotonically increase over time. Thus, the number of observed craters (ideally) can be interpreted as evidence to bound the age of the terrain. In general, a higher density of shallow craters mostly involves much older surfaces, as craters usually get shallower or even get wiped out with age due to erosion, infilling, or seismic shaking. However, under different circumstances, craters can evolve in different ways. For instance, two identical craters of similar age and size have different depths on different terrain types (i.e., materials). Therefore, making a one-to-one correlation between fresh-eroded and young-old is not always accurate. As a result, the temporal evolution of cratered terrain is rather complex and, its correct interpretation involves considering various special mechanisms. However, in theory, the production of craters is primarily due to the flux of incoming impactors and is susceptible to temporal variations reflecting variation in the impact rates or impactor SFD (or both). Crater degradation is a natural result of the topography over time. For the sake of simplicity, we only consider the natural degradation of craters and the monotonic increase of craters on the surface to deduce the age, and regardless of complex physical events such as impact erosion, filling, or seismic shaking. Here, we use prior research on the mechanism of crater degradation numerically, with the aim to derive a correlation between the surface evolution and the parameters in our kernel functions.

4.1. Crater size-frequency distribution

Once the coarse shape of the asteroid has been created, craters of different sizes can be assigned to each cell. We use inversion sampling to simulate the crater’s size-frequency distribution based on the assumption of the equal possibility of collision events within each cell of the surface (see Fig. 4). The SFD of craters can serve to invert the distribution of impactor’s sizes, densities, and velocities, and, furthermore, reconstruct the target’s properties such as surface gravity and strength (Marchi et al. 2015). The high resolution images from spacecraft not only manifest the real morphology of craters, especially their sizes, but also reveal another fundamental property of a group of craters: the ratio between the number of small and large craters, namely, the size-frequency distribution (SFD; Marchi et al. 2015). The crater’s SFD is commonly approximated as a power law, N(D) = k ⋅ D^−B, where N(D) is the cumulative number of craters greater than or equal to a particular diameter D, and k, B are empirically determined (Neukum et al. 2001).

There exist various mechanisms that will affect the measurement of the SFD at a given surface. For instance, the ejecta from the impactor will form secondary craters and occasionally erase existing craters or alter their morphology (i.e., mass wasting). Moreover, SFDs are also intensely affected by the physical properties of the target surface (i.e., hard rock, rubbly megaregolith, and icy or volatile-rich material) (Marchi et al. 2015). Therefore, the slope of an SFD usually varies spatially across the target asteroid’s surface, suggesting lower (lower slope) or higher (higher slope) energy-density events (Belton et al. 1994). As a result, the synthesis of different SFDs may eventually help crater studies to expose many properties of asteroid interiors, surfaces, or even the history of geological processes.

The early observation of König et al. (1977) and subsequent laboratory experiments suggest that the fragments from the collisional breakup or impact cratering events have more complex SFD curves than a single power law (Fassett et al. 2017; Strom et al. 2005). Therefore, the SFD curve can be fitted either by segments of individual power laws or more complex polynomial curves instead of a straight line in the logarithmic coordinate system. For sake of simplicity, we will still use the single power-law function, namely, N(D) = k ⋅ D^−B, as the crater production function instead of the accurate global model by Neukum et al. (2001), but more complex SFD’s can be integrated into our simulation just as well.

In our asteroid model, we have subdivided the surface into equally sized cells (see Fig. 4) and within each cell the upper bound of the size of the craters is given by s. Therefore, in our method the total number of craters is determined by the number of cells and the parameter s. We partition the total range of crater diameters for a given asteroid (based on the scaling factor s) into a number of sub-ranges; within each range, we generate a number of craters, as specified by the production function (power-law). For instance, in Fig. 9 we assume that on the given asteroid, Didymos B, there exists one large crater with a diameter of 100 m (this basically determines parameter s). We also defined three bins of crater sizes [10 m, 30 m], [30 m, 70 m], [70 m, 100 m]. For each bin, we then evaluated the corresponding cumulative number of craters corresponding to the cumulative power-law distribution (B = 3).

Fig. 9. Example of crater populations generated for an asteroid of size comparable to Dimorphos. We constrain the distribution by setting exactly one crater at the largest size (here 100 m). We then estimate the number of craters at smaller sizes with a cumulative power-law distribution (dotted line). Discrete points show several crater populations generated from our inversion sampling method. Deviations from the ideal power law are due to the nature of fitting a discrete population to a continuous cumulative function. The colored area represents the 95% confidence interval.

In order to reconstruct the power-law distribution from the inversion sampling (Devroye 1986), we need to compute the probability density pdf(D) from the cumulative power-law distribution (Clauset et al. 2009), where:

(17)

This gives the probability density distribution as pdf(D) = c ⋅ D^{−B − 1}, where c is a constant. In each crater size range [D_L, D_R], we can do inversion sampling on the variable D ∈ [D_L, D_R] and the variable D follows the power law distribution (D ∼ cdf⁻¹(u), where u∼ uniform distribution:

(18)

so

(19)

(20)

(21)

where k is a constant, B is the power law index, D_L, D_R represent the lower and upper bound of the crater size, and u is the uniform distribution. Overall, we can thus generate corresponding diameters of craters that follow the power law distribution from the uniform distribution u.

For the simple crater, the diameter (see Eq. (5)), and the power law distribution of the diameter D follow (see Eq. (21)), which can be rewritten with the parameter a as . As a result, we can use the kernel parameter a to generate the power law distribution of craters within each crater size range.

As mentioned above, a single power-law approximation hardly matches the SFD curve seen on the surface of a real small body. To accomodate this, we can simply assign different power-law indices to different crater size ranges to approximate a more complex SFD curve (consist into a more complex piece-wise SFD curve).

Burchell & Leliwa-Kopystynski (2010) showed that, for small solar system bodies (in a size range of 1–500 km), the ratio of the largest crater’s diameter, D, to their mean radius, R, is a nearly constant value: D/R = 0.9 ± 0.05, averaged over rocky and icy bodies. However, whilst a shape model may exist for the body, it is difficult to determine its mass, so the overall porosity can hardly be obtained. From a group of small rocky bodies with clear densities and porosities, we can observe the relationship between the ratio D/R and the total porosity which could be interpreted using a weighted linear model (Burchell & Leliwa-Kopystynski 2010). Therefore, the maximum size of a crater on a particular small body will “indirectly” ascertain the porosity of the parent body. Usually, the size of the largest crater can be determined by a camera very well, but the huge number of smaller craters cannot be observed well given the limited resolution of cameras. Our method can close this gap very well by setting the size of the largest crater, then synthesizing a massive number of smaller craters on the small body’s shape from different power laws through inversion sampling (see Fig. 10).

Fig. 10. Application of the size-frequency distribution (SFD) defined in Fig. 9 to the asteroid Dimorphos. Both images have been generated using the same crater production function (cumulative power-law distribution with B = 3). Due to randomization, we can generate an infinite number of similar “look-alike” shapes.

Our fitting of the power-law was done in log space. Counting the precise number of craters is, up to now, difficult, because usually the surface is covered with boulders, rocks, and dust. Also, many craters are very shallow and hard to distinguish. Overall, the difference on the final SFD curve is not that huge (no more than 5%).

4.2. Crater aging

Marchi et al. (2015) show that the simple crater’s depth-to-diameter ratio follows a near-linear law, , where k for asteroids is approximately in [0.2, 0.05] and the variation reflects the different states of degradation. The newly formed simple craters usually own higher depth-to-diameter ratios (in [0.2, 0.25]), then suffer erosion and evolve into much shallower craters (k ≈ 0.05). Small craters have typically a lower depth-to-diameter ratio. Similar morphometric relationships exist for the type of fresh craters where the decay of their ejecta forms the rim at the border. The fresh crater with sharper rims usually contains a higher depth-to-diameter ratio (0.14) than small craters. These relationships exist because most fresh, primary, simple craters in the size range considered here have consistent morphologies.

Furthermore, Marchi et al. (2015) find the overall distribution of ratio usually peaks near a particular value, and more complete statistics reveal that the average varies from region to region. For instance, on Vesta a double peak value (at 0.15 and 0.19) is observed in the global distribution of with about 25% of difference between peaks (Vincent et al. 2014). The same effect can be seen on Lutetia, with up to 40% of difference in the average ratio between different regions (Vincent et al. 2012). The morphology of craters on Eros is measured through the NEAR Laser Rangefinder data, which indicates that many of the freshest craters typically possess near 0.2 for D > 100 m. For D < 100 m, fresh craters on Eros are shallower than those on the Moon. The average of all craters on Eros indicates that most are degraded and follow a pattern akin to what is observed from Vesta, Steins, and Lutetia. This variation means that different surface regions suffer different collisions and possess non-homogeneous physical properties. Here, we assume the physical property of an asteroid will remain stable and the modification of the morphology of a crater is only due to its degradation rate.

Consequently, we can design a function to link the degradation state with the crater’s age and, furthermore, link the ratio with age. For instance, the dispersion of for small craters on regions of similar physical properties on Vesta implies a degradation rate of ∼3 × 10⁻⁷ m ⋅ yr⁻¹ (Vincent et al. 2014). Schmedemann et al. (2020) estimate an erosion rate of 6 × 10⁻⁸ to 3 × 10⁻⁶ m ⋅ yr⁻¹ for kilometer-sized asteroids. We can estimate that over each Gyr, the crater’s diameter will erode about 3km (i.e., its depth will decrease) at the rate of ∼3 × 10⁻⁷ m ⋅ yr⁻¹. We can easily approximate this rate with our degradation vector α: each Gyr, α = (Δk = −0.1k, Δγ = −0.1r, Δa = −0.1a) (see Sect. 3.3, degradation model). We also estimate the degradation rate from the degradation model of Fassett & Thomson (2014) in Table A.1. Their model simulates morphometric variation of newly formed craters like fresh craters from real measurements, and we approximate their results with our kernel₃. Figure 11 depicts the morphometry of simple, fresh craters of specific diameter (from 300m to 3km) with our kernels, which allowed for a database of crater profiles in various states of diffusive evolution to be established.

Fig. 11. Crater aging: We applied the crater degradation model of Fassett & Thomson (2014) to three craters of different diameters (D = 300 m, 1 km, 3 km) from an age of 0 to 3 Gyr. Each line in those figures corresponds to a certain age interval ∼1 Gyr (0 Gyr, 1.0 Gyr, 2.0 Gyr, and 3.0 Gyr). After 3 Gyr of evolution, the smallest crater is almost indiscernible, the 1 km crater is flattened but still has distinct relief, and the 3 km crater only shows apparent changes near its rim and interior floor.

In conclusion, we can use the ratio of depth-to-diameter () and the degradation vector(α) to control the morphology of simple craters, complex, and fresh craters. From those, we can derive a spatial distribution function (crater SFD) to manipulate the density of craters of different sizes on the surface. In that way, we can synthesize an asteroid surface which corresponds to different ages. More specifically, we can easily synthesize “young” (see Fig. 12, low-density crater surface) and “old” (see Fig. 13, i.e., high-density crater surface) asteroid surfaces, respectively.

Fig. 12. Simulated Dimorphos-like asteroid, with low crater density. We add several layers of craters on the surface, e.g., simple crater types 1 and 2, and fresh craters. The total saturation is nearly 20% and the slope of power law is B = 3. The depth-to-diameter, , ratio of simple craters is in the range [0.15, 0.25].

Fig. 13. Simulated Dimorphos-like asteroid, with high crater density. As in Fig. 12, we add several layers of craters (simple craters 1 and 2, fresh craters). The total saturation is nearly 30% and the slope of the power-law is B = 4. The depth-to-diameter ratio, , of simple crater is in [0.12, 0.19]. In fresh craters, we set the degradation vector as: α = (Δk = −0.15, Δγ = −0.1, Δa = −0.25).

We would like to point out that we also added a dust details layer on the final result in order to improve the appearance. Therefore, some of the very small craters could disappear under the dust layer. Since modeling dust is not the focus of this paper, we just used conventional lattice-based random noise, as described in Li et al. (2018). We do not that other, more complex noise models such as Gabor noise, could produce more realistic dust patterns.

5. Results

We implemented our asteroid modeling system in C++ using OpenGL and GLSL. All examples in this paper were created on a desktop computer equipped with Intel 3.8 GHz Core i7 with 64GB RAM and an NVIDIA TITAN V. In order to generate a watertight mesh, we extracted the isovalue points from the 3D grid and used a screened Poisson surface reconstruction algorithm (Kazhdan & Hoppe 2013) to generate the final mesh from the extracted point cloud. The output mesh was then loaded into Blender to render photo-realistic asteroid surfaces using ray-tracing.

Figures A.1–A.4, 14, and 15 show example results of our method for synthetic shapes tuned to resemble Eros, Bennu, and Didymos A. In these examples, we have set parameters such that the diameter of the craters varies between 10–100 m (following the power-law). When accumulating the layers, the weights should be chosen carefully to avoid unwanted morphologies such as overhangs.

Fig. 14. Simulated Didymos-A-like asteroid, with low crater density. We added several layers of craters on the surface with different types (i.e., simple crater types 1 & 2, complex, and fresh craters). The slope of power-law is B = 3. The depth-to-diameter ratio, , of simple craters and fresh craters is in [0.1, 0.2].

Fig. 15. Simulated Didymos-A-like asteroid, with high crater density. We add several layers of craters, like simple crater 1 & 2, complex and fresh crater, on the surface. The slope of power law is B = 3. The depth-to-diameter ratio, , of simple craters is in [0.08, 0.16], and for fresh craters is in [0.12, 0.19]. For fresh craters, we set the degradation vector α = (Δk = −0.2, Δγ = −0.2, Δa = −0.38). For complex craters, we set the degradation vector α = (Δk = −0.04, Δγ = −0.4, Δa = −0.2).

The 3D meshes generated with our algorithm are well suited for scientific analysis: we can generate models with specific sets of parameters to describe the surface under different test hypotheses, for example, by selecting particular crater size distributions that are representative of what we expect for a young or evolved surface. Once the mesh is generated, we can export it in standard formats (e.g., OBJ, PLY) that can be processed by other analysis tools. To illustrate such an approach, we show in Fig. 16 the distribution of shadows on an Eros-like virtual asteroid for different phase angles. The shadow fraction in each image is an important parameter when trying to understand the effect of large-scale roughness on the thermal evolution of asteroids, and must be accounted for in the photometric analysis (Vincent 2019). Figure 17 shows the modeled effective gravity (gravitational acceleration plus centrifugal force on the surface) calculated with the approach described in Cheng et al. (2012) and 433 Eros’ physical parameters, as well as the gravitational slopes on the surface. We can use such calculations to evaluate the effects of regolith migration through local ejections and mass wasting. Another application of our method is training and testing of algorithms for navigation on spacecraft: using different synthesized asteroid surfaces, we can provide realistic synthetic sensor input, which correctly simulates shadowing of craters depending on a given solar elevation, in order to train and test the robustness and resilience of the navigation algorithms.

Fig. 16. Simulated Eros-like asteroid, observed at different phase angles (solar elevation). The shading is Lambertian and does not attempt to describe complex processes such as the opposition effect, but rather aims to display the extent of shadowed areas as a function of the phase.

Fig. 17. Effective gravity and gravitational slopes on a simulated Eros-like asteroid.

We also investigated the performance of our method with respect to the number of layers, namely, the number of different ranges of crater sizes. Figure 18 shows the mean average computation time for asteroids with respect to the particular layers of craters. Computation times include all steps of our method, including converting the implicit surface into a mesh. (Here, the resolution of the meshing was set to around 200 000 polygons.) The time increases almost linearly with an increasing number of layers; it varies a great deal among different asteroids, which is mostly due to the different number of spheres needed in the metaballs representation. The creation of craters can be parallelized on the GPU thanks to the partitioning with the 3D grid (one thread per cell). Obviously, the analytical expression of the ellipsoid allows for much better performance than the metaballs.

Fig. 18. Computational time of our algorithm to generate different asteroid models with different layers (crater bins). For each asteroid (Ellisoid, Eros, Bennu) we tested 1, …, 5 layers. The bulk of the time is spent on extracting the polygonal mesh.

Obviously, these computation times are two high for interactive experimentation with the parameters of our method. On the one hand, however, for interactive usage, we think that rendering by simple ray casting will be feasible and much faster. On the other hand, extracting meshes is useful for other purposes, once a good set of parameters has been found, such as importing the shape into virtual testbeds.

6. Conclusions

In this paper, we present a novel technique to synthesize realistic 3D shapes of small celestial bodies such as asteroids. Such realistic surface features as craters are defined by kernel functions and their spatial distribution is governed by laws well-known from the literature, such as the power law. Those laws are initially parameterized by results from the extensive literature, mission data, and existing numerical simulations, but can easily be modified by the user to be adapted to specific use cases.

We can generate shapes of asteroids and some types of comets that are suited to the planning and preparatory data analysis of future missions, for which the shape and morphology of the target are essentially unknown until arrival. For instance, the upcoming Hera mission (Michel et al. 2018) will measure the morphology of the impact crater and trajectory deflection produced by the DART mission (Sarli et al. 2017) and assess the hazard and effect of the impact on the binary asteroid system. They could use our tool to test operational scenarios and data analysis techniques under realistic conditions (see Fig. 19). In addition, the integrated photometric signal of a surface or the thermal inertia are parameters strongly affected by the roughness of the surface and may require different types of orbit for a spacecraft to acquire conclusive data (i.e., images of the terminator orbit could be compromised by too many shadows). Our model can simulate surfaces with different types of crater shapes and size distributions, generating meshes that can be used to test operations under those different conditions. Although this paper focuses on craters as the major type of topographic features, our method is also suitable for generating rubble pile asteroids that are covered with boulders of various sizes, which could be described by a different group of kernel functions; this will be implemented in a future version of our tool.

Fig. 19. Simulated views of asteroids Didymos and Dimorphos, as they will be observed by the Hera mission (ESA). Left panel: using the object’s shape according to current best knowledge, obtained from radar observations (Naidu et al. 2020). Right panel: same models, upscaled to high resolution topography with the technique described in this paper.

Of course, it is possible to perform such analyzes using shape models derived and re-scaled from previous missions. But that approach would not allow for blind testing of potential scenarios in which the surface may radically differ from expectations. Thus, it seems critical to be able to rapidly generate new shape models that allow for testing varieties of observed conditions.

Our discrete global grid system (DGGS) is able to digitize the small body and provide a framework to visualize, analyze and combine various types of datasets of given small bodies. For instance, geospatial imagery datasets (such as aerial or satellite photographs) could be used to texture the surface.

Our method opens up a number of avenues for future works. One main avenue is to improve the method itself, for instance, by introducing other crater and erosion models applicable to comets, by simulating more erosion factors (such as thermal erosion or boulder mass wasting), or by adding more realistic noise functions to model dust on the surface.

Another avenue for future work are a number of potential applications of our method and tool, which would enable future studies of asteroids and comets that require a high-resolution shape model. Such studies are beyond the scope of this paper, but here we introduce a couple of examples that highlight direct potential applications in asteroid and comet science:

• Disk-integrated photometry and phase function: We mentioned in the introduction that the phase darkening of comets can be partly attributed to topographic roughness. Vincent et al. (2017) and Vincent (2019) have shown that this topography is a function of the dynamic age of a comet, and can be parameterized with a power law. The study was conducted by generating virtual comets using an algorithm that is akin to a simplified version of what we presented in this manuscript. For each virtual object, we generated a phase curve using ray-tracing modeling to accurately account for all the shadows. This effect has not yet been studied on asteroids, but we expect that it might also play a role and could be used to better assess the evolutionary state of the objects (e.g., how eroded the craters are). Figure 16 shows an example of a virtual asteroid generated with our model, using the phase- or ray-tracing modeling introduced in Vincent (2019).

• Gravity and slopes: From the shape and internal density, one can easily calculate the gravitational field on the surface of a small body, using for instance the algorithm by Cheng et al. (2012). Figure 17 shows how the slope distribution would look like on an Eros-like asteroid generated with our tool. The angular distribution of slopes is a characteristic measure of the roughness and can be used to predict areas more likely to experience mass-wasting or ponding. It is also a necessary parameter to assess the safety of areas under consideration for landing or sampling. Using a virtual asteroid allows operations experts to prepare for any possible scenario.

• Thermal modeling: High-resolution topography and multiscale roughness are necessary for precisely modeling the distribution and evolution of temperatures on an asteroid. Figure 19 shows that pre-encounter shape models provide a very low-resolution description of the topography at best, which is not suited for thermal studies. Features added by the algorithm presented in this paper provide a realistic surface that allows in-depth investigation of the asteroid properties. We can, for instance, apply different crater size distributions and morphology to test the influence of topography on the regional thermal inertia and to model the signal that would be observed from ground- or space-based instruments.

Another issue is boulders, which ought to be generated with a method similar to the one we have presented here. It remains, however, an open question as to whether our approach could yield a similar (easy) method, since boulders are formed by rather different astrophysical processes. Some boulders are formed by impacts and are subsequently distributed in the ejecta blanket (Enos & Lauretta 2019). Others may have been created when the parent asteroid breaks into pieces and re-accrets (a common process for small asteroids) (Hergenrother et al. 2019). The spatial distribution of these boulders must account for the complex gravitational field of the asteroid and require much more complex modeling than what we typically do in the case of craters. Even after formation, boulders will continue to fragment and “migrate” due to thermal stresses and gravity (Hergenrother et al. 2019). These effects should also be accounted for, going beyond what our current method for craters can achieve. Recently, several works of research from the missions to Ryugu and Bennu revealed that the boulders are moving on the surface (with low gravity) and we plan to find new methods to simulate those effects.

Acknowledgments

This work was partially funded by DLR grant 50NA1916. J.B.V. acknowledges funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 870377 (project NEO-MAPP). X.L. acknowledges funding from the China National Space Administration’s Pre-research Project on Civil Aerospace Technologies (Grant No. D020303). This research has made use of NASA’s Astrophysics Data System and the scientific software shapeViewer (www.comet-toolbox.com).

References

Appendix A: More results for different types of asteroids

In this appendix, we provide more results from our algorithm. More specifically, we have applied our algorithm to two big families of asteroids: the ones similar to Eros and those similar to Bennu. Figures A.1, A.2 show simulations of the ”young” or ”old” Eros-like asteroid surface, while Figs. A.3, A.4 show the ”young” or ”old” Bennu-like asteroid surfaces. In all cases, we added several layers of craters on the surface, consisting of the types simple crater 1 & 2, complex crater, and fresh crater. We also list the kernel parameters in Table A.1 (kernel in Fig. 11) to help us inverse the crater’s age by using our kernel parameters.

Fig. A.1. Simulated Eros-like asteroid, with low crater density. Just as in Figs. 14 and 14, we added several layers of craters (the surface partitioned into nearly 1500 patches). The total saturation is nearly 11% and the slope of power law B = 3. For both simple and fresh craters, .

Fig. A.2. Simulated Eros-like asteroid, with high crater density (the surface partitioned into nearly 1500 patches). The total saturation is nearly 18% and the power law index B = 3. For simple craters, , while for fresh craters . For fresh craters, we set the degradation vector α = (Δk = −0.15, Δγ = −0.2, Δa = −0.38). In complex crater, we set the degradation vector α = (Δk = −0.2, Δγ = −0.3, Δa = −0.25).

Fig. A.3. Simulated Bennu-like asteroid, with low crater density. The total saturation is nearly 19% (the surface partitioned into nearly 2000 patches) and the slope of power law B = 3. For both simple and fresh craters, .

Fig. A.4. Simulated Bennu-like asteroid, with high crater density. The total saturation is nearly 25% (the surface partitioned into nearly 2000 patches) and the slope of power law B = 3. For simple craters, we have , while for fresh craters . For fresh craters, we set the degradation vector as α = (Δk = −0.2, Δγ = −0.2, Δa = −0.38). In complex crater, we set the degradation vector as α = (Δk = −0.04, Δγ = −0.4, Δa = −0.2).

All Tables

All Figures

	Fig. 1. Rough estimate of the shape of Dimorphos (Didymos B), which is the target of the Double Asteroid Redirection Test (DART) mission (Sarli et al. 2017). The estimate of the physical size of Dimorphos is: long axis 164 m, second axis 126 m, and third axis 105 m, corresponding to the parameters a, b, c in Eq. (1). These are all the details that are known until a few minutes before the impact.
In the text

Fig. 2. Metaballs (Blinn 1982) modeling pipeline. The upper left image shows a given constraint shape, which can be an extreme low polygon mesh. The upper right image shows the result of packing (non-uniform) spheres into the constraint shape. The bottom left shows the isosurface of the implicitly defined surface, and the bottom right shows the final surface after the blending operators. This can then be turned into a high-resolution, high-detail mesh.

In the text

	Fig. 3. Results of our shape synthetization system. From the constraint shape, namely, the real low-poly coarse shape of the asteroid Didymos A (Naidu et al. 2020) in the top left image, our system (AstroGen) is able to generate a series of look-alike shapes, here still without craters.
In the text

	Fig. 4. Visualization of our 3D spot noise pattern on an asteroid surface. We created a 3D grid that partitions the virtual asteroid (Eros, in this example) into cells, which allows us to distribute 3D noise patterns across the shape and align them to the surface.
In the text

	Fig. 5. Asteroid Eros as observed by NASA’s NEAR/Shoemaker mission in 2000 (Thomas et al. 2002). The surface is covered with craters of different sizes and stages of degradation.
In the text

Fig. 6. Profile lines of all kernels used to generate the different morphologies of different types of craters. The simple crater 1 (blue line), simple crater 2 (red line), fresh crater (green line), and complex crater (yellow line) correspond to kernel1, 2, 3, 4, resp. The kernel1 is a standard Gaussian function with the true height Ht = 1, the apparent diameter . In kernel2, Ht = 1 and (see Eqs. (4) and (5)). The fresh crater is created using kernel3 with Ht = 1,  Da = 1.31 and extra parameters Drim ≈ 2,  Hrim ≈ 0.309, which represent the rim-to-rim diameter and the rim height, respectively (see Eqs. (7)–(11)). The new degradation length, the rim-to-rim diameter and the rim height of the complex crater are defined as Drim = 2,  Hrim = 0.471,  Hc = 0.2 (see Eqs. (13)–(16)).

In the text

	Fig. 7. Examples of different crater types that can be generated on a small body. The uppermost image shows the simple crater (type 1 and 2) using kernel1, 2. The middle image represents fresh craters generated from the kernel3, and the bottom image depicts the morphology of the complex crater from kernel4.
In the text

Fig. 8. Profile lines showing crater degradation. We can model the degradation status of craters by changing the parameters in our kernel functions. In simple craters 1 & 2, the solid line represents the original crater, while the dashed line represents the same crater but with a lower (degraded) depth-to-diameter ratio (ratio of solid/dashed line = 0.4/0.1). In fresh craters, we set the degradation vector α = (Δk = −0.3, Δγ = −0.4, Δa = −0.25) (defined in Sect. 3.3). For complex craters, we set the degradation vector α = (Δk = −0.05, Δγ = −0.2, Δa = −0.1).

In the text

Fig. 9. Example of crater populations generated for an asteroid of size comparable to Dimorphos. We constrain the distribution by setting exactly one crater at the largest size (here 100 m). We then estimate the number of craters at smaller sizes with a cumulative power-law distribution (dotted line). Discrete points show several crater populations generated from our inversion sampling method. Deviations from the ideal power law are due to the nature of fitting a discrete population to a continuous cumulative function. The colored area represents the 95% confidence interval.

In the text

	Fig. 10. Application of the size-frequency distribution (SFD) defined in Fig. 9 to the asteroid Dimorphos. Both images have been generated using the same crater production function (cumulative power-law distribution with B = 3). Due to randomization, we can generate an infinite number of similar “look-alike” shapes.
In the text

Fig. 11. Crater aging: We applied the crater degradation model of Fassett & Thomson (2014) to three craters of different diameters (D = 300 m, 1 km, 3 km) from an age of 0 to 3 Gyr. Each line in those figures corresponds to a certain age interval ∼1 Gyr (0 Gyr, 1.0 Gyr, 2.0 Gyr, and 3.0 Gyr). After 3 Gyr of evolution, the smallest crater is almost indiscernible, the 1 km crater is flattened but still has distinct relief, and the 3 km crater only shows apparent changes near its rim and interior floor.

In the text

	Fig. 12. Simulated Dimorphos-like asteroid, with low crater density. We add several layers of craters on the surface, e.g., simple crater types 1 and 2, and fresh craters. The total saturation is nearly 20% and the slope of power law is B = 3. The depth-to-diameter, , ratio of simple craters is in the range [0.15, 0.25].
In the text

	Fig. 13. Simulated Dimorphos-like asteroid, with high crater density. As in Fig. 12, we add several layers of craters (simple craters 1 and 2, fresh craters). The total saturation is nearly 30% and the slope of the power-law is B = 4. The depth-to-diameter ratio, , of simple crater is in [0.12, 0.19]. In fresh craters, we set the degradation vector as: α = (Δk = −0.15, Δγ = −0.1, Δa = −0.25).
In the text

	Fig. 14. Simulated Didymos-A-like asteroid, with low crater density. We added several layers of craters on the surface with different types (i.e., simple crater types 1 & 2, complex, and fresh craters). The slope of power-law is B = 3. The depth-to-diameter ratio, , of simple craters and fresh craters is in [0.1, 0.2].
In the text

Fig. 15. Simulated Didymos-A-like asteroid, with high crater density. We add several layers of craters, like simple crater 1 & 2, complex and fresh crater, on the surface. The slope of power law is B = 3. The depth-to-diameter ratio, , of simple craters is in [0.08, 0.16], and for fresh craters is in [0.12, 0.19]. For fresh craters, we set the degradation vector α = (Δk = −0.2, Δγ = −0.2, Δa = −0.38). For complex craters, we set the degradation vector α = (Δk = −0.04, Δγ = −0.4, Δa = −0.2).

In the text

	Fig. 16. Simulated Eros-like asteroid, observed at different phase angles (solar elevation). The shading is Lambertian and does not attempt to describe complex processes such as the opposition effect, but rather aims to display the extent of shadowed areas as a function of the phase.
In the text

	Fig. 17. Effective gravity and gravitational slopes on a simulated Eros-like asteroid.
In the text

	Fig. 18. Computational time of our algorithm to generate different asteroid models with different layers (crater bins). For each asteroid (Ellisoid, Eros, Bennu) we tested 1, …, 5 layers. The bulk of the time is spent on extracting the polygonal mesh.
In the text

	Fig. 19. Simulated views of asteroids Didymos and Dimorphos, as they will be observed by the Hera mission (ESA). Left panel: using the object’s shape according to current best knowledge, obtained from radar observations (Naidu et al. 2020). Right panel: same models, upscaled to high resolution topography with the technique described in this paper.
In the text

	Fig. A.1. Simulated Eros-like asteroid, with low crater density. Just as in Figs. 14 and 14, we added several layers of craters (the surface partitioned into nearly 1500 patches). The total saturation is nearly 11% and the slope of power law B = 3. For both simple and fresh craters, .
In the text

Fig. A.2. Simulated Eros-like asteroid, with high crater density (the surface partitioned into nearly 1500 patches). The total saturation is nearly 18% and the power law index B = 3. For simple craters, , while for fresh craters . For fresh craters, we set the degradation vector α = (Δk = −0.15, Δγ = −0.2, Δa = −0.38). In complex crater, we set the degradation vector α = (Δk = −0.2, Δγ = −0.3, Δa = −0.25).

In the text

	Fig. A.3. Simulated Bennu-like asteroid, with low crater density. The total saturation is nearly 19% (the surface partitioned into nearly 2000 patches) and the slope of power law B = 3. For both simple and fresh craters, .
In the text

Fig. A.4. Simulated Bennu-like asteroid, with high crater density. The total saturation is nearly 25% (the surface partitioned into nearly 2000 patches) and the slope of power law B = 3. For simple craters, we have , while for fresh craters . For fresh craters, we set the degradation vector as α = (Δk = −0.2, Δγ = −0.2, Δa = −0.38). In complex crater, we set the degradation vector as α = (Δk = −0.04, Δγ = −0.4, Δa = −0.2).

In the text