Add function and CLI flag to check mesh for holes

Author: w1th0utnam3

splashsurf/src/reconstruction.rs

@@ -3,7 +3,7 @@ use anyhow::{anyhow, Context};
 use arguments::{
     ReconstructionRunnerArgs, ReconstructionRunnerPathCollection, ReconstructionRunnerPaths,
 };
-use log::info;
+use log::{error, info};
 use rayon::prelude::*;
 use splashsurf_lib::coarse_prof::profile;
 use splashsurf_lib::mesh::PointCloud3d;
@@ -95,6 +95,9 @@ pub struct ReconstructSubcommandArgs {
     /// Set the number of threads for the worker thread pool
     #[structopt(long, short = "-n")]
     num_threads: Option<usize>,
+    /// Whether to check the final mesh for problems such as holes (note that when stitching is disabled this will lead to a lot of reported problems)
+    #[structopt(long, default_value = "off", possible_values = &["on", "off"], case_insensitive = true)]
+    check_mesh: Switch,
 }
 
 arg_enum! {
@@ -124,7 +127,7 @@ pub fn reconstruct_subcommand(cmd_args: &ReconstructSubcommandArgs) -> Result<()
 
     let result = if cmd_args.parallelize_over_files.into_bool() {
         paths.par_iter().try_for_each(|path| {
-            reconstruction_entry_point(path, &args)
+            reconstruction_pipeline(path, &args)
                 .with_context(|| {
                     format!(
                         "Error while processing input file '{}' from a file sequence",
@@ -140,7 +143,7 @@ pub fn reconstruct_subcommand(cmd_args: &ReconstructSubcommandArgs) -> Result<()
     } else {
         paths
             .iter()
-            .try_for_each(|path| reconstruction_entry_point(path, &args))
+            .try_for_each(|path| reconstruction_pipeline(path, &args))
     };
 
     if result.is_ok() {
@@ -166,6 +169,7 @@ mod arguments {
     pub struct ReconstructionRunnerArgs {
         pub params: splashsurf_lib::Parameters<f64>,
         pub use_double_precision: bool,
+        pub check_mesh: bool,
         pub io_params: io::FormatParameters,
     }
 
@@ -239,6 +243,7 @@ mod arguments {
             Ok(ReconstructionRunnerArgs {
                 params,
                 use_double_precision: args.use_double_precision.into_bool(),
+                check_mesh: args.check_mesh.into_bool(),
                 io_params: io::FormatParameters::default(),
             })
         }
@@ -473,32 +478,39 @@ mod arguments {
 }
 
 /// Calls the reconstruction pipeline for single or double precision depending on the runtime parameters
-pub(crate) fn reconstruction_entry_point(
+pub(crate) fn reconstruction_pipeline(
     paths: &ReconstructionRunnerPaths,
     args: &ReconstructionRunnerArgs,
 ) -> Result<(), anyhow::Error> {
     if args.use_double_precision {
         info!("Using double precision (f64) for surface reconstruction.");
-        reconstruction_entry_point_generic::<i64, f64>(paths, &args.params, &args.io_params)?;
+        reconstruction_pipeline_generic::<i64, f64>(
+            paths,
+            &args.params,
+            &args.io_params,
+            args.check_mesh,
+        )?;
     } else {
         info!("Using single precision (f32) for surface reconstruction.");
-        reconstruction_entry_point_generic::<i64, f32>(
+        reconstruction_pipeline_generic::<i64, f32>(
             paths,
             &args.params.try_convert().ok_or(anyhow!(
                 "Unable to convert surface reconstruction parameters from f64 to f32."
             ))?,
             &args.io_params,
+            args.check_mesh,
         )?;
     }
 
     Ok(())
 }
 
 /// Wrapper for the reconstruction pipeline: loads input file, runs reconstructions, stores output files
-pub(crate) fn reconstruction_entry_point_generic<I: Index, R: Real>(
+pub(crate) fn reconstruction_pipeline_generic<I: Index, R: Real>(
     paths: &ReconstructionRunnerPaths,
     params: &splashsurf_lib::Parameters<R>,
     io_params: &io::FormatParameters,
+    check_mesh: bool,
 ) -> Result<(), anyhow::Error> {
     profile!("surface reconstruction cli");
 
@@ -518,6 +530,14 @@ pub(crate) fn reconstruction_entry_point_generic<I: Index, R: Real>(
     let grid = reconstruction.grid();
     let mesh = reconstruction.mesh();
 
+    if check_mesh {
+        if let Err(err) = splashsurf_lib::marching_cubes::check_mesh_consistency(grid, mesh) {
+            error!("{}", err);
+        } else {
+            info!("Checked mesh for problems (holes, etc.), no problems were found.");
+        }
+    }
+
     // Store the surface mesh
     {
         profile!("write surface mesh to file");

splashsurf_lib/src/marching_cubes.rs

@@ -10,6 +10,7 @@ use crate::marching_cubes::triangulation::{
 use crate::mesh::TriMesh3d;
 use crate::uniform_grid::{DummySubdomain, OwningSubdomainGrid, Subdomain};
 use crate::{new_map, DensityMap, Index, MapType, Real, UniformGrid};
+use nalgebra::Vector3;
 
 pub mod marching_cubes_lut;
 mod narrow_band_extraction;
@@ -180,6 +181,87 @@ pub(crate) fn triangulate_density_map_to_surface_patch<I: Index, R: Real>(
     }
 }
 
+/// Checks the consistency of the mesh (currently only checks for holes) and returns a string with debug information in case of problems
+pub fn check_mesh_consistency<I: Index, R: Real>(
+    grid: &UniformGrid<I, R>,
+    mesh: &TriMesh3d<R>,
+) -> Result<(), String> {
+    let boundary_edges = mesh.find_boundary_edges();
+
+    if boundary_edges.is_empty() {
+        return Ok(());
+    }
+
+    let mut error_string = String::new();
+    error_string += &format!("Mesh is not closed. It has {} boundary edges (edges that are connected to only one triangle):", boundary_edges.len());
+    for (edge, tri_idx, _) in boundary_edges {
+        let v0 = mesh.vertices[edge[0]];
+        let v1 = mesh.vertices[edge[1]];
+        let center = (v0 + v1) / (R::one() + R::one());
+        let cell = grid.enclosing_cell(&center);
+        if let Some(cell_index) = grid.get_cell(cell) {
+            let point_index = grid
+                .get_point(*cell_index.index())
+                .expect("Unable to get point index of cell");
+            let cell_center = grid.point_coordinates(&point_index)
+                + &Vector3::repeat(grid.cell_size().times_f64(0.5));
+
+            error_string += &format!("\n\tTriangle {}, boundary edge {:?} is located in cell with {:?} with center coordinates {:?} and edge length {}.", tri_idx, edge, cell_index, cell_center, grid.cell_size());
+        } else {
+            error_string += &format!(
+                "\n\tCannot get cell index for boundary edge {:?} of triangle {}",
+                edge, tri_idx
+            );
+        }
+    }
+
+    Err(error_string)
+}
+
+/// Same as [`check_mesh_consistency`] but also adds debug information taken from the marching cubes input
+#[allow(unused)]
+fn check_mesh_with_cell_data<I: Index, R: Real>(
+    grid: &UniformGrid<I, R>,
+    marching_cubes_data: &MarchingCubesInput<I>,
+    mesh: &TriMesh3d<R>,
+) -> Result<(), String> {
+    let boundary_edges = mesh.find_boundary_edges();
+
+    if boundary_edges.is_empty() {
+        return Ok(());
+    }
+
+    let mut error_string = String::new();
+    error_string += &format!("Mesh is not closed. It has {} boundary edges (edges that are connected to only one triangle):", boundary_edges.len());
+    for (edge, tri_idx, _) in boundary_edges {
+        let v0 = mesh.vertices[edge[0]];
+        let v1 = mesh.vertices[edge[1]];
+        let center = (v0 + v1) / (R::one() + R::one());
+        let cell = grid.enclosing_cell(&center);
+        if let Some(cell_index) = grid.get_cell(cell) {
+            let point_index = grid
+                .get_point(*cell_index.index())
+                .expect("Unable to get point index of cell");
+            let cell_center = grid.point_coordinates(&point_index)
+                + &Vector3::repeat(grid.cell_size().times_f64(0.5));
+
+            let cell_data = marching_cubes_data
+                .cell_data
+                .get(&grid.flatten_cell_index(&cell_index))
+                .expect("Unabel to get cell data of cell");
+
+            error_string += &format!("\n\tTriangle {}, boundary edge {:?} is located in cell with {:?} with center coordinates {:?} and edge length {}. {:?}", tri_idx, edge, cell_index, cell_center, grid.cell_size(), cell_data);
+        } else {
+            error_string += &format!(
+                "\n\tCannot get cell index for boundary edge {:?} of triangle {}",
+                edge, tri_idx
+            );
+        }
+    }
+
+    Err(error_string)
+}
+
 #[allow(unused)]
 #[inline(never)]
 fn assert_cell_data_point_data_consistency<I: Index, R: Real>(

splashsurf_lib/tests/integration_tests/test_full.rs

@@ -1,4 +1,5 @@
 use nalgebra::Vector3;
+use splashsurf_lib::marching_cubes::check_mesh_consistency;
 use splashsurf_lib::{
     reconstruct_surface, Parameters, Real, SpatialDecompositionParameters, SubdivisionCriterion,
 };
@@ -107,11 +108,11 @@ macro_rules! generate_test {
 
             if test_for_boundary(&parameters) {
                 // Ensure that the mesh does not have a boundary
-                assert_eq!(
-                    reconstruction.mesh().find_boundary_edges(),
-                    vec![],
-                    "Mesh boundary is not empty"
-                );
+                if let Err(e) = check_mesh_consistency(reconstruction.grid(), reconstruction.mesh())
+                {
+                    eprintln!("{}", e);
+                    panic!("Mesh contains boundary edges");
+                }
             }
         }
     };