Issue 
A&A
Volume 659, March 2022



Article Number  A176  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202140709  
Published online  25 March 2022 
Numerical approach to synthesizing realistic asteroid surfaces from morphological parameters
^{1}
University of Bremen, Bremen, Germany
email: lixizhi@unibremen.de
^{2}
DLR Institute for Planetary Research, Berlin, Germany
^{3}
Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai, PR China
Received:
26
April
2021
Accepted:
29
November
2021
Context. The complex shape of asteroids and comets is a critical parameter in many scientific and operational studies. From the global irregular shape down to the local surface details, these topographies reflect the formation and evolutionary processes that remould the celestial body. Furthermore, these processes control how the surface will continue to evolve: from mass wasting on high slopes to spinup due to anisotropic reemission of thermal radiation. In addition, for space missions, the irregular coarse shape and complex landscape are a hazard to navigation, which must be accounted for in the planning phase.
Aims. In this paper, we propose a novel method to synthesize physically correct 3D shape models of small celestial bodies, such as asteroids, to support the testing of a wide range of parameters in scientific and operational studies.
Methods. We modeled virtual asteroid shapes using nonuniform sphere packings to represent the coarse shape, define an implicit surface, and then synthesize highresolution topography with userdefined, locally controlled spot noise models. This effectively replaces the random noise model (e.g., Perlin noise) used in traditional approaches and allows us to construct a morphology based on actual physical shapes of the most common features observed on asteroids and comets. As an example of such a feature, we propose several kernel functions to add virtual craters to the coarse shape of the asteroid, of which the spatial distribution is controlled by typical crater production functions (e.g., a power law).
Results. We demonstrate how this technique can be used to generate a variety of asteroid shapes and topographies using different cratering parameters and distributions. We apply our technique to artificially increase the resolution of existing models of the DidymosDimorphos system, the target of the Double Asteroid Redirection Test, and Hera missions. We show that our approach generates models that are suitable for typical analysis relying on detailed asteroid shapes, as well as operational scenarios for space missions. The meshes created with our algorithm can be directly used with existing visualization software and operations or science pipelines to generate data suitable for mission planning and to validate data analysis techniques.
Key words: methods: numerical / minor planets / asteroids: general
© 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, microimpacts 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 lowresolution 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 spinup 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 YarkovskyO’KeefeRadzievskiiPaddack effect (YORP; Rubincam 2000) changes the spin rate of an object due to the anisotropic reemission 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 “Vshape” dispersion pattern in collisional families (Bolin et al. 2017).
The diskintegrated 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 phasedarkening 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 (ElkinsTanton & 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 fuelefficient 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 shapefromshading 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 powerlaw 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^{1} 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 parameterfree texturing method to add features through displacement mapping.
Hypertextures (Ebert et al. 2003) as a typical parameterfree 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 selfsimilar 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 polygonlike 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 microstructural 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 spectraloriented functions into spatialoriented 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 microdetails 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 selfsimilar 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 largescale landforms with mesoscale or smallscale 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 flyby 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 multiellipsoids 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 groundbased 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 silicaterich chondrules, metalrich matrix) and produces a variety of morphologies (from bowlshaped 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 threedimensional space ℝ^{3} as:
where r represents the radius, 0 is the isovalue 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 twostep 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 nonuniform 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 watertight constraint shape. Hence, it is possible to adapt this implicit representation to any existing polygonal shape like the lowpoly 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 lookalike 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 (nonuniform) 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 highresolution, highdetail mesh. 
Fig. 3.
Results of our shape synthetization system. From the constraint shape, namely, the real lowpoly 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 lookalike 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., bowlshaped craters, polygonal craters, floorfractured craters, secondary craters, and crater chains) and diverse sizefrequency 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^{3} cells of the next more finegrained 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 sizefrequency 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:
with ω ∈ [0, 1] a usercontrollable 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 bowlshaped 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 bowlshaped craters, polygonal craters, floorfractured 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 concaveupward 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 simpletocomplex transition diameter, most small asteroids display simple craters, namely, bowlshaped depressions produced by an impact. Simple craters can exhibit a benchlike 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 bowlshape normal craters:
where kernel_{1} is the standard Gaussian function (q_{1} = 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 bowlshape 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 craterfill 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} = 2 and obtain kernel_{2} 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).
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 kernel_{1, 2, 3, 4}, resp. The kernel_{1} is a standard Gaussian function with the true height H_{t} = 1, the apparent diameter . In kernel_{2}, H_{t} = 1 and (see Eqs. (4) and (5)). The fresh crater is created using kernel_{3} with H_{t} = 1, D_{a} = 1.31 and extra parameters D_{rim} ≈ 2, H_{rim} ≈ 0.309, which represent the rimtorim diameter and the rim height, respectively (see Eqs. (7)–(11)). The new degradation length, the rimtorim diameter and the rim height of the complex crater are defined as D_{rim} = 2, H_{rim} = 0.471, H_{c} = 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_{2} contains the central floor) and walls (Marchi et al. 2015). Consequently, another type of crater we are aiming at is typically bowlshaped 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 bowlshaped 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_{3} to generate younger fresh craters with rim and flanks:
where κ and b control the height above the 0plane; β and γ define the radii of the rim above the 0plane; and a controls the depth of the crater under the 0plane.
Similarly, we can define five variables H_{rim}, H_{t}, D_{rim}, D_{a}, L_{d} to represent the accurate morphology of a fresh bowlshape crater (see Fig. 6, green line). Since the parameters of kernel_{3} are controlled by our method, we let kernel_{3} = 0 and in that way, we obtain the crater morphology (W is the Lambert W function; Corless et al. 1996):
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 “simpletocomplex transition”, varies according to different planetary bodies. In order to synthesize such large complex craters with central uplift, we propose kernel_{4} to be:
where κ and b control the height above the 0plane, β and γ define the radii of bulge above the 0plane, a controls the depth of the crater under the 0plane. From kernel_{4}, 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_{4} = 0.01 as the border of the complex crater to evaluate the parameter L_{d}:
The definition of our four kernels enables us to accurately control the morphology of different crater types. The easytouse 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 kernel_{1, 2}. The middle image represents fresh craters generated from the kernel_{3}, and the bottom image depicts the morphology of the complex crater from kernel_{4}. 
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 depthtodiameter 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^{n}, 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 depthtodiameter 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 depthtodiameter 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) depthtodiameter 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_{3}), 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_{4}), 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_{3}, kernel_{4}). 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 timeintegrated 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 onetoone correlation between fresheroded and youngold 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 sizefrequency 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 sizefrequency 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 sizefrequency 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 volatilerich 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) energydensity 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 powerlaw 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 subranges; within each range, we generate a number of craters, as specified by the production function (powerlaw). 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 powerlaw 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 powerlaw 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 powerlaw distribution from the inversion sampling (Devroye 1986), we need to compute the probability density pdf(D) from the cumulative powerlaw distribution (Clauset et al. 2009), where:
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^{−1}(u), where u∼ uniform distribution:
so
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 powerlaw approximation hardly matches the SFD curve seen on the surface of a real small body. To accomodate this, we can simply assign different powerlaw indices to different crater size ranges to approximate a more complex SFD curve (consist into a more complex piecewise SFD curve).
Burchell & LeliwaKopystynski (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 & LeliwaKopystynski 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 sizefrequency distribution (SFD) defined in Fig. 9 to the asteroid Dimorphos. Both images have been generated using the same crater production function (cumulative powerlaw distribution with B = 3). Due to randomization, we can generate an infinite number of similar “lookalike” shapes. 
Our fitting of the powerlaw 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 depthtodiameter ratio follows a nearlinear 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 depthtodiameter ratios (in [0.2, 0.25]), then suffer erosion and evolve into much shallower craters (k ≈ 0.05). Small craters have typically a lower depthtodiameter 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 depthtodiameter 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 nonhomogeneous 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^{−7} m ⋅ yr^{−1} (Vincent et al. 2014). Schmedemann et al. (2020) estimate an erosion rate of 6 × 10^{−8} to 3 × 10^{−6} m ⋅ yr^{−1} for kilometersized 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^{−7} m ⋅ yr^{−1}. 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_{3}. 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 depthtodiameter () 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, lowdensity crater surface) and “old” (see Fig. 13, i.e., highdensity crater surface) asteroid surfaces, respectively.
Fig. 12.
Simulated Dimorphoslike 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 depthtodiameter, , ratio of simple craters is in the range [0.15, 0.25]. 
Fig. 13.
Simulated Dimorphoslike 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 powerlaw is B = 4. The depthtodiameter 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 latticebased 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 photorealistic asteroid surfaces using raytracing.
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 powerlaw). When accumulating the layers, the weights should be chosen carefully to avoid unwanted morphologies such as overhangs.
Fig. 14.
Simulated DidymosAlike 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 powerlaw is B = 3. The depthtodiameter ratio, , of simple craters and fresh craters is in [0.1, 0.2]. 
Fig. 15.
Simulated DidymosAlike 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 depthtodiameter 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 Eroslike virtual asteroid for different phase angles. The shadow fraction in each image is an important parameter when trying to understand the effect of largescale 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 Eroslike 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 Eroslike 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 wellknown 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 rescaled 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 highresolution 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:

