Fix bug with edges going from above threshold to outside support

See #215

splashsurf_lib/src/density_map.rs

@@ -528,7 +528,7 @@ pub(crate) fn compute_kernel_evaluation_radius<I: Index, R: Real>(
     // The number of points corresponding to the number of supported cells
     let supported_points: I = I::one() + supported_cells;
 
-    // Evaluate kernel in a smaller domain, points outside of this radius have to be assumed to be outside of the iso-surface
+    // Evaluate kernel in a smaller domain, points outside of this radius have to be assumed to be outside the iso-surface
     let kernel_evaluation_radius =
         cube_size * half_supported_cells_real * (R::one() + R::default_epsilon().sqrt());
 
@@ -599,7 +599,7 @@ impl<I: Index, R: Real> SparseDensityMapGenerator<I, R> {
         particle: &Vector3<R>,
         particle_density: R,
     ) {
-        // Skip particles outside of allowed domain
+        // Skip particles outside allowed domain
         if !self.allowed_domain.contains_point(particle) {
             return;
         }

splashsurf_lib/src/lib.rs

@@ -412,7 +412,7 @@ pub fn grid_for_reconstruction<I: Index, R: Real>(
             } else {
                 Aabb3d::from_points(particle_positions)
             };
-            // TODO: Is this really necessary?
+            // TODO: Is this really necessary? This seems unnecessary purely for the density map...
             aabb.grow_uniformly(particle_radius);
             aabb
         };

splashsurf_lib/src/marching_cubes/narrow_band_extraction.rs

@@ -45,7 +45,7 @@ pub(crate) fn construct_mc_input<I: Index, R: Real>(
 /// Note: The threshold flags in the resulting cell data are not complete and still have to be updated after
 /// this procedure using the [update_cell_data_threshold_flags] function.
 ///
-/// Note: This functions assumes that the default value for missing point data is below the iso-surface threshold.
+/// Note: This functions assumes that the default value for missing point data is zero.
 #[inline(never)]
 fn interpolate_points_to_cell_data_generic<I: Index, R: Real>(
     grid: &UniformGrid<I, R>,
@@ -66,15 +66,13 @@ fn interpolate_points_to_cell_data_generic<I: Index, R: Real>(
         density_map.for_each(|flat_point_index, point_value| {
             let point = grid.try_unflatten_point_index(flat_point_index).unwrap();
 
-            // We want to find edges that cross the iso-surface,
-            // therefore we can choose to either skip all points above or below the threshold.
-            //
-            // In most scenes, the sparse density map should contain more entries above than
-            // below the threshold, as it contains the whole fluid interior, whereas areas completely
-            // devoid of fluid are not part of the density map.
-            //
-            // Skip points with densities above the threshold to improve efficiency
-            if point_value > iso_surface_threshold {
+            // We want to find edges that cross the iso-surface, i.e. edges where one point is above
+            // and the other below the threshold. To not process an edge twice, we skip all points
+            // that are already below the iso-surface threshold. Note that we cannot do it the other
+            // way around, as some points below the threshold might not be part of the density map
+            // at all (e.g. points outside the kernel evaluation radius). This could lead to missing
+            // edges that go directly from above the threshold to e.g. zero.
+            if point_value < iso_surface_threshold {
                 return;
             }
 
@@ -89,15 +87,13 @@ fn interpolate_points_to_cell_data_generic<I: Index, R: Real>(
                 let neighbor_value = if let Some(v) = density_map.get(flat_neighbor_index) {
                     v
                 } else {
-                    // Neighbors that are not in the point-value map were outside of the kernel evaluation radius.
-                    // This should only happen for cells that are completely outside of the compact support of a particle.
-                    // The point-value map has to be consistent such that for each cell, where at least one point-value
-                    // is missing like this, the cell has to be completely below the iso-surface threshold.
-                    continue;
+                    // Neighbors that are not in the point-value map were outside the kernel evaluation radius.
+                    // Assume zero density for these points.
+                    R::zero()
                 };
 
                 // Skip edges that don't cross the iso-surface
-                if !(neighbor_value > iso_surface_threshold) {
+                if !(neighbor_value < iso_surface_threshold) {
                     continue;
                 }
 
@@ -132,8 +128,8 @@ fn interpolate_points_to_cell_data_generic<I: Index, R: Real>(
                     );
                     cell_data_entry.iso_surface_vertices[local_edge_index] = Some(vertex_index);
 
-                    // Mark the neighbor as above the iso-surface threshold
-                    let local_vertex_index = cell.local_point_index_of(neighbor.index()).unwrap();
+                    // Mark the corner of the current point as above the iso-surface threshold
+                    let local_vertex_index = cell.local_point_index_of(point.index()).unwrap();
                     cell_data_entry.corner_above_threshold[local_vertex_index] =
                         RelativeToThreshold::Above;
                 }

splashsurf_lib/tests/integration_tests/mod.rs

@@ -3,3 +3,4 @@ pub mod test_full;
 #[cfg(feature = "io")]
 pub mod test_mesh;
 pub mod test_neighborhood_search;
+pub mod test_simple;

splashsurf_lib/tests/integration_tests/test_simple.rs

@@ -0,0 +1,141 @@
+use nalgebra::Vector3;
+use splashsurf_lib::marching_cubes::check_mesh_consistency;
+use splashsurf_lib::{
+    reconstruct_surface, Aabb3d, GridDecompositionParameters, Parameters, Real,
+    SpatialDecomposition,
+};
+
+enum Strategy {
+    Global,
+    SubdomainGrid,
+}
+
+fn params_with_aabb<R: Real>(
+    particle_radius: R,
+    compact_support_radius: R,
+    cube_size: R,
+    iso_surface_threshold: R,
+    domain_aabb: Option<Aabb3d<R>>,
+    strategy: Strategy,
+) -> Parameters<R> {
+    let compact_support_radius = particle_radius * compact_support_radius;
+    let cube_size = particle_radius * cube_size;
+
+    let mut parameters = Parameters {
+        particle_radius,
+        rest_density: R::from_f64(1000.0).unwrap(),
+        compact_support_radius,
+        cube_size,
+        iso_surface_threshold,
+        particle_aabb: domain_aabb,
+        enable_multi_threading: false,
+        spatial_decomposition: None,
+        global_neighborhood_list: false,
+    };
+
+    match strategy {
+        Strategy::Global => {}
+        Strategy::SubdomainGrid => {
+            parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                GridDecompositionParameters {
+                    subdomain_num_cubes_per_dim: 64,
+                },
+            ))
+        }
+    }
+
+    parameters
+}
+
+fn params<R: Real>(
+    particle_radius: R,
+    compact_support_radius: R,
+    cube_size: R,
+    iso_surface_threshold: R,
+    strategy: Strategy,
+) -> Parameters<R> {
+    params_with_aabb(
+        particle_radius,
+        compact_support_radius,
+        cube_size,
+        iso_surface_threshold,
+        None,
+        strategy,
+    )
+}
+
+/// This tests ensures that a surface is reconstructed properly if a single edge is above the threshold
+/// on one side but outside the compact support on the other side.
+#[test]
+fn test_edge_above_threshold_to_outside_of_compact_support_global() {
+    let particle_positions: Vec<Vector3<f32>> = vec![Vector3::new(0.01, 0.0, 0.0)];
+    let parameters = params(1.0, 1.0, 1.0, 0.1, Strategy::Global);
+    let reconstruction =
+        reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap();
+
+    println!(
+        "Reconstructed mesh from {} particles has {} triangles.",
+        particle_positions.len(),
+        reconstruction.mesh().triangles.len()
+    );
+
+    assert_eq!(
+        reconstruction.mesh().vertices.len(),
+        6,
+        "Number of vertices"
+    );
+    assert_eq!(
+        reconstruction.mesh().triangles.len(),
+        8,
+        "Number of triangles"
+    );
+
+    // Ensure that the mesh does not have a boundary
+    if let Err(e) = check_mesh_consistency(
+        reconstruction.grid(),
+        reconstruction.mesh(),
+        true,
+        true,
+        true,
+    ) {
+        eprintln!("{}", e);
+        panic!("Mesh contains topological/manifold errors");
+    }
+}
+
+#[test]
+fn test_edge_above_threshold_to_outside_of_compact_support_subdomains() {
+    let particle_positions: Vec<Vector3<f32>> = vec![Vector3::new(0.01, 0.0, 0.0)];
+    let parameters = params(1.0, 1.0, 1.0, 0.1, Strategy::SubdomainGrid);
+    let reconstruction =
+        reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap();
+
+    println!(
+        "Reconstructed mesh from {} particles has {} triangles.",
+        particle_positions.len(),
+        reconstruction.mesh().triangles.len()
+    );
+
+    assert_eq!(
+        reconstruction.mesh().vertices.len(),
+        6,
+        "Number of vertices"
+    );
+    assert_eq!(
+        reconstruction.mesh().triangles.len(),
+        8,
+        "Number of triangles"
+    );
+
+    // Ensure that the mesh does not have a boundary
+    if let Err(e) = check_mesh_consistency(
+        reconstruction.grid(),
+        reconstruction.mesh(),
+        true,
+        true,
+        true,
+    ) {
+        eprintln!("{}", e);
+        panic!("Mesh contains topological/manifold errors");
+    }
+}