# ![splashsurf logo](https://raw.githubusercontent.com/InteractiveComputerGraphics/splashsurf/main/logos/logo_small.svg "splashsurf") [![On crates.io](https://img.shields.io/crates/v/splashsurf)](https://crates.io/crates/splashsurf) [![On docs.rs](https://docs.rs/splashsurf_lib/badge.svg)](https://docs.rs/splashsurf_lib) [![Commits since last release](https://img.shields.io/github/commits-since/InteractiveComputerGraphics/splashsurf/latest)](https://github.com/InteractiveComputerGraphics/splashsurf) [![License: MIT](https://img.shields.io/crates/l/splashsurf)](https://github.com/InteractiveComputerGraphics/splashsurf/blob/main/LICENSE) [![Dependency status](https://deps.rs/repo/github/InteractiveComputerGraphics/splashsurf/status.svg)](https://deps.rs/repo/github/InteractiveComputerGraphics/splashsurf) [![Build and test GitHub Actions workflow](https://github.com/InteractiveComputerGraphics/splashsurf/workflows/Build%20and%20test/badge.svg)](https://github.com/InteractiveComputerGraphics/splashsurf/actions/workflows/build.yml) CLI for surface reconstruction of particle data from SPH simulations, written in Rust. For a the library used by the CLI see the [`splashsurf_lib`](https://crates.io/crates/splashsurf_lib) crate.
`splashsurf` is a tool to reconstruct surfaces meshes from SPH particle data. The first image shows the visualization of a set of particles from an SPH fluid simulation from [SPlisHSPlasH](https://github.com/InteractiveComputerGraphics/SPlisHSPlasH). The particle radius is `0.025`. As the rendering of a fluid should not look like a ball pit, a surface mesh has to be reconstructed from this particle data. The next image shows a reconstructed surface mesh of the fluid produced by `splashsurf` with a "smoothing length" of `2.2` times the particles radius and a cell size of `1.1` times the particle radius. The third image shows a finer reconstruction with a cell size of `0.45` times the particle radius. These surface meshes can then be fed into 3D rendering software such as [Blender](https://www.blender.org/) to generate beautiful water animations. The result might look something like this:
Note: This animation does not show the recently added smoothing features of the tool, for more recent rendering see [this video](https://youtu.be/2bYvaUXlBQs). --- **Contents** - [The `splashsurf` CLI](#the-splashsurf-cli) - [Introduction](#introduction) - [Domain decomposition](#domain-decomposition) - [Octree-based decomposition](#octree-based-decomposition) - [Subdomain grid-based decomposition](#subdomain-grid-based-decomposition) - [Notes](#notes) - [Installation](#installation) - [Usage](#usage) - [Recommended settings](#recommended-settings) - [Weighted surface smoothing](#weighted-surface-smoothing) - [Benchmark example](#benchmark-example) - [Sequences of files](#sequences-of-files) - [Input file formats](#input-file-formats) - [VTK](#vtk) - [VTU](#vtu) - [BGEO](#bgeo) - [PLY](#ply) - [XYZ](#xyz) - [JSON](#json) - [Output file formats](#output-file-formats) - [All command line options](#all-command-line-options) - [The `reconstruct` command](#the-reconstruct-command) - [The `convert` subcommand](#the-convert-subcommand) - [License](#license) # The `splashsurf` CLI The following sections mainly focus on the CLI of `splashsurf`. For more information on the library, see the [corresponding readme](https://github.com/InteractiveComputerGraphics/splashsurf/blob/main/splashsurf_lib) in the `splashsurf_lib` subfolder or the [`splashsurf_lib` crate](https://crates.io/crates/splashsurf_lib) on crates.io. ## Introduction This is CLI to run a fast marching cubes based surface reconstruction for SPH fluid simulations (e.g. performed with [SPlisHSPlasH](https://github.com/InteractiveComputerGraphics/SPlisHSPlasH)). The output of this tool is the reconstructed triangle surface mesh of the fluid. At the moment it supports computing normals on the surface using SPH gradients and interpolating scalar and vector particle attributes to the surface. To get rid of the typical bumps from SPH simulations, it supports a weighted Laplacian smoothing approach [detailed below](#weighted-surface-smoothing). As input, it supports reading particle positions from `.vtk`/`.vtu`, `.bgeo`, `.ply`, `.json` and binary `.xyz` (i.e. files containing a binary dump of a particle position array) files. Required parameters to perform a reconstruction are the kernel radius and particle radius (to compute the volume of particles) used for the original SPH simulation as well as the marching cubes resolution (a default iso-surface threshold is pre-configured). ## Domain decomposition A naive dense marching cubes reconstruction allocating a full 3D array over the entire fulid domain quickly becomes infeasible for larger simulations. Instead, one could use a global hashmap where only cubes that contain non-zero fluid density values are allocated. This approach is used in `splashsurf` if domain decomposition is disabled completely. However, the global hashmap approach does not lead to good cache locality and is not well suited for parallelization (even specialized parallel map implementations like [`dashmap`](https://github.com/xacrimon/dashmap) have their performance limitations). To improve on this situation `splashsurf` currently implements two domain decomposition approaches. ### Octree-based decomposition The octree-based decomposition is currently the default approach if no other option is specified but will probably be replaced by the grid-based approach described below. For the octree-based decomposition an octree is built over all particles with an automatically determined target number of particles per leaf node. For each leaf node, a hashmap is used like outlined above. As each hashmap is smaller, cache locality is improved and due to the decomposition, each thread can work on its own local hashmap. Finally, all surface patches are stitched together by walking the octree back up, resulting in a closed surface. Downsides of this approach are that the octree construction starting from the root and stitching back towards the root limit the amount of paralleism during some stages. ### Subdomain grid-based decomposition Since version 0.10.0, `splashsurf` implements a new domain decomposition approach called the "subdomain grid" approach, toggeled with the `--subdomain-grid=on` flag. Here, the goal is to divide the fluid domain into subdomains with a fixed number of marching cubes cells, by default `64x64x64` cubes. For each subdomain a dense 3D array is allocated for the marching cubes cells. Of course, only subdomains that contain fluid particles are actually allocated. For subdomains that contain only a very small number of fluid particles (less th 5% of the largest subdomain) a hashmap is used instead to not waste too much storage. As most domains are dense however, the marching cubes triangulation per subdomain is very fast as it can make full use of cache locality and the entire procedure is trivially parallelizable. For the stitching we ensure that we perform floating point operations in the same order at the subdomain boundaries (this can be ensured without synchronization). If the field values on the subdomain boundaries are identical from both sides, the marching cubes triangulations will be topologically compatible and can be merged in a post-processing step that is also parallelizable. Overall, this approach should almost always be faster than the previous octree-based aproach. ## Notes For small numbers of fluid particles (i.e. in the low thousands or less) the domain decomposition implementation may have worse performance due to the task based parallelism and the additional overhead of domain decomposition and stitching. In this case, you can try to disable the domain decomposition. The reconstruction will then use a global approach that is parallelized using thread-local hashmaps. For larger quantities of particles the decomposition approach is expected to be always faster. Due to the use of hash maps and multi-threading (if enabled), the output of this implementation is not deterministic. As shown below, the tool can handle the output of large simulations. However, it was not tested with a wide range of parameters and may not be totally robust against corner-cases or extreme parameters. If you experience problems, please report them together with your input data. ## Installation The command-line tool can be built from this repository but is also available on crates.io. If you have a [Rust toolchain installed](https://www.rust-lang.org/tools/install) you can install `splashsurf` with the command ``` cargo install splashsurf ``` ## Usage ### Recommended settings Good settings for the surface reconstruction depend on the original simulation and can be influenced by different conventions of different simulators. The following settings appear to work well with [SPlisHSPlasH](https://github.com/InteractiveComputerGraphics/SPlisHSPlasH): - `particle-radius`: should be a bit larger than the particle radius used for the actual simulation. A radius around 1.4 to 1.6 times larger than the original SPH particle radius seems to be appropriate. - `smoothing-length`: should be set around `1.2`. Larger values smooth out the iso-surface more but also artificially increase the fluid volume. - `surface-threshold`: a good value depends on the selected `particle-radius` and `smoothing-length` and can be used to counteract a fluid volume increase e.g. due to a larger particle radius. In combination with the other recommended values a threshold of `0.6` seemed to work well. - `cube-size` usually should not be chosen larger than `1.0` to avoid artifacts (e.g. single particles decaying into rhomboids), start with a value in the range of `0.75` to `0.5` and decrease/increase it if the result is too coarse or the reconstruction takes too long. ### Weighted surface smoothing The CLI implements the paper ["Weighted Laplacian Smoothing for Surface Reconstruction of Particle-based Fluids" (Löschner, Böttcher, Jeske, Bender; 2023)](https://animation.rwth-aachen.de/publication/0583/) which proposes a fast smoothing approach to avoid typical bumpy surfaces while preventing loss of volume that typically occurs with simple smoothing methods. The following images show a rendering of a typical surface reconstruction (on the right) with visible bumps due to the particles compared to the same surface reconstruction with weighted smoothing applied (on the left):
You can see this rendering in motion in [this video](https://youtu.be/2bYvaUXlBQs). To apply this smoothing, we recommend the following settings: - `--mesh-smoothing-weights=on`: This enables the use of special weights during the smoothing process that preserve fluid details. For more information we refer to the [paper](https://animation.rwth-aachen.de/publication/0583/). - `--mesh-smoothing-iters=25`: This enables smoothing of the output mesh. The individual iterations are relatively fast and 25 iterations appeared to strike a good balance between an initially bumpy surface and potential over-smoothing. - `--mesh-cleanup=on`/`--decimate-barnacles=on`: On of the options should be used when applying smoothing, otherwise artifacts can appear on the surface (for more details see the paper). The `mesh-cleanup` flag enables a general purpose marching cubes mesh cleanup procedure that removes small sliver triangles everywhere on the mesh. The `decimate-barnacles` enables a more targeted decimation that only removes specific triangle configurations that are problematic for the smoothing. The former approach results in a "nicer" mesh overall but can be slower than the latter. - `--normals-smoothing-iters=10`: If normals are being exported (with `--normals=on`), this results in an even smoother appearance during rendering. For the reconstruction parameters in conjunction with the weighted smoothing we recommend parameters close to the simulation parameters. That means selecting the same particle radius as in the simulation, a corresponding smoothing length (e.g. for SPlisHSPlasH a value of `2.0`), a surface-threshold between `0.6` and `0.7` and a cube size usually between `0.5` and `1.0`. A full invocation of the tool might look like this: ``` splashsurf reconstruct particles.vtk -r=0.025 -l=2.0 -c=0.5 -t=0.6 --subdomain-grid=on --mesh-cleanup=on --mesh-smoothing-weights=on --mesh-smoothing-iters=25 --normals=on --normals-smoothing-iters=10 ``` ### Benchmark example For example: ``` splashsurf reconstruct canyon_13353401_particles.xyz -r=0.011 -c=1.5 -l=2.0 -t=0.6 --subdomain-grid=on ``` With these parameters, a scene with 13353401 particles is reconstructed in less than 3 seconds on a Ryzen 9 5950X. The output is a mesh with 6069576 triangles. ``` [23:44:58.432][INFO] splashsurf v0.10.0 (splashsurf) [23:44:58.432][INFO] Called with command line: splashsurf reconstruct canyon_13353401_particles.xyz -r=0.011 -c=1.5 -l=2.0 -t=0.6 --subdomain-grid=on [23:44:58.432][INFO] Using single precision (f32) for surface reconstruction. [23:44:58.432][INFO] Reading particle dataset from "canyon_13353401_particles.xyz"... [23:44:58.515][INFO] Successfully read dataset with 13353401 particle positions. [23:44:58.520][INFO] Minimal enclosing bounding box of particles was computed as: AxisAlignedBoundingBox { min: [-25.0060978, -5.0146289, -40.0634613], max: [24.4994926, 18.3062096, 39.7757950] } [23:44:58.520][INFO] The ghost margin volume is 42.38% of the subdomain volume [23:44:58.520][INFO] The ghost margin is 3.03 MC cells or 0.05 subdomains thick [23:44:58.520][INFO] Number of subdomains: 82156 (47x23x76) [23:44:58.520][INFO] Number of MC cells per subdomain: 262144 (64x64x64) [23:44:58.520][INFO] Number of MC cells globally: 21536702464 (3008x1472x4864) [23:44:58.520][INFO] Starting classification of particles into subdomains. [23:44:58.601][INFO] Starting computation of global density vector. [23:44:59.548][INFO] Largest subdomain has 167861 particles. [23:44:59.548][INFO] Subdomains with 3358 or less particles will be considered sparse. [23:44:59.548][INFO] Starting reconstruction (level-set evaluation and local triangulation). [23:45:00.876][INFO] Starting stitching of subdomains to global mesh. [23:45:00.946][INFO] Global mesh has 3038116 vertices and 6069576 triangles. [23:45:00.996][INFO] Writing surface mesh to "canyon_surface.vtk"... [23:45:00.996][INFO] Writing mesh with 3038116 vertices and 6069576 cells to "canyon_surface.vtk"... [23:45:01.175][INFO] Successfully wrote mesh to file. [23:45:01.175][INFO] Done. [23:45:01.188][INFO] Successfully finished processing all inputs. [23:45:01.188][INFO] Timings: [23:45:01.188][INFO] reconstruct subcommand: 100.00%, 2756.58ms avg, 1 call (total: 2.757s) [23:45:01.188][INFO] surface reconstruction: 100.00%, 2756.54ms avg, 1 call (total: 2.757s) [23:45:01.188][INFO] loading particle positions: 3.00%, 82.68ms avg, 1 call (total: 0.083s) [23:45:01.188][INFO] compute minimum enclosing aabb: 0.21%, 5.91ms avg, 1 call (total: 0.006s) [23:45:01.188][INFO] surface reconstruction subdomain-grid: 88.17%, 2430.35ms avg, 1 call (total: 2.430s) [23:45:01.188][INFO] decomposition: 3.29%, 80.06ms avg, 1 call (total: 0.080s) [23:45:01.188][INFO] classifying particles: 21.22%, 16.99ms avg, 1 call (total: 0.017s) [23:45:01.188][INFO] merging TL per cell particle counters: 0.18%, 0.14ms avg, 1 call (total: 0.000s) [23:45:01.188][INFO] initializing flat subdomain data and index mapping: 0.08%, 0.06ms avg, 1 call (total: 0.000s) [23:45:01.188][INFO] copying particles to subdomains: 64.65%, 51.76ms avg, 1 call (total: 0.052s) [23:45:01.188][INFO] sort subdomain particles: 13.74%, 11.00ms avg, 1 call (total: 0.011s) [23:45:01.188][INFO] compute_global_density_vector: 39.00%, 947.82ms avg, 1 call (total: 0.948s) [23:45:01.188][INFO] subdomain density computation: ≈100.00%, 20.58ms avg, 1275 calls (total: 26.235s) [23:45:01.188][INFO] collect subdomain data: 0.67%, 0.14ms avg, 1275 calls (total: 0.175s) [23:45:01.188][INFO] initialize particle filter: 0.27%, 0.06ms avg, 1275 calls (total: 0.072s) [23:45:01.188][INFO] neighborhood_search_spatial_hashing_flat_filtered: 88.48%, 18.21ms avg, 1275 calls (total: 23.214s) [23:45:01.188][INFO] sequential_generate_cell_to_particle_map: 4.07%, 0.74ms avg, 1275 calls (total: 0.946s) [23:45:01.188][INFO] write particle neighbors: 95.09%, 17.31ms avg, 1275 calls (total: 22.074s) [23:45:01.188][INFO] sequential_compute_particle_densities_filtered: 10.30%, 2.12ms avg, 1275 calls (total: 2.702s) [23:45:01.188][INFO] update global density values: 0.26%, 0.05ms avg, 1275 calls (total: 0.068s) [23:45:01.188][INFO] reconstruction: 54.63%, 1327.61ms avg, 1 call (total: 1.328s) [23:45:01.188][INFO] subdomain reconstruction (sparse): ≈1.95%, 0.85ms avg, 899 calls (total: 0.768s) [23:45:01.188][INFO] density grid loop: 65.87%, 0.56ms avg, 899 calls (total: 0.506s) [23:45:01.188][INFO] mc triangulation loop: 25.02%, 0.21ms avg, 899 calls (total: 0.192s) [23:45:01.188][INFO] subdomain reconstruction (dense): ≈98.05%, 102.70ms avg, 376 calls (total: 38.617s) [23:45:01.188][INFO] density grid loop: 94.47%, 97.03ms avg, 376 calls (total: 36.482s) [23:45:01.188][INFO] mc triangulation loop: 5.19%, 5.33ms avg, 376 calls (total: 2.005s) [23:45:01.188][INFO] stitching: 2.84%, 69.04ms avg, 1 call (total: 0.069s) [23:45:01.188][INFO] surface patch offset scan: 0.01%, 0.01ms avg, 1 call (total: 0.000s) [23:45:01.188][INFO] copy interior verts/tris and deduplicate exterior verts: 83.51%, 57.65ms avg, 1 call (total: 0.058s) [23:45:01.188][INFO] write surface mesh to file: 6.50%, 179.17ms avg, 1 call (total: 0.179s) [23:45:01.188][INFO] writing mesh: 99.98%, 179.14ms avg, 1 call (total: 0.179s) ``` ### Sequences of files You can either process a single file or let the tool automatically process a sequence of files. A sequence of files is indicated by specifying a filename with a `{}` placeholder pattern in the name. The tool will treat the placeholder as a `(\d+)` regex, i.e. a group matching to at least one digit. This allows for any zero padding as well as non-zero padded incrementing indices. All files in the input path matching this pattern will then be processed in natural sort order (i.e. silently skipping missing files in the sequence). Note that the tool collects all existing filenames as soon as the command is invoked and does not update the list while running. The first and last file of a sequences that should be processed can be specified with the `-s`/`--start-index` and/or `-e`/`--end-index` arguments. By specifying the flag `--mt-files=on`, several files can be processed in parallel. If this is enabled, you should also set `--mt-particles=off` as enabling both will probably degrade performance. The combination of `--mt-files=on` and `--mt-particles=off` can be faster if many files with only few particles have to be processed. The number of threads can be influenced using the `--num-threads`/`-n` argument or the `RAYON_NUM_THREADS` environment variable **NOTE:** Currently, some functions do not have a sequential implementation and always parallelize over the particles or the mesh/domain. This includes: - the new "subdomain-grid" domain decomposition approach, as an alternative to the previous octree-based approach - some post-processing functionality (interpolation of smoothing weights, interpolation of normals & other fluid attributes) Using the `--mt-particles=off` argument does not have an effect on these parts of the surface reconstruction. For now, it is therefore recommended to not parallelize over multiple files if this functionality is used. ## Input file formats ### VTK Legacy VTK files with the "`.vtk`" extension are loaded using [`vtkio`](https://crates.io/crates/vtkio). The VTK file is loaded as a big endian binary file and has to contain an "Unstructured Grid" with either `f32` or `f64` vertex coordinates. Any other data or attributes are ignored except for those attributes that were specified with the ` --interpolate-attributes` command line argument. Currently supported attribute data types are scalar integers, floats and three-component float vectors. Only the first "Unstructured Grid" is loaded, other entities are ignored. Not that currently only the "pure" v4.2 legacy format is supported as documented on [here](https://kitware.github.io/vtk-examples/site/VTKFileFormats/#simple-legacy-formats). This corresponds to the `--output-format vtk42` flag of the [`meshio convert`](https://github.com/nschloe/meshio) tool. ### VTU VTK XML files with the "`.vtu`" extension are loaded using [`vtkio`](https://crates.io/crates/vtkio). Currently only VTU files using ASCII or encoded binary are supported. Files using "raw" binary sections (i.e. a `