Make SurfaceReconstruction an inout param (#7)

Allow passing an existing `SurfaceReconstruction` to the reconstruction procedure to reuse memory allocated for the surface mesh.

Author: whiterabbit42k

splashsurf/src/reconstruction.rs

@@ -53,7 +53,9 @@ pub(crate) fn entry_point_generic<I: Index, R: Real>(
         splashsurf_lib::reconstruct_surface::<I, R>(particle_positions.as_slice(), &params)?;
 
     let grid = reconstruction.grid();
-    let density_map = reconstruction.density_map();
+    let density_map = reconstruction
+        .density_map()
+        .ok_or_else(|| anyhow::anyhow!("No density map was created during reconstruction"))?;
     let mesh = reconstruction.mesh();
 
     {

splashsurf_lib/src/lib.rs

@@ -146,20 +146,32 @@ pub struct SurfaceReconstruction<I: Index, R: Real> {
     /// The background grid that was used as a basis for generating the density map for marching cubes
     grid: UniformGrid<I, R>,
     /// The point-based density map generated from the particles that was used as input to marching cubes
-    density_map: DensityMap<I, R>,
+    density_map: Option<DensityMap<I, R>>,
     /// The actual mesh that is the result of the surface reconstruction
     mesh: TriMesh3d<R>,
 }
 
+impl<I: Index, R: Real> Default for SurfaceReconstruction<I, R> {
+    fn default() -> Self {
+        let grid = UniformGrid::new_zero();
+        let mesh = TriMesh3d::default();
+        Self {
+            grid,
+            density_map: None,
+            mesh,
+        }
+    }
+}
+
 impl<I: Index, R: Real> SurfaceReconstruction<I, R> {
     /// Returns a reference to the actual triangulated surface mesh that is the result of the reconstruction
     pub fn mesh(&self) -> &TriMesh3d<R> {
         &self.mesh
     }
 
     /// Returns a reference to the sparse density map (discretized on the vertices of the background grid) that is used as input for marching cubes
-    pub fn density_map(&self) -> &DensityMap<I, R> {
-        &self.density_map
+    pub fn density_map(&self) -> Option<&DensityMap<I, R>> {
+        self.density_map.as_ref()
     }
 
     /// Returns a reference to the virtual background grid that was used as a basis for discretization of the density map for marching cubes, can be used to convert the density map to a hex mesh (using [sparse_density_map_to_hex_mesh](density_map::sparse_density_map_to_hex_mesh))
@@ -205,6 +217,17 @@ pub fn reconstruct_surface<I: Index, R: Real>(
     parameters: &Parameters<R>,
 ) -> Result<SurfaceReconstruction<I, R>, ReconstructionError<I, R>> {
     profile!("reconstruct_surface");
+    let mut surface = SurfaceReconstruction::default();
+    reconstruct_surface_inplace(particle_positions, parameters, &mut surface)?;
+    Ok(surface)
+}
+
+pub fn reconstruct_surface_inplace<'a, I: Index, R: Real>(
+    particle_positions: &[Vector3<R>],
+    parameters: &Parameters<R>,
+    surface: &'a mut SurfaceReconstruction<I, R>,
+) -> Result<(), ReconstructionError<I, R>> {
+    profile!("reconstruct_surface_inplace");
 
     let Parameters {
         particle_radius,
@@ -217,12 +240,13 @@ pub fn reconstruct_surface<I: Index, R: Real>(
         enable_multi_threading,
     } = parameters.clone();
 
-    let grid = grid_for_reconstruction(
+    surface.grid = grid_for_reconstruction(
         particle_positions,
         particle_radius,
         cube_size,
         domain_aabb.as_ref(),
     )?;
+    let grid = &surface.grid;
 
     info!(
         "Using a grid with {:?}x{:?}x{:?} points and {:?}x{:?}x{:?} cells of edge length {}.",
@@ -291,14 +315,16 @@ pub fn reconstruct_surface<I: Index, R: Real>(
         enable_multi_threading,
     );
 
-    let mesh =
-        marching_cubes::triangulate_density_map::<I, R>(&grid, &density_map, iso_surface_threshold);
+    marching_cubes::triangulate_density_map::<I, R>(
+        &grid,
+        &density_map,
+        iso_surface_threshold,
+        &mut surface.mesh,
+    );
+
+    surface.density_map = Some(density_map);
 
-    Ok(SurfaceReconstruction {
-        grid,
-        density_map,
-        mesh,
-    })
+    Ok(())
 }
 
 /// Constructs the background grid for marching cubes based on the parameters supplied to the surface reconstruction

splashsurf_lib/src/marching_cubes.rs

@@ -1,5 +1,4 @@
 use log::info;
-use nalgebra::Vector3;
 
 use crate::marching_cubes_lut::get_marching_cubes_triangulation;
 use crate::mesh::TriMesh3d;
@@ -10,11 +9,12 @@ pub fn triangulate_density_map<I: Index, R: Real>(
     grid: &UniformGrid<I, R>,
     density_map: &DensityMap<I, R>,
     iso_surface_threshold: R,
-) -> TriMesh3d<R> {
+    mesh: &mut TriMesh3d<R>,
+) {
     profile!("triangulate_density_map");
     let marching_cubes_data =
-        interpolate_points_to_cell_data::<I, R>(&grid, &density_map, iso_surface_threshold);
-    triangulate::<I, R>(marching_cubes_data)
+        interpolate_points_to_cell_data::<I, R>(&grid, &density_map, iso_surface_threshold, mesh);
+    triangulate::<I, R>(marching_cubes_data, mesh)
 }
 
 /// Flag indicating whether a vertex is above or below the iso-surface
@@ -73,9 +73,7 @@ impl Default for CellData {
 
 /// Input for the marching cubes algorithm
 #[derive(Clone, Debug)]
-pub(crate) struct MarchingCubesInput<I: Index, R: Real> {
-    /// All interpolated vertices on the iso-surface
-    vertices: Vec<Vector3<R>>,
+pub(crate) struct MarchingCubesInput<I: Index> {
     /// Data for all cells that marching cubes has to visit
     cell_data: MapType<I, CellData>,
 }
@@ -91,7 +89,8 @@ pub(crate) fn interpolate_points_to_cell_data<I: Index, R: Real>(
     grid: &UniformGrid<I, R>,
     density_map: &DensityMap<I, R>,
     iso_surface_threshold: R,
-) -> MarchingCubesInput<I, R> {
+    mesh: &mut TriMesh3d<R>,
+) -> MarchingCubesInput<I> {
     profile!("interpolate_points_to_cell_data");
 
     // Note: This functions assumes that the default value for missing point data is below the iso-surface threshold
@@ -100,7 +99,8 @@ pub(crate) fn interpolate_points_to_cell_data<I: Index, R: Real>(
     // Map from flat cell index to all data that is required per cell for the marching cubes triangulation
     let mut cell_data: MapType<I, CellData> = new_map();
     // Storage for vertices that are created on edges crossing the iso-surface
-    let mut vertices = Vec::new();
+    let vertices = &mut mesh.vertices;
+    vertices.clear();
 
     // Generate iso-surface vertices and identify affected cells & edges
     {
@@ -227,31 +227,25 @@ pub(crate) fn interpolate_points_to_cell_data<I: Index, R: Real>(
     );
     info!("Interpolation done.");
 
-    MarchingCubesInput {
-        vertices,
-        cell_data,
-    }
+    MarchingCubesInput { cell_data }
 }
 
 /// Converts the marching cubes input cell data into a triangle surface mesh
 #[inline(never)]
-pub(crate) fn triangulate<I: Index, R: Real>(input: MarchingCubesInput<I, R>) -> TriMesh3d<R> {
+pub(crate) fn triangulate<I: Index, R: Real>(
+    input: MarchingCubesInput<I>,
+    mesh: &mut TriMesh3d<R>,
+) {
     profile!("triangulate");
 
-    let MarchingCubesInput {
-        vertices,
-        cell_data,
-    } = input;
+    let MarchingCubesInput { cell_data } = input;
 
     info!(
         "Starting marching cubes triangulation of {} cells...",
         cell_data.len()
     );
 
-    let mut mesh = TriMesh3d {
-        vertices,
-        triangles: Vec::new(),
-    };
+    mesh.triangles.clear();
 
     // Triangulate affected cells
     for (&_flat_cell_index, cell_data) in &cell_data {
@@ -273,8 +267,6 @@ pub(crate) fn triangulate<I: Index, R: Real>(input: MarchingCubesInput<I, R>) ->
         mesh.vertices.len()
     );
     info!("Triangulation done.");
-
-    mesh
 }
 
 #[allow(unused)]
@@ -328,9 +320,11 @@ fn assert_cell_data_point_data_consistency<I: Index, R: Real>(
 
 #[test]
 fn test_interpolate_cell_data() {
+    use nalgebra::Vector3;
     let iso_surface_threshold = 0.25;
     //let default_value = 0.0;
 
+    let mut trimesh = crate::TriMesh3d::default();
     let origin = Vector3::new(0.0, 0.0, 0.0);
     let n_cubes_per_dim = [1, 1, 1];
     let cube_size = 1.0;
@@ -341,10 +335,14 @@ fn test_interpolate_cell_data() {
 
     let mut sparse_data = new_map();
 
-    let marching_cubes_data =
-        interpolate_points_to_cell_data(&grid, &sparse_data.clone().into(), iso_surface_threshold);
+    let marching_cubes_data = interpolate_points_to_cell_data(
+        &grid,
+        &sparse_data.clone().into(),
+        iso_surface_threshold,
+        &mut trimesh,
+    );
 
-    assert_eq!(marching_cubes_data.vertices.len(), 0);
+    assert_eq!(trimesh.vertices.len(), 0);
     assert_eq!(marching_cubes_data.cell_data.len(), 0);
 
     let points = vec![
@@ -362,12 +360,16 @@ fn test_interpolate_cell_data() {
         sparse_data.insert(grid.flatten_point_index_array(&ijk), val);
     }
 
-    let marching_cubes_data =
-        interpolate_points_to_cell_data(&grid, &sparse_data.clone().into(), iso_surface_threshold);
+    let marching_cubes_data = interpolate_points_to_cell_data(
+        &grid,
+        &sparse_data.clone().into(),
+        iso_surface_threshold,
+        &mut trimesh,
+    );
 
     assert_eq!(marching_cubes_data.cell_data.len(), 1);
     // Check that the correct number of vertices was created
-    assert_eq!(marching_cubes_data.vertices.len(), 6);
+    assert_eq!(trimesh.vertices.len(), 6);
 
     let cell = &marching_cubes_data.cell_data[&0];
 
@@ -389,6 +391,6 @@ fn test_interpolate_cell_data() {
     assert!(cell.iso_surface_vertices[11].is_some());
 
     // TODO: Continue writing test
-    let _mesh = triangulate(marching_cubes_data);
+    let _mesh = triangulate(marching_cubes_data, &mut trimesh);
     //println!("{:?}", mesh)
 }

splashsurf_lib/src/mesh.rs

@@ -13,8 +13,17 @@ pub struct TriMesh3d<R: Real> {
     pub triangles: Vec<[usize; 3]>,
 }
 
+impl<R: Real> Default for TriMesh3d<R> {
+    fn default() -> Self {
+        Self {
+            vertices: Vec::new(),
+            triangles: Vec::new(),
+        }
+    }
+}
+
 /// A hexahedral (volumetric) mesh in 3D
-#[derive(Clone, Debug)]
+#[derive(Clone, Debug, Default)]
 pub struct HexMesh3d<R: Real> {
     /// Coordinates of all vertices of the mesh
     pub vertices: Vec<Vector3<R>>,
@@ -23,14 +32,14 @@ pub struct HexMesh3d<R: Real> {
 }
 
 /// A point cloud in 3D
-#[derive(Clone, Debug)]
+#[derive(Clone, Debug, Default)]
 pub struct PointCloud3d<R: Real> {
     /// Coordinates of all points in the point cloud
     pub points: Vec<Vector3<R>>,
 }
 
 /// A mesh with attached vertex or point data
-#[derive(Clone, Debug)]
+#[derive(Clone, Debug, Default)]
 pub struct MeshWithPointData<MeshT, DataT> {
     /// The mesh geometry itself
     pub mesh: MeshT,

splashsurf_lib/src/uniform_grid.rs

@@ -200,6 +200,16 @@ impl<I: Index, R: Real> UniformCartesianCubeGrid3d<I, R> {
         })
     }
 
+    /// Create a new zeroed grid
+    pub(crate) fn new_zero() -> Self {
+        Self {
+            aabb: AxisAlignedBoundingBox3d::new(Vector3::zeros(), Vector3::zeros()),
+            cell_size: R::zero(),
+            n_points_per_dim: [I::zero(); 3],
+            n_cells_per_dim: [I::zero(); 3],
+        }
+    }
+
     /// Returns the bounding box of the grid
     #[inline(always)]
     pub fn aabb(&self) -> &AxisAlignedBoundingBox3d<R> {