From 26699bbec6874145065cbe429e1d17b919b6d211 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20L=C3=B6schner?= Date: Thu, 27 Jul 2023 17:45:10 +0200 Subject: [PATCH] Implement filtering of particles according to AABB --- splashsurf/src/reconstruction.rs | 2 +- splashsurf_lib/benches/benches/bench_full.rs | 6 +- splashsurf_lib/benches/benches/bench_mesh.rs | 2 +- .../benches/benches/bench_subdomain_grid.rs | 2 +- splashsurf_lib/src/lib.rs | 88 ++++++++++++++----- splashsurf_lib/src/workspace.rs | 7 ++ .../tests/integration_tests/test_full.rs | 2 +- 7 files changed, 82 insertions(+), 27 deletions(-) diff --git a/splashsurf/src/reconstruction.rs b/splashsurf/src/reconstruction.rs index a93f459..33619ff 100644 --- a/splashsurf/src/reconstruction.rs +++ b/splashsurf/src/reconstruction.rs @@ -411,7 +411,7 @@ mod arguments { compact_support_radius, cube_size, iso_surface_threshold: args.surface_threshold, - domain_aabb, + particle_aabb: domain_aabb, enable_multi_threading: args.parallelize_over_particles.into_bool(), spatial_decomposition, }; diff --git a/splashsurf_lib/benches/benches/bench_full.rs b/splashsurf_lib/benches/benches/bench_full.rs index 3fe6573..1550c21 100644 --- a/splashsurf_lib/benches/benches/bench_full.rs +++ b/splashsurf_lib/benches/benches/bench_full.rs @@ -101,7 +101,7 @@ pub fn surface_reconstruction_dam_break(c: &mut Criterion) { compact_support_radius, cube_size, iso_surface_threshold: 0.6, - domain_aabb: None, + particle_aabb: None, enable_multi_threading: true, spatial_decomposition: None, }; @@ -199,7 +199,7 @@ pub fn surface_reconstruction_double_dam_break(c: &mut Criterion) { compact_support_radius, cube_size, iso_surface_threshold: 0.6, - domain_aabb: None, + particle_aabb: None, enable_multi_threading: true, spatial_decomposition: None, }; @@ -297,7 +297,7 @@ pub fn surface_reconstruction_double_dam_break_inplace(c: &mut Criterion) { compact_support_radius, cube_size, iso_surface_threshold: 0.6, - domain_aabb: None, + particle_aabb: None, enable_multi_threading: true, spatial_decomposition: None, }; diff --git a/splashsurf_lib/benches/benches/bench_mesh.rs b/splashsurf_lib/benches/benches/bench_mesh.rs index baff8d6..36951e2 100644 --- a/splashsurf_lib/benches/benches/bench_mesh.rs +++ b/splashsurf_lib/benches/benches/bench_mesh.rs @@ -22,7 +22,7 @@ fn reconstruct_particles>(particle_file: P) -> SurfaceReconstruct compact_support_radius, cube_size, iso_surface_threshold: 0.6, - domain_aabb: None, + particle_aabb: None, enable_multi_threading: true, spatial_decomposition: Some(SpatialDecomposition::Octree( OctreeDecompositionParameters { diff --git a/splashsurf_lib/benches/benches/bench_subdomain_grid.rs b/splashsurf_lib/benches/benches/bench_subdomain_grid.rs index e14ec3f..24dec38 100644 --- a/splashsurf_lib/benches/benches/bench_subdomain_grid.rs +++ b/splashsurf_lib/benches/benches/bench_subdomain_grid.rs @@ -18,7 +18,7 @@ fn parameters_canyon() -> Parameters { compact_support_radius, cube_size, iso_surface_threshold: 0.6, - domain_aabb: None, + particle_aabb: None, enable_multi_threading: true, spatial_decomposition: Some(SpatialDecomposition::UniformGrid( GridDecompositionParameters { diff --git a/splashsurf_lib/src/lib.rs b/splashsurf_lib/src/lib.rs index 5d3fc21..1c53368 100644 --- a/splashsurf_lib/src/lib.rs +++ b/splashsurf_lib/src/lib.rs @@ -27,6 +27,7 @@ use log::info; /// Re-export the version of `nalgebra` used by this crate pub use nalgebra; use nalgebra::Vector3; +use std::borrow::Cow; use std::hash::Hash; use thiserror::Error as ThisError; /// Re-export the version of `vtkio` used by this crate, if vtk support is enabled @@ -258,7 +259,7 @@ pub struct Parameters { /// The surface reconstruction always results in a closed mesh around the particles. /// The final mesh can extend beyond this AABB due to the smoothing of the kernel. /// If not provided, the smallest AABB enclosing all particles is computed instead. - pub domain_aabb: Option>, + pub particle_aabb: Option>, /// Whether to allow multi threading within the surface reconstruction procedure pub enable_multi_threading: bool, /// Parameters for the spatial decomposition of the surface reconstruction @@ -275,7 +276,7 @@ impl Parameters { compact_support_radius: self.compact_support_radius.try_convert()?, cube_size: self.cube_size.try_convert()?, iso_surface_threshold: self.iso_surface_threshold.try_convert()?, - domain_aabb: map_option!(&self.domain_aabb, aabb => aabb.try_convert()?), + particle_aabb: map_option!(&self.particle_aabb, aabb => aabb.try_convert()?), enable_multi_threading: self.enable_multi_threading, spatial_decomposition: map_option!(&self.spatial_decomposition, sd => sd.try_convert()?), }) @@ -287,12 +288,14 @@ impl Parameters { pub struct SurfaceReconstruction { /// Background grid that was used as a basis for generating the density map for marching cubes grid: UniformGrid, - /// Octree constructed for domain decomposition + /// Octree constructed for domain decomposition (if octree-based spatial decomposition was enabled) octree: Option>, /// Point-based density map generated from the particles that was used as input to marching cubes density_map: Option>, - /// Per particle densities + /// Per particle densities (contains only data of particles inside the domain) particle_densities: Option>, + /// If an AABB was specified to restrict the reconstruction, this stores per input particle whether they were inside + particle_inside_aabb: Option>, /// Surface mesh that is the result of the surface reconstruction mesh: TriMesh3d, /// Workspace with allocated memory for subsequent surface reconstructions @@ -307,6 +310,7 @@ impl Default for SurfaceReconstruction { octree: None, density_map: None, particle_densities: None, + particle_inside_aabb: None, mesh: TriMesh3d::default(), workspace: ReconstructionWorkspace::default(), } @@ -409,13 +413,52 @@ pub fn reconstruct_surface_inplace<'a, I: Index, R: Real>( // Clear the existing mesh output_surface.mesh.clear(); + // Filter out particles + let filtered_particle_positions = if let Some(particle_aabb) = ¶meters.particle_aabb { + profile!("filtering particles"); + + use rayon::prelude::*; + let mut particle_inside = output_surface + .particle_inside_aabb + .take() + .unwrap_or_default(); + utils::reserve_total(&mut particle_inside, particle_positions.len()); + particle_positions + .par_iter() + .map(|p| particle_aabb.contains_point(p)) + .collect_into_vec(&mut particle_inside); + let particle_inside_count = particle_inside.par_iter().copied().filter(|i| *i).count(); + + // Take temporary storage for filtered particles from workspace + let mut filtered_particles = + std::mem::take(output_surface.workspace.filtered_particles_mut()); + filtered_particles.clear(); + utils::reserve_total(&mut filtered_particles, particle_inside_count); + + // Collect filtered particles + filtered_particles.extend( + particle_positions + .iter() + .zip(particle_inside.iter().copied()) + .filter(|(_, is_inside)| *is_inside) + .map(|(p, _)| p) + .cloned(), + ); + + output_surface.particle_inside_aabb = Some(particle_inside); + Cow::Owned(filtered_particles) + } else { + Cow::Borrowed(particle_positions) + }; + let particle_positions = filtered_particle_positions.as_ref(); + // Initialize grid for the reconstruction output_surface.grid = grid_for_reconstruction( particle_positions, parameters.particle_radius, parameters.compact_support_radius, parameters.cube_size, - parameters.domain_aabb.as_ref(), + parameters.particle_aabb.as_ref(), parameters.enable_multi_threading, )?; @@ -443,6 +486,12 @@ pub fn reconstruct_surface_inplace<'a, I: Index, R: Real>( )?, } + // Put back temporary storage for filtered particles for next reconstruction + if let Cow::Owned(mut filtered_particles) = filtered_particle_positions { + filtered_particles.clear(); + *output_surface.workspace.filtered_particles_mut() = filtered_particles; + } + Ok(()) } @@ -458,39 +507,38 @@ pub fn grid_for_reconstruction( particle_radius: R, compact_support_radius: R, cube_size: R, - domain_aabb: Option<&Aabb3d>, + particle_aabb: Option<&Aabb3d>, enable_multi_threading: bool, ) -> Result, ReconstructionError> { - let domain_aabb = if let Some(domain_aabb) = domain_aabb { - domain_aabb.clone() + let mut particle_aabb = if let Some(particle_aabb) = particle_aabb { + particle_aabb.clone() } else { profile!("compute minimum enclosing aabb"); - let mut domain_aabb = { + let particle_aabb = { let mut aabb = if enable_multi_threading { Aabb3d::par_from_points(particle_positions) } else { Aabb3d::from_points(particle_positions) }; + // TODO: Is this really necessary? aabb.grow_uniformly(particle_radius); aabb }; info!( "Minimal enclosing bounding box of particles was computed as: {:?}", - domain_aabb + particle_aabb ); - // Ensure that we have enough margin around the particles such that the every particle's kernel support is completely in the domain - let kernel_margin = density_map::compute_kernel_evaluation_radius::( - compact_support_radius, - cube_size, - ) - .kernel_evaluation_radius; - domain_aabb.grow_uniformly(kernel_margin); - - domain_aabb + particle_aabb }; - Ok(UniformGrid::from_aabb(&domain_aabb, cube_size)?) + // Ensure that we have enough margin around the particles such that the every particle's kernel support is completely in the domain + let kernel_margin = + density_map::compute_kernel_evaluation_radius::(compact_support_radius, cube_size) + .kernel_evaluation_radius; + particle_aabb.grow_uniformly(kernel_margin); + + Ok(UniformGrid::from_aabb(&particle_aabb, cube_size)?) } diff --git a/splashsurf_lib/src/workspace.rs b/splashsurf_lib/src/workspace.rs index 0b7f173..eac3a2e 100644 --- a/splashsurf_lib/src/workspace.rs +++ b/splashsurf_lib/src/workspace.rs @@ -12,6 +12,8 @@ use thread_local::ThreadLocal; #[derive(Default)] pub struct ReconstructionWorkspace { global_densities: Vec, + /// Temporary storage for storing a filtered set of the user provided particles + filtered_particles: Vec>, local_workspaces: ThreadLocal>>, } @@ -26,6 +28,11 @@ impl ReconstructionWorkspace { &mut self.global_densities } + /// Returns a mutable reference to the global filtered particles vector + pub(crate) fn filtered_particles_mut(&mut self) -> &mut Vec> { + &mut self.filtered_particles + } + /// Returns a reference to a thread local workspace pub(crate) fn get_local(&self) -> &RefCell> { self.local_workspaces.get_or_default() diff --git a/splashsurf_lib/tests/integration_tests/test_full.rs b/splashsurf_lib/tests/integration_tests/test_full.rs index d8c570b..4133f67 100644 --- a/splashsurf_lib/tests/integration_tests/test_full.rs +++ b/splashsurf_lib/tests/integration_tests/test_full.rs @@ -36,7 +36,7 @@ fn params_with_aabb( compact_support_radius, cube_size, iso_surface_threshold, - domain_aabb, + particle_aabb: domain_aabb, enable_multi_threading: false, spatial_decomposition: None, };