Diskintegrated 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 raytracing 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 raytracing 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 Eroslike 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 masswasting 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: Highresolution topography and multiscale roughness are necessary for precisely modeling the distribution and evolution of temperatures on an asteroid. Figure 19 shows that preencounter shape models provide a very lowresolution 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 indepth 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 spacebased 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 reaccrets (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.
The executable version is available at: https://cgvr.informatik.unibremen.de/research/procedural_asteroid/
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 NEOMAPP). X.L. acknowledges funding from the China National Space Administration’s Preresearch 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.comettoolbox.com).
References
 Becher, M., Krone, M., Reina, G., & Ertl, T. 2017, IEEE Trans. Visualization Comput. Graphics, 25, 1283 [Google Scholar]
 Belton, M., Chapman, C., Veverka, J., et al. 1994, Science, 265, 1543 [NASA ADS] [CrossRef] [Google Scholar]
 Blinn, J. F. 1982, ACM Trans. Graphics (TOG), 1, 235 [Google Scholar]
 Bolin, B. T., Delbo, M., Morbidelli, A., & Walsh, K. J. 2017, Icarus, 282, 290 [Google Scholar]
 Bottke, W. F., Vokrouhlický, D., Rubincam, D. P., & Nesvorný, D. 2006, Ann. Rev. Earth Planet. Sci., 34, 157 [Google Scholar]
 Bovik, A. C., Clark, M., & Geisler, W. S. 1990, IEEE Trans. Pattern Anal. Mach. Intell., 12, 55 [Google Scholar]
 Burchell, M., & LeliwaKopystynski, J. 2010, Icarus, 210, 707 [Google Scholar]
 Cheng, A. F., Barnouin, O. S., Ernst, C. M., & Kahn, E. G. 2012, in Asteroids, Comets, Meteors 2012, 1667, 6447 [Google Scholar]
 Cheng, A. F., Rivkin, A. S., Michel, P., et al. 2018, Planet. Space Sci., 157, 104 [CrossRef] [Google Scholar]
 Clauset, A., Shalizi, C. R., & Newman, M. E. 2009, SIAM Rev., 51, 661 [NASA ADS] [CrossRef] [Google Scholar]
 Corless, R. M., Gonnet, G. H., Hare, D. E., Jeffrey, D. J., & Knuth, D. E. 1996, Adv. Comput. Math., 5, 329 [Google Scholar]
 de Leeuw, W., & van Liere, R. 1997, in Proceedings of the 1997 ACM/IEEE conference on Supercomputing (CDROM)  Supercomputing ’97 (San Jose, CA: ACM Press), 1 [Google Scholar]
 Devroye, L. 1986, in Proceedings of the 18th Conference on Winter Simulation, 260 [CrossRef] [Google Scholar]
 Dombard, A. J., Barnouin, O. S., Prockter, L. M., & Thomas, P. C. 2010, Icarus, 210, 713 [Google Scholar]
 Ebert, D. S., Musgrave, F. K., Peachey, D., Perlin, K., & Worley, S. 2003, Texturing & Modeling: A Procedural Approach (Morgan Kaufmann) [Google Scholar]
 ElkinsTanton, L. T., & Bell,, J. F. I. 2017, in European Planetary Science Congress, EPSC2017384 [Google Scholar]
 Enos, H. L., & Lauretta, D. S. 2019, Nat. Astron., 3, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Fassett, C. I., & Thomson, B. J. 2014, J. Geophys. Res.: Planets, 119, 2255 [NASA ADS] [CrossRef] [Google Scholar]
 Fassett, C. I., Crowley, M. C., Leight, C., et al. 2017, Geophys. Res. Lett., 44, 5326 [Google Scholar]
 Fujimoto, M., & Tasker, E. J. 2019, Nat. Astron., 3, 284 [NASA ADS] [CrossRef] [Google Scholar]
 Galerne, B., Lagae, A., Lefebvre, S., & Drettakis, G. 2012, ACM Trans. Graphics (TOG), 31, 73 [Google Scholar]
 Gilet, G., Sauvage, B., Vanhoey, K., Dischler, J.M., & Ghazanfarpour, D. 2014, ACM Trans. Graphics (TOG), 33, 195 [Google Scholar]
 Hergenrother, C. W., Maleszewski, C. K., Nolan, M. C., Li, J. Y., & Bernacki, M. 2019, Nat. Commun., 10, 1291 [NASA ADS] [CrossRef] [Google Scholar]
 Hiesinger, H., Marchi, S., Schmedemann, N., et al. 2016, Science, 353, 4756 [CrossRef] [Google Scholar]
 Jorda, L., Lamy, P., & Besse, S., et al. 2010, European Planetary Science Congress 2010, 200 [Google Scholar]
 Kazhdan, M., & Hoppe, H. 2013, ACM Trans. Graphics (ToG), 32, 29 [Google Scholar]
 König, B., Neukum, G., & Fechtig, H. 1977, in Lunar and Planetary Science Conference, 8 [Google Scholar]
 Lamy, P. L., Toth, I., Weaver, H. A., et al. 2006, A&A, 458, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lamy, P., Kaasalainen, M., Lowry, S., et al. 2008, A&A, 487, 1179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levison, H. F., & Lucy Science Team 2016, in Lunar and Planetary Science Conference, Lunar and Planetary Science Conference, 2061 [Google Scholar]
 Lewis, J. P. 1987, ACM Trans. Graphics, 6, 167 [Google Scholar]
 Li, X. Z., Weller, R., & Zachmann, G. 2018, in 2018 15th International Conference on Control, Automation, Robotics and Vision (ICARCV), IEEE, 1771 [Google Scholar]
 Li, X. Z., Weller, R., & Zachmann, G. 2020, in Eurographics 2020  Short Papers, eds. A. Wilkie, & F. Banterle (The Eurographics Association) [Google Scholar]
 Lorensen, W. E., & Cline, H. E. 1987, ACM Siggraph Comput. Graphics, 21, 163 [Google Scholar]
 Mainzer, A., Usui, F., & Trilling, D. E. 2015, Asteroids IV, eds. P. Michel, F. E. DeMeo, & W. F. Bottke (Tucson: University of Arizona Press), 89 [Google Scholar]
 Marchi, S., McSween, H. Y., O’Brien, D. P., et al. 2012, Science, 336, 690 [NASA ADS] [CrossRef] [Google Scholar]
 Marchi, S., Chapman, C. R., Barnouin, O. S., Richardson, J. E., & Vincent, J. B. 2015, Asteroids IV, eds. P. Michel, F. E. DeMeo, & W. F. Bottke (Tucson: University of Arizona Press), 725 [Google Scholar]
 Martin, I., Parkes, S., Dunstan, M., & Rowell, N. 2014, IEEE Comput. Graphics Appl., 34, 52 [Google Scholar]
 McMahon, S. K. 1996, Planet. Space Sci., 44, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Michel, P., Küppers, M., & Carnelli, I. 2018, Cosp, 42, B1 [Google Scholar]
 Naidu, S., Benner, L., Brozovic, M., et al. 2020, Icarus, 348, 113777 [Google Scholar]
 Neukum, G., Ivanov, B. A., & Hartmann, W. K. 2001, Chronology and Evolution of Mars (Springer), 55 [Google Scholar]
 Paris, A., Galin, E., Peytavie, A., Guérin, E., & Gain, J. 2019, ACM Trans. Graphics (TOG), 38, 1 [Google Scholar]
 Paris, A., Peytavie, A., Guérin, E., Dischler, J.M., & Galin, E. 2020, Visual Comput., 36, 2251 [CrossRef] [Google Scholar]
 Pavie, N., Gilet, G., Dischler, J.M., & Ghazanfarpour, D. 2016, Comput. Sci. Res. Notes, 2601, 71 [Google Scholar]
 Pelletier, J. D. 2008, Quantitative Modeling of Earth Surface Processes (Cambridge University Press) [CrossRef] [Google Scholar]
 Raducan, S. D., Davison, T. M., & Collins, G. S., 2019, LPI Contribution No. 2132, 1799 [NASA ADS] [Google Scholar]
 Rubincam, D. P. 2000, Icarus, 148, 2 [Google Scholar]
 Sarli, B. V., Ozimek, M. T., Atchison, J. A., Englander, J. A., & Barbee, B. W. 2017, AAS Paper, 17 [Google Scholar]
 Schmedemann, N., Hiesinger, H., & Michael, G. 2020, in Lunar and Planetary Science Conference, 1450 [Google Scholar]
 Snodgrass, C., & Jones, G. H. 2019, Nat. Commun., 10, 5418 [NASA ADS] [CrossRef] [Google Scholar]
 Srinivas, A., Weller, R., & Zachmann, G. 2017, in ICATEGVE 2017  International Conference on Artificial Reality and Telexistence and Eurographics Symposium on Virtual Environments (The Eurographics Association) [Google Scholar]
 Strom, R., Malhotra, R., Ito, T., Yoshida, F., & Kring, D. 2005, Science, 309, 1847 [NASA ADS] [CrossRef] [Google Scholar]
 Tanga, P., Comito, C., Paolicchi, P., et al. 2009, ApJ, 706, L197 [Google Scholar]
 Thomas, P., Joseph, J., Carcich, B., et al. 2002, Icarus, 155, 18 [Google Scholar]
 Tricard, T., Efremov, S., Zanni, C., et al. 2019, ACM Trans. Graphics (TOG), 38, 1 [Google Scholar]
 van Wijk, J. J. 1991, ACM SIGGRAPH Comput. Graphics, 25, 309 [Google Scholar]
 Vincent, J. B. 2018, in Lunar and Planetary Science Conference, 1281 [Google Scholar]
 Vincent, J.B. 2019, A&A, 624, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vincent, J.B., Besse, S., Marchi, S., et al. 2012, Planet. Space Sci., 66, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Vincent, J.B., Schenk, P., Nathues, A., et al. 2014, Planet. Space Sci., 103, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Vincent, J.B., Hviid, S. F., Mottola, S., et al. 2017, MNRAS, 469, S329 [NASA ADS] [CrossRef] [Google Scholar]
 Weller, R., & Zachmann, G. 2010, SIGGRAPH ASIA [Google Scholar]
 Worley, S. 1996, in Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques, 291 [Google Scholar]
 Wünnemann, K., Collins, G., & Melosh, H. 2006, Icarus, 180, 514 [Google Scholar]
 Zhang, X., Huang, J., Wang, T., & Huo, Z. 2019, LPI Contribution No. 2132, 1045 [Google Scholar]
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” Eroslike asteroid surface, while Figs. A.3, A.4 show the ”young” or ”old” Bennulike 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 Eroslike 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 Eroslike 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 Bennulike 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 Bennulike 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). 
Parameters used in kernel_{3} in Fig. 11.
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 (nonuniform) 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 highresolution, highdetail mesh. 

