SISPO: Space Imaging Simulator for Proximity Operations

This paper describes the architecture and demonstrates the capabilities of a newly developed, physically-based imaging simulator environment called SISPO, developed for small solar system body fly-by and terrestrial planet surface mission simulations. The image simulator utilises the open-source 3-D visualisation system Blender and its Cycles rendering engine, which supports physically based rendering capabilities and procedural micropolygon displacement texture generation. The simulator concentrates on realistic surface rendering and has supplementary models to produce realistic dust- and gas-environment optical models for comets and active asteroids. The framework also includes tools to simulate the most common image aberrations, such as tangential and sagittal astigmatism, internal and external comatic aberration, and simple geometric distortions. The model framework’s primary objective is to support small-body space mission design by allowing better simulations for characterisation of imaging instrument performance, assisting mission planning, and developing computer-vision algorithms. SISPO allows the simulation of trajectories, light parameters and camera’s intrinsic parameters.


Introduction
A versatile image-simulation environment is required in order to design advanced deep-space missions, to simulate large sets of mission scenarios in parallel, and to develop and validate algorithms for semi-autonomous operations, visual navigation, localisation and image processing. This is especially true in the case of Small Solar System Body (SSSB) mission scenarios, where the mission has to be designed with either very limited information about the target (i.e., precise size, shape, exact composition and activity) or the targets can remain a near-complete mystery before their close encounter (i.e., as in the case of interstellar objects [1,2]). Some publicly known cosmic-synthetic-image generators for space missions are available, and they are briefly described in the next paragraph. Agency (ESA) announcement of opportunity for "New Science Ideas" [27] and it was not chosen despite reaching the final three [28]. Currently, the SISPO simulator is actively used for Optical Periscopic Imager for Comets (OPIC) development [29], which will be hosted on one of three spacecraft making up the Comet Interceptor mission (https://www.cometinterceptor.space). Comet Interceptor is ESA's first F-class mission (in cooperation with Japan Aerospace Exploration Agency (JAXA)) to fly by either a dynamically new comet approaching the Sun for the first time from the Ö pik-Oort cloud, an interstellar object, or a long-period comet as a backup target [30]. OPIC will use automatic image capturing-algorithms will be developed and tested using photorealistic 3-D renderings and a reconstruction pipeline of SISPO using a set of possible encounter velocities and geometries, as well as cometary and camera properties. From a scientific point of view, the simulator will also be used to develop image-prioritisation algorithms, which are required due to the limited data budget, short fly-by timeline and the possibility of the probe being damaged by high-velocity dust impact. The EnVisS coma mapper [31] of the Comet Interceptor mission also uses SISPO for algorithm development.

General architecture
The general structure of SISPO comprises the following parts: 1. Core features: a. Image rendering using Blender Cycles (see Subsection: Rendering in Blender Cycles) or OpenGL (see Subsection: Lightweight OpenGL-based rendering); b. Simulation of on-board image processing. Currently the image processing is primarily related to compression; however, inclusion of both cropping and image prioritisation is planned; c. 3-D reconstruction (see Section: Usability of images produced for 3D reconstruction).

Additional models:
a. Gas and dust environment (see Subsection: Gas and dust); b. Camera distortions (see Subsection: Camera); c. Attitude dynamics in the initial stage.
The core functionalities are split into three subpackages. The first subpackage uses Keplerian orbit data for an SSSB and a simplified definition of the encounter geometry so the spacecraft can propagate realistic trajectories using the Orekit library [32] and render an image series of the encounter. The second subpackage provides various algorithms for image compression and decompression. The third subpackage uses images to reconstruct a textured 3-D model using the Structure from Motion (SfM) technique. The three subpackages combined provide a processing pipeline from an initial 3-D model to a reconstructed 3-D model via rendered and compressed images. SISPO is a Python software package that is hosted on a public GitHub repository under a GPL v3.0 licence and is maintained by the authors [33] among other contributors (e.g., authors of this paper). A description of the core functionality and the additional models is provided in Fig 2.

Rendering system
The most crucial part of SISPO is image synthesis. Two separate rendering modes are implemented: Blender Cycles and OpenGL.

Rendering in Blender Cycles
Blender Cycles uses path tracing, which is a type of ray tracing. In classic ray tracing [34], each camera's pixel shoots one or multiple rays, which interact with a surface and then encounters light sources directly or interacts with other surfaces before reaching the light. For realistic reproduction, various optimisations can be used. The ray-surface and ray-volume interactions can be modelled by shaders, which model or approximate the interacting objects' properties.
The Cycles rendering engine supports various shaders, the most relevant of those for this paper being: (i) diffuse bidirectional scattering distribution function, which provides access to Lambertian and Oren-Nayal shaders based on surface roughness; (ii) emission shader that allows surfaces or volumes to emit light; (iii) subsurface scattering that supports cubic, Gaussian and Christensen-Burley models; (iii) glossy shader that supports Sharp, Beckmann, GGX, Asikhmin-Shirley and multi-scatter GGX models; (iv) volume scattering shader that allows simulating light scattering in volumes [15]. These shaders can be combined with both procedurally generated or pre-existing texture maps to change their parameters and mix various shaders on and within the scene models. Cycles also supports Open Shading Language, which allows the use of arbitrarily defined shaders.
Path tracing is, in some way, an improvement on ray tracing; it produces multiple rays from the same pixel in random directions, then each ray keeps bouncing (without producing new rays) until it reaches the light source or user-defined bounce limit. In Blender Cycles, one can limit the following bounce limits (ideally it should be infinite; however, typically, a limited amount is sufficient): total, diffuse, glossy, transparency, transmission and volume scattering. Then, the amount of light per pixel is calculated for each ray along with surface colour properties; afterwards, the value is averaged and assigned to that specific pixel. Moreover, Blender has a branched path tracer, which can be used for volumetric scattering (see Subsection: Case 2: volumetric particle effects). The branched path tracer is similar to path tracing, but at the first ray interaction it will split the path for different surface components, and for shading, it will take all lighting parameters into account instead of just one [15]. The branched path tracer can also be useful for solving light-related problems such as caustics. All the surface features can be generated procedurally (e.g., by applying mathematical equations instead of externally-captured image texture) inside the Blender software [35]. The SISPO surface model includes sand flats, rock formations and craters. The size of these features, their number and their distribution can be adjusted by modifying the corresponding parameters.
The elevation details are simulated by using height maps for craters, rocks and sand inside the rendering engine (Fig 3). These consist of textures converted into surface displacement distances and are used to displace parts of the model mesh and provide surface normal changes.
Height maps can be mixed together to simulate different types of terrestrial bodies (see Section: Simulations of realistic asteroid imagery). Each map's weight in the mix can be modified, or even obliterated, and they can be made to affect only specific areas of the mesh by using masks. The shader adds a texture whose albedo corresponds to the albedo measured in real asteroids in addition to the procedural elevation. The final procedural shader example is shown in Fig 4.

Lightweight OpenGL-based rendering
There is an option to use a custom OpenGL-based renderer instead of Cycles and Blender. It was added because it is desirable to generate images fast in certain situations even though some fidelity (e.g., softer shadows or procedurally generated details) might be lost. For instance, during the development of image-processing algorithms, it is useful to generate large datasets for training data in a relatively short time. Validation data can be generated with Cycles to give better confidence in the characterised algorithm performance. Fast image generation also enables a closed-loop simulation of a guidance, navigation and control system that depends on the navigation camera input.
Three different BRDFs have been implemented as OpenGL shaders: Lambertian [4], Lunar-Lambertian [36] and Hapke [37]. Textures are supported and correspond to the single scattering albedo of the corresponding shape model region.
Efficient shadowing is achieved by shadow mapping [38]. Shadow mapping works by first rendering the scene orthographically from the direction of the Sun. The resulting depth buffer contents are imported to the actual rendering pass as a texture. The object vertices are again projected orthographically towards the Sun; the corresponding fragment is considered to be in the shadow if the depth of that projection is greater what is in the shadow texture. The output of the OpenGL render is a floating-point value giving the irradiance [W � m −2 ] incident on each theoretical pixel. This "irradiance image" can be further passed to the default camera model used by SISPO, which then applies appropriate imaging effects to get the final synthetic image. There is currently no procedural shape or texture generation due to the simplified nature of the OpenGL renderer. Dust and gas coma rendering is not fully integrated with the OpenGL renderer. Instead, emission-based rendering of a previously generated coma is done separately in Python after OpenGL rendering (see Subsection: Case 2: volumetric

PLOS ONE
particle effects), and it does not take into account shadows cast by objects in the scene. One fundamental limitation to OpenGL-based rendering is that it does not support extended light sources. Thus, for instance, earthshine in low Earth orbit-or in the case of a binary asteroid, the reflected light from the primary body to the secondary body-cannot be reasonably modelled. As the light from reflecting surfaces is not considered, the shadow of rocks, cliffs and other geological formations will be darker than in reality.
The OpenGL rendering engine was used for a preliminary guidance, navigation and control analysis done for the preliminary design review of the now-defunct Asteroid Prospection Explorer project. For this purpose, the rendering engine was interfaced from a Simulink model designed to handle the system dynamics. The same model then forwarded the generated images to the visual navigation algorithms for position estimates. The generation of one image took around 0.7-1.0 seconds. The simulation time for one mission week was around 12 hours on a laptop with an Intel i7-7700HQ (2.8 GHz) processor, with most of the processing time spent on executing the visual navigation algorithms. Fig 5 demonstrates the comparison between the image taken by the AMICA instrument of the Hayabusa spacecraft and rendered image using OpenGL. The main rendering discrepancies are induced by the inaccuracy of the 3D model [14] and local variations of albedo and roughness.

Gas and dust
In the inner Solar System, the solar radiation heats nuclei of comets, which causes them to release gas together with dust. This gas and dust surround a comet forming a coma, which is blown by the solar wind resulting in ion and dust tails that can extend over 100 million kilometres and might be visible from Earth. Sometimes the comet releases strong outbursts of these gasses and dust to produce distinctive coma-derived features called jets [39]. The gas and dust environment is visible to the camera instrumentation, and hence, in order to provide realistic space-scene simulation, the visible effects should be included in modelling.
The dust and gas environment causes noise that can affect the 3-D shape reconstruction. SISPO offers models for comae, jets and tails with various details. The problem can be divided into two parts: modelling the environment and rendering it. The modelling parts can contain, for instance, the description of the environment, such as jets represented as geometric cones from the surface of the comet, mathematical models from [40], or gas and dust simulations such as those from [41,42].
The rendering problem comes down to the level of desired details and accuracy, but on the other hand, to reasonable computing time. In reality, the volume is a mixture of gas and dust, which are composed of different particle sizes with different densities, and this produces various scattering characteristics [43]. Realistic and accurate representation of the effects mentioned above in the rendering pipeline, which is not specialised in volume scattering, is a challenge; a visually realistic approximation should be sought. In the current version of SISPO, volume-scattering effects from the coma and jets are computed using a simple volume-scattering shader from Cycle that uses the voxel presentation of the comet's surroundings as an input for the density.

Camera
The Optical Aberrations for Still Images Simulator (OASIS) provides simulation tools for optical aberrations that are usually not implemented in popular 3-D software (e.g., Blender). Since it uses two-dimensional images as input, motion blur effects, which require spatial awareness of a scene, are not modelled. Optical Aberrations for Still Images Simulator (OASIS) is used in a complementary manner with SISPO to further enhance rendered two-dimensional output images.
Tangential and sagittal astigmatism, as well as an internal and external comatic aberration, are modelled with distinct PSFs, which vary with the field height and orientation of the sensor centre-image point vector. By default, aberration intensity increases with the field height. However, a custom lens file can be provided to model any desired lens array.
Lateral chromatic aberration is modelled by rescaling individual colour channels and simulating wavelength-dependent refraction of light rays. While wavelengths are continuous in the real world, the digital Red, Green and Blue (RGB) triplet only distinguishes between three discrete primitive colours, which introduces sharp edges at the boundaries of colour separation. These can be avoided by blending chromatic aberration with tangential astigmatism (more information with an example can be found in Section 5.7 in [44]).
Dark-current noise is currently drawn from a folded normal distribution with a zero mean and user-defined standard deviation. Readout noise is modelled by the addition and subtraction of random values at the subpixel level, governed by a Gaussian custom standard deviation process. The projection of light rays with random origin positions generates realistic shot noise and follows the user-defined average sample size per pixel.
Lens distortion is simulated with a sophisticated Brown-Conrady model [45], that corrects for both radial and tangential distortion. It is implemented in the computer-vision library OpenCV [46] and supports up to six radial distortion coefficients-k 1 to k 6 -and two tangential distortion coefficients-p 1 and p 2 .
Monochrome sensors are modelled by averaging an RGB colour triplet of a virtual light ray. Additionally, a weighted-average model can be selected that reflects the indeed perceived luminosity. The generation of dark current and readout noise is adjusted accordingly to maintain a specified standard deviation on monochromatic detectors. The simulation of sagittal astigmatism, coma, chromatic aberration, shot noise and readout noise is demonstrated in more detail in [44] (Section 5.7).
Currently, OASIS is limited to simulating one type of PSF at a time, as a realistic convolution of multiple PSFs is not yet provided, with the exception of small aberrations. Also, the generation of dark noise, which should follow a Poisson distribution, is subject to change once physical units are implemented for photon flux and for the conversion between light-ray energy and digital-sensor value.

Orbital simulation
The trajectory simulation within SISPO is handled by the Orekit library [32]. A simple Keplerian orbit propagator is used to propagate both the SSSB and the spacecraft. Additionally, it is possible to rotate the SSSB around a single axis at a constant spin rate.
The implemented Orekit Python bindings run a virtual machine to execute the underlying Java code. During propagation, Orekit determines state information of the SSSB and the spacecraft for each sample along the trajectory. The state information includes the date, position and the rotation angles of the SSSB.
In SISPO, the spacecraft trajectory is normally defined by its Keplerian elements, but the user does not explicitly enter these elements. Instead, the elements are calculated from the target body's Keplerian elements and the expected encounter geometry relative to the SSSB at the closest approach. The parameters presented in Table 1 are used to calculate the state vector of the spacecraft at the encounter, which defines the spacecraft trajectory. The sssb_state is calculated based on the SSSB input data and the encounter data.
The propagation process is defined by the duration of a fly-by, the number of frames to be rendered, timesampler mode and a slowmotion factor as presented in Table 2. The timesampler mode determines whether the steps are distributed linearly in time (mode 1, default) or whether an exponential model (mode 2) is used, which increases the number of frames around the encounter. The number of additional samples can be controlled with the slowmotion factor.

PLOS ONE
Mode 2 is especially helpful when simulating a long fly-by, since, when the spacecraft is far from the SSSB nucleus, only minor changes are visible in rendered images.

Attitude dynamics
Along with the orbital simulation, a simplified attitude dynamics portion is already built into the Orekit framework, also accessible via the SISPO library. Attitude information in Orekit is handled essentially as another frame transformation. It contains the rotation from the reference frame to the satellite frame, and then the angular velocity (i.e., spin) and angular acceleration (i.e., rotation acceleration) of the spacecraft in its frame. This makes Orekit's spacecraft_state a great construct to hold attitude-related and orbital data, which can then successfully be propagated in time. The Orekit library also allows SISPO to build user-defined "attitude law". This attitude law can be selected from several pre-existing and common attitude modes (e.g., pointing, spin stabilised) or recompiled entirely in a way necessary for the mission simulation. It is, however, essential to note that the Orekit library mainly focuses on orbital mechanics and propagation. The attitude model, and any attitude laws it follows, presumes the existence of a "perfect attitude control system" without necessarily considering the physical limitations or perturbations of the satellite attitude by external forces. The Orekit library lacks the necessary constructs to convey and utilise the moments of inertia of the satellite. Since this information is crucial for the inclusion of external forces and their effect on the spacecraft's attitude or other celestial bodies, this will be implemented as a secondary layer into the current framework; this is explained in Section: Discussion and future work.

Simulations of realistic asteroid imagery
This section demonstrates five different use cases of the SISPO software.  Fig 7 shows a comparison between a plain mesh of the asteroid 25143 Itokawa and a mesh recreation using the procedural shader in Blender. The initial mesh is a 3-D reconstruction of the actual asteroid as described in [47]. The shader then adds procedural surface features on top.

Case 1: Asteroid 25143 Itokawa
A small area of the model mesh was rendered using OpenGL and Cycles to show the ultimate advantage of procedural texturing. Both renders are shown in Fig 8, and the one with procedural texturing indicates the capability to generate almost arbitrary level of detail. This feature would be required, for example, to train probes that can land on the surface.

Case 2: Volumetric particle effects
Coma and dust jets are produced by particles emanating from the surface, and these are used to produce a volumetric density distribution. For the preliminary proof of concept, the jets were modelled using simplified versions of the gas and dust models from [41,42]. First, the gas source strengths and gas velocities are generated for each face of the mesh from the noise function. The gas source data are then used as input for the gas model [41], which computes the gas field densities and velocities around the comet. This data is then used in a particle simulation following [42] from which the particle position data is further processed into a threedimensional number-density map around the comet, which is encoded to a single OpenEXR file [48]. Before the render, the number-density map is smoothed with tricubic interpolation [49]. It can then be loaded by Cycles and rendered using either volumetric scattering or emission (volumetric emission, although less accurate, is orders of magnitude faster).
Volumetric particle effects on the comet 67P/C-G. In Fig 9 (the coma might appear differently to the article's viewer depending on the monitor or the print quality), the dust environment capabilities of SISPO are presented by simulating and rendering images of the comet 67P/C-G. This image is then compared to the actual image, and with an OpenGL rendered  image. All the shader features (rocks, craters and sand flats) can be used simultaneously for complex SSSBs if needed. The masks that separate the "sand flats" from the rest of the surface features are created with procedural noise, but they could be painted by hand if necessary.

Case 3: Larger bodies
The reproduction of larger terrestrial bodies in SISPO is demonstrated in Fig 10, which is based on the narrow-angle-camera digital terrain model of the Apollo 15 landing site and operations area [50], obtained by Lunar Reconnaissance Orbiter Narrow Angle Camera (http://wms.lroc.asu.edu/lroc/view_rdr_product/NAC_DTM_APOLLO15_M111571816_ 50CM, accessed 28.01.2021). The original terrain model has a 2 m resolution for elevation maps and 0.5 m resolution for orthographic photos, which is not enough to produce images from the surface that would be useful for navigation testing. This case demonstrates the image scalability produced in SISPO. The lunar rendering was made by implementing a digital terrain model and procedural texturing over the whole 5.1 × 28.6 km 2 area; the farthest point visible on the render is roughly 11 km away from the camera.

Case 4: Subsurface exploration
Lava caves, the isolated underground environments, exist on Moon [51] and Mars [52]. These caves have been studied by remote sensing and have not been explored by dedicated missions, although some were proposed and developed, such as Moon Diver [53] and rock climber Lemur [54]. Lava tubes have been morphologically related to ones formed on Earth in volcanic rock by a volcanic eruption [55]. SISPO could be utilised for synthesising versatile sets of lavacaves images for (i) developing navigation algorithms and sampling spots detection by image processing (e.g., biological mats or their preserved remains on the geological substrate) or (ii) mapping the caves by photogrammetry for potential human settlements. An example of SISPO cave rendering is demonstrated in Fig 11.

Case 5: Fly-by of a spacecraft
The SISPO environment is suitable for spacecraft fly-by rendering. This can be useful for multi-spacecraft missions, in-orbit servicing and fleets. The Cycles rendering of the spacecraft fly-by is demonstrated in Fig 12. The model uses a principled bidirectional scattering distribution function shader including multiple layers to create spacecraft materials. The set of produced images can be used for developing and training formation-flying proximity algorithms demanded by attitude control and orbit determination subsystem.

Usability of images produced for 3D reconstruction
A set of synthetically generated images can be used for photogrammetry-based 3-D surface reconstruction within the SISPO environment. The steps executed in the reconstruction pipeline are described in [16]. The reconstruction uses two libraries: 1. Open Multiple View Geometry (OpenMVG) C++ library [56], which creates a sparse point cloud based on two different SfM techniques:

Open Multiple View Stereo Reconstruction (OpenMVS), which uses input from
OpenMVG, creates a dense point cloud, a faceted surface (mesh) or a set of planes [59].

Case 1: SSSB fly-by
A set of 25 images was generated in SISPO and then used to reconstruct the 3-D model using the pipeline. The images were generated using a pinhole camera (i.e., no optical aberrations,  the fly-by (i.e., during the fly-by only one illuminated part of the target is visible). The reconstructed model was compared with the input 3-D model externally using the CloudCompare software [60]. Generally, the reconstructed model is close to the input model, except for the edges (e.g., the black area in Fig 15), where the error deviation was very high (up to 600 m). These highly deviated parts could be removed in the point cloud or mesh post-processing, but have not been done in this analysis. The Gaussian distribution mean error is 14.7 m and the standard deviation is 63.9 m. The visualisation of deviation and the error analysis are shown in Fig 15.

Case 2: SSSB orbiting
During SSSB orbiting, it is possible to observe every side of the target because as the target spins, each part will be illuminated at a certain point. Three different input-image views (out of 53 used in total) and a fully reconstructed 3-D model are shown in Fig 16. Fifty-three images were generated from different angles using a pinhole camera.
The input and reconstructed 3-D model were compared using the same method as in the fly-by case. Because of the larger set of observational angles, the deviation is decreased in the orbiting case and reaches 60-89.7 m at some crater spots, but most of the surface is aligned with the input model as shown in Fig 17. The Gaussian mean error distribution is 1.2 m and the standard deviation is 6.6 m. The accuracy can be increased by using a more significant set of images and higher resolution.

Discussion and future work
SISPO is a sophisticated tool for terrestrial-cosmic-scenery rendering. Its most significant advantage is the physically based, high-quality, realistic renderings, which can reproduce submillimeter surface resolution and beyond for SSSBs and other terrestrial bodies thanks to micropolygon procedural texturing and Blender's Cycles path-tracing rendering engine. The produced renderings are mature for deep-space mission design and algorithm development for semi-autonomous operations, visual navigation, localisation and image processing. There are a few more functions that are considered for future implementation: • Attitude dynamics auxiliary package (see below);

PLOS ONE
• Image compression algorithms for efficient data storage and transmission evaluation, which was preliminarily assessed by [61]; • Reconstruction of various scenery; • Algorithms for spacecraft photogrammetry-based localisation; • Capability to simulate measurements of the target spectral reflectance in certain wavelength channels; • Improved shader and user-defined parameterisation of Cycles in order to minimise the approximation; • Solar System ephemeris integration for historical and upcoming events. This can include both terrestrial bodies and spacecraft via the possible adaptation of Spacecraft Planet Instrument Camera-matrix Events (SPICE) data.

Attitude dynamics auxiliary package
The SISPO simulation environment will include optional packages for calculating the effects of external forces such as dust particle impacts, atmospheric drag, gravity gradient and various other effects necessary to account for. Each of these functionalities will be built as an add-on to the current Orekit-based attitude framework. This is done by adding an optional extra calculation step between the propagations of the spacecraft state through time by Orekit to account for the perturbations. With the addition of the second layer on top of Orekit orbital simulations, for accurately simulating the attitude dynamics of the spacecraft or a celestial body, the SISPO environment will have the necessary groundwork for later expansions into different active and passive attitude control systems present in any future mission simulations.

Conclusions
In this paper, the Space Imaging Simulator for Proximity Operations (SISPO) architecture and capabilities were described and several features demonstrated by case simulations. The tool generates physically based, photorealistic images from an input 3-D model using Blender's Cycles path-tracing rendering engine. For regolith surface reconstruction it uses procedural texturing. SISPO has supplementary models for optical aberrations as well as for gas and dust of small bodies. The description and usage of these models were discussed in this paper. Other use cases have demonstrated renderings of the Moon and a spacecraft. The set of produced images can be implemented for 3-D surface reconstruction. The tool is open access and currently requires basic programming skills to be used [33]. The level of detail currently produced by SISPO is suitable for the design of advanced deep-space missions, the simulation of large sets of scenarios, and the development and validation of algorithms for (semi-)autonomous operations, vision-based navigation, localisation and image processing. SISPO has already supported the development of deep-space mission concepts and is currently being used for the ESA-JAXA Comet Interceptor mission (more details in Section: Application to space mission designs). Some further improvements have been considered and are listed in Section: Discussion and future work. Other teams are welcome to use SISPO, implement new functions and contact the authors of this paper.
Supporting information S1 Script. Image comparison. The file contains python algorithm used to compare images in this manuscript. (ZIP)