Fix the marching-cubes mesh extraction from the poisson reconstruction

CHANGELOG.md

@@ -0,0 +1,8 @@
+# Changelog
+
+## v0.3.1
+
+- Fix the extraction of a mesh from the poisson reconstruction.
+- Add `PoissonReconstruction::reconstruct_trimesh` and `::reconstruct_mesh_buffers` for extracting
+  a triangle mesh with properly wired-up topology.
+

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "poisson_reconstruction"
 repository = "https://github.com/ForesightMiningSoftwareCorporation/PoissonReconstruction"
-version = "0.3.0"
+version = "0.3.1"
 license = "MIT OR Apache-2.0"
 description = "Screened Poisson Reconstruction algorithm in Rust"
 authors = ["Sébastien Crozet <sebcrozet@dimforge.com>"]
@@ -22,5 +22,6 @@ itertools = "0.10"
 fnv = "1"
 
 [dev-dependencies]
-bevy = "0.9"
+bevy = "0.15"
+bevy_panorbit_camera = "0.22"
 ply-rs = "0.1"

examples/reconstruction_demo.rs

@@ -1,19 +1,22 @@
+use bevy::asset::RenderAssetUsages;
 use bevy::pbr::wireframe::{Wireframe, WireframePlugin};
 use bevy::prelude::*;
-use bevy::render::mesh::PrimitiveTopology;
+use bevy::render::mesh::{Indices, PrimitiveTopology};
+use bevy_panorbit_camera::{PanOrbitCamera, PanOrbitCameraPlugin};
 use nalgebra::{Point3, Vector3};
 use ply_rs::{parser, ply};
+use poisson_reconstruction::marching_cubes::MeshBuffers;
 use poisson_reconstruction::{PoissonReconstruction, Real};
 use std::io::BufRead;
 use std::path::Path;
 use std::str::FromStr;
 
 fn main() {
     App::new()
-        .add_plugins(DefaultPlugins)
-        .add_plugin(WireframePlugin)
-        .add_startup_system(setup_camera_and_light)
-        .add_startup_system(setup_scene)
+        .add_plugins((DefaultPlugins, PanOrbitCameraPlugin))
+        .add_plugins(WireframePlugin)
+        .add_systems(Startup, setup_camera_and_light)
+        .add_systems(Startup, setup_scene)
         .run();
 }
 
@@ -38,32 +41,40 @@ fn setup_camera_and_light(mut commands: Commands) {
         transform: Transform::from_xyz(-100.0, 50.0, -100.0),
         ..default()
     });
-    commands.spawn(Camera3dBundle {
-        transform: Transform::from_xyz(-100.0, 25.0, -100.0)
-            .looking_at(Vec3::new(0., 0., 0.), Vec3::Y),
-        ..default()
-    });
+    commands.spawn((
+        Camera3dBundle {
+            transform: Transform::from_xyz(-100.0, 25.0, -100.0)
+                .looking_at(Vec3::new(0., 0., 0.), Vec3::Y),
+            ..default()
+        },
+        PanOrbitCamera::default(),
+    ));
 }
 
 fn spawn_mesh(
     commands: &mut Commands,
     meshes: &mut Assets<Mesh>,
     materials: &mut Assets<StandardMaterial>,
-    points: Vec<Point3<f64>>,
+    points: MeshBuffers,
 ) {
     // Create the bevy mesh.
     let vertices: Vec<_> = points
+        .vertices()
         .iter()
-        .map(|pt| [pt.x as f32, -pt.y as f32, pt.z as f32])
+        .map(|pt| [pt.x as f32, pt.y as f32, pt.z as f32])
         .collect();
-    let mut mesh = Mesh::new(PrimitiveTopology::TriangleList);
+    let mut mesh = Mesh::new(
+        PrimitiveTopology::TriangleList,
+        RenderAssetUsages::default(),
+    );
     mesh.insert_attribute(Mesh::ATTRIBUTE_POSITION, vertices);
-    mesh.compute_flat_normals();
+    mesh.insert_indices(Indices::U32(points.indices().to_vec()));
 
     commands
         .spawn(PbrBundle {
-            mesh: meshes.add(mesh),
-            material: materials.add(Color::rgb(1.0, 1.0, 0.0).into()),
+            mesh: Mesh3d(meshes.add(mesh)),
+            material: MeshMaterial3d(materials.add(Color::srgb(0.0, 1.0, 0.0))),
+            transform: Transform::from_rotation(Quat::from_rotation_x(180.0f32.to_radians())),
             ..default()
         })
         .insert(Wireframe);
@@ -133,12 +144,12 @@ fn parse_file(path: impl AsRef<Path>, ply: bool) -> Vec<VertexWithNormal> {
     }
 }
 
-fn reconstruct_surface(vertices: &[VertexWithNormal]) -> Vec<Point3<Real>> {
+fn reconstruct_surface(vertices: &[VertexWithNormal]) -> MeshBuffers {
     let points: Vec<_> = vertices.iter().map(|v| v.pos).collect();
     let normals: Vec<_> = vertices.iter().map(|v| v.normal).collect();
 
     dbg!("Running poisson.");
     let poisson = PoissonReconstruction::from_points_and_normals(&points, &normals, 0.0, 6, 6, 10);
     dbg!("Extracting vertices.");
-    poisson.reconstruct_mesh()
+    poisson.reconstruct_mesh_buffers()
 }

src/marching_cubes.rs

@@ -2,6 +2,52 @@
 
 use crate::Real;
 use na::Point3;
+use parry::bounding_volume::Aabb;
+use parry::shape::{TriMesh, TriMeshFlags};
+use parry::utils::SortedPair;
+use std::collections::HashMap;
+
+type MarchingCubesCellKey = [i32; 3];
+
+/// Represents an index and vertex buffer of a mesh for incremental construction.
+#[derive(Default)]
+pub struct MeshBuffers {
+    vertices: Vec<Point3<f64>>,
+    indices: Vec<u32>,
+    edge_to_index: HashMap<SortedPair<MarchingCubesCellKey>, u32>,
+}
+
+impl MeshBuffers {
+    /// The mesh’s index buffer.
+    pub fn indices(&self) -> &[u32] {
+        &self.indices
+    }
+
+    /// The mesh’s vertex buffer.
+    pub fn vertices(&self) -> &[Point3<f64>] {
+        &self.vertices
+    }
+
+    /// Return the results as a soup of triangle, with duplicated vertices.
+    pub fn result_as_triangle_soup(&self) -> Vec<Point3<f64>> {
+        self.indices
+            .iter()
+            .map(|i| self.vertices[*i as usize])
+            .collect()
+    }
+
+    /// Constructs a `TriMesh` from this buffer.
+    ///
+    /// The result is `None` if the index buffer of `self` is empty.
+    pub fn result(&self, flags: TriMeshFlags) -> Option<TriMesh> {
+        let idx: Vec<_> = self
+            .indices
+            .chunks_exact(3)
+            .map(|i| [i[0], i[1], i[2]])
+            .collect();
+        (!idx.is_empty()).then(|| TriMesh::with_flags(self.vertices.clone(), idx, flags))
+    }
+}
 
 /* The cube vertex and edge indices for base rotation:
  *
@@ -104,6 +150,69 @@ pub fn march_cube(
     }
 }
 
+// The triangle table gives us the mapping from index to actual
+// triangles to return for this configuration
+// v0 assumed at 0.0, 0.0, 0.0 & v6 at 1.0, 1.0, 1.0
+pub(crate) fn march_cube_idx(
+    aabb: &Aabb,
+    corner_values: &[f64; 8],
+    // Grid coordinates of v0.
+    first_corner_cell_key: [i32; 3],
+    iso_value: f64,
+    out: &mut MeshBuffers,
+) {
+    // Compute the index for MC_TRI_TABLE
+    let mut index = 0;
+    let old_indices_len = out.indices.len();
+
+    for (v, value) in corner_values.iter().enumerate() {
+        if *value < iso_value {
+            index |= 1 << v;
+        }
+    }
+
+    for t in MC_TRI_TABLE[index].iter().take_while(|t| **t >= 0) {
+        let v_idx = *t as usize;
+        let [v0, v1] = EDGE_VERTICES[v_idx];
+
+        let local_corner_0 = INDEX_TO_VERTEX[v0 as usize];
+        let local_corner_1 = INDEX_TO_VERTEX[v1 as usize];
+
+        let eid0 = [
+            first_corner_cell_key[0] + local_corner_0[0] as i32,
+            first_corner_cell_key[1] + local_corner_0[1] as i32,
+            first_corner_cell_key[2] + local_corner_0[2] as i32,
+        ];
+        let eid1 = [
+            first_corner_cell_key[0] + local_corner_1[0] as i32,
+            first_corner_cell_key[1] + local_corner_1[1] as i32,
+            first_corner_cell_key[2] + local_corner_1[2] as i32,
+        ];
+
+        let edge_key = SortedPair::new(eid0, eid1);
+        let vid = *out.edge_to_index.entry(edge_key).or_insert_with(|| {
+            // The normalized_vert will have components
+            // in 0..1.
+            let normalized_vert = lerp_vertices(
+                &INDEX_TO_VERTEX[v0 as usize],
+                &INDEX_TO_VERTEX[v1 as usize],
+                corner_values[v0 as usize],
+                corner_values[v1 as usize],
+                iso_value,
+            );
+
+            // Convert the normalized_vert into an Aabb vert.
+            let vertex = aabb.mins + aabb.extents().component_mul(&normalized_vert.coords);
+            out.vertices.push(vertex);
+            (out.vertices.len() - 1) as u32
+        });
+
+        out.indices.push(vid);
+    }
+
+    out.indices[old_indices_len..].reverse();
+}
+
 /// Interpolates linearly between two weighted integer points.
 ///
 /// # Parameters

src/poisson.rs

@@ -1,12 +1,13 @@
 use crate::hgrid::HGrid;
-use crate::marching_cubes::march_cube;
+use crate::marching_cubes::{march_cube_idx, MeshBuffers};
 use crate::poisson_layer::PoissonLayer;
 use crate::poisson_vector_field::PoissonVectorField;
 use crate::polynomial::{eval_bspline, eval_bspline_diff};
 use crate::Real;
 use na::{vector, Point3, Vector3};
 use parry::bounding_volume::{Aabb, BoundingVolume};
 use parry::partitioning::IndexedData;
+use parry::shape::{TriMesh, TriMeshFlags};
 use std::collections::HashMap;
 use std::ops::{AddAssign, Mul};
 
@@ -167,23 +168,83 @@ impl PoissonReconstruction {
 
     /// Reconstructs a mesh from this implicit function using a simple marching-cubes, extracting
     /// the isosurface at 0.
+    #[deprecated = "use `reconstruct_mesh_buffers` or `reconstruct_trimesh` instead"]
     pub fn reconstruct_mesh(&self) -> Vec<Point3<Real>> {
-        let mut vertices = vec![];
+        self.reconstruct_mesh_buffers().result_as_triangle_soup()
+    }
+
+    /// Reconstructs a `TriMesh` from this implicit function using a simple marching-cubes, extracting
+    /// the isosurface at 0.
+    pub fn reconstruct_trimesh(&self, flags: TriMeshFlags) -> Option<TriMesh> {
+        self.reconstruct_mesh_buffers().result(flags)
+    }
+
+    /// Reconstructs a mesh from this implicit function using a simple marching-cubes, extracting
+    /// the isosurface at 0.
+    pub fn reconstruct_mesh_buffers(&self) -> MeshBuffers {
+        let mut result = MeshBuffers::default();
+        let mut visited = HashMap::new();
 
         if let Some(last_layer) = self.layers.last() {
-            for cell in last_layer.cells_qbvh.raw_proxies() {
-                let aabb = last_layer.cells_qbvh.node_aabb(cell.node).unwrap();
+            // Check all the existing leaves.
+            let mut eval_cell = |key: Point3<i64>, visited: &mut HashMap<Point3<i64>, bool>| {
+                let cell_center = last_layer.grid.cell_center(&key);
+                let cell_width = Vector3::repeat(last_layer.grid.cell_width() / 2.0);
+                let aabb = Aabb::from_half_extents(cell_center, cell_width);
                 let mut vertex_values = [0.0; 8];
 
                 for (pt, val) in aabb.vertices().iter().zip(vertex_values.iter_mut()) {
                     *val = self.eval(pt);
                 }
 
-                march_cube(&aabb.mins, &aabb.maxs, &vertex_values, 0.0, &mut vertices);
+                let len_before = result.indices().len();
+                march_cube_idx(
+                    &aabb,
+                    &vertex_values,
+                    key.cast::<i32>().into(),
+                    0.0,
+                    &mut result,
+                );
+                let has_sign_change = result.indices().len() != len_before;
+                visited.insert(key, has_sign_change);
+                has_sign_change
+            };
+
+            for cell in last_layer.cells_qbvh.raw_proxies() {
+                // let aabb = last_layer.cells_qbvh.node_aabb(cell.node).unwrap();
+                eval_cell(cell.data.cell, &mut visited);
+            }
+
+            // Checking only the leaves isn’t enough, isosurfaces might escape leaves through levels
+            // at a coarser level. So we also check adjacent leaves that experienced a sign change.
+            // PERF: instead of traversing ALL the adjacent leaves, only traverse the ones adjacent
+            //       to an edge that actually experienced a sign change.
+            // PERF: don’t re-evaluate vertices that were already evaluated.
+            let mut stack: Vec<_> = visited
+                .iter()
+                .filter(|(_key, sign_change)| **sign_change)
+                .map(|e| *e.0)
+                .collect();
+
+            while let Some(cell) = stack.pop() {
+                for i in -1..=1 {
+                    for j in -1..=1 {
+                        for k in -1..=1 {
+                            let new_cell = cell + Vector3::new(i, j, k);
+
+                            if !visited.contains_key(&new_cell) {
+                                let has_sign_change = eval_cell(new_cell, &mut visited);
+                                if has_sign_change {
+                                    stack.push(new_cell);
+                                }
+                            }
+                        }
+                    }
+                }
             }
         }
 
-        vertices
+        result
     }
 }