In the text 
Fig. 3.
Results of our shape synthetization system. From the constraint shape, namely, the real lowpoly 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 lookalike 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 kernel_{1, 2, 3, 4}, resp. The kernel_{1} is a standard Gaussian function with the true height H_{t} = 1, the apparent diameter . In kernel_{2}, H_{t} = 1 and (see Eqs. (4) and (5)). The fresh crater is created using kernel_{3} with H_{t} = 1, D_{a} = 1.31 and extra parameters D_{rim} ≈ 2, H_{rim} ≈ 0.309, which represent the rimtorim diameter and the rim height, respectively (see Eqs. (7)–(11)). The new degradation length, the rimtorim diameter and the rim height of the complex crater are defined as D_{rim} = 2, H_{rim} = 0.471, H_{c} = 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 kernel_{1, 2}. The middle image represents fresh craters generated from the kernel_{3}, and the bottom image depicts the morphology of the complex crater from kernel_{4}. 

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) depthtodiameter 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 powerlaw 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 sizefrequency distribution (SFD) defined in Fig. 9 to the asteroid Dimorphos. Both images have been generated using the same crater production function (cumulative powerlaw distribution with B = 3). Due to randomization, we can generate an infinite number of similar “lookalike” 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 Dimorphoslike 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 depthtodiameter, , ratio of simple craters is in the range [0.15, 0.25]. 

In the text 
Fig. 13.
Simulated Dimorphoslike 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 powerlaw is B = 4. The depthtodiameter 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 DidymosAlike 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 powerlaw is B = 3. The depthtodiameter ratio, , of simple craters and fresh craters is in [0.1, 0.2]. 

In the text 
Fig. 15.
Simulated DidymosAlike 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 depthtodiameter 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 Eroslike 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 Eroslike 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 Eroslike 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 Eroslike 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 Bennulike 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 Bennulike 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.