diff --git a/.github/workflows/pysplashsurf_CI.yml b/.github/workflows/pysplashsurf_CI.yml index 319b29b7..e288857a 100644 --- a/.github/workflows/pysplashsurf_CI.yml +++ b/.github/workflows/pysplashsurf_CI.yml @@ -1,12 +1,10 @@ -# This file is autogenerated by maturin v1.8.2 -# To update, run -# -# maturin generate-ci github -# name: Python bindings on: push: + branches: [main] + pull_request: + branches: [main] workflow_dispatch: release: types: [published] diff --git a/pysplashsurf/pysplashsurf/__init__.py b/pysplashsurf/pysplashsurf/__init__.py index f936c9d9..c817f0a4 100644 --- a/pysplashsurf/pysplashsurf/__init__.py +++ b/pysplashsurf/pysplashsurf/__init__.py @@ -523,11 +523,15 @@ def convert_tris_to_quads( def reconstruction_pipeline( particles, *, attributes_to_interpolate={}, particle_radius, rest_density=1000.0, smoothing_length=2.0, cube_size, - iso_surface_threshold=0.6, enable_multi_threading=True, mesh_smoothing_weights=False, sph_normals=False, + iso_surface_threshold=0.6, enable_multi_threading=True, + check_mesh_closed=False, check_mesh_manifold=False, + check_mesh_orientation=False, check_mesh_debug=False, + mesh_smoothing_weights=False, sph_normals=False, mesh_smoothing_weights_normalization=13.0, mesh_smoothing_iters=None, normals_smoothing_iters=None, - mesh_cleanup=False, decimate_barnacles=False, keep_vertices=False, - compute_normals=False, output_raw_normals=False, output_mesh_smoothing_weights=False, mesh_aabb_clamp_vertices=False, - subdomain_grid=True, subdomain_num_cubes_per_dim=64, aabb_min=None, aabb_max=None, mesh_aabb_min=None, mesh_aabb_max=None + mesh_cleanup=False, mesh_cleanup_snap_dist=None, decimate_barnacles=False, keep_vertices=False, + compute_normals=False, output_raw_normals=False, output_raw_mesh=False, output_mesh_smoothing_weights=False, mesh_aabb_clamp_vertices=False, + subdomain_grid=True, subdomain_num_cubes_per_dim=64, aabb_min=None, aabb_max=None, mesh_aabb_min=None, mesh_aabb_max=None, + generate_quads=False, quad_max_edge_diag_ratio=1.75, quad_max_normal_angle=10.0, quad_max_interior_angle=135.0 ): """Surface reconstruction based on particle positions with subsequent post-processing @@ -558,6 +562,18 @@ def reconstruction_pipeline( enable_multi_threading: bool Multi-threading flag + + check_mesh_closed: bool + Enable checking the final mesh for holes + + check_mesh_manifold: bool + Enable checking the final mesh for non-manifold edges and vertices + + check_mesh_orientation: bool + Enable checking the final mesh for inverted triangles (compares angle between vertex normals and adjacent face normals) + + check_mesh_debug: bool + Enable additional debug output for the check-mesh operations (has no effect if no other check-mesh option is enabled) sph_normals: bool Flag to compute normals using SPH interpolation instead of geometry-based normals. @@ -578,6 +594,9 @@ def reconstruction_pipeline( mesh_cleanup: bool Flag to perform mesh cleanup\n This implements the method from “Compact isocontours from sampled data” (Moore, Warren; 1992) + + mesh_cleanup_snap_dist: float + If MC mesh cleanup is enabled, vertex snapping can be limited to this distance relative to the MC edge length (should be in range of [0.0,0.5]) decimate_barnacles: bool Flag to perform barnacle decimation\n @@ -590,11 +609,14 @@ def reconstruction_pipeline( Flag to compute normals\n If set to True, the normals will be computed and stored in the mesh. + output_mesh_smoothing_weights: bool + Flag to store the mesh smoothing weights if smoothing weights are computed. + output_raw_normals: bool Flag to output the raw normals in addition to smoothed normals if smoothing of normals is enabled - output_mesh_smoothing_weights: bool - Flag to store the mesh smoothing weights if smoothing weights are computed. + output_raw_mesh: bool + When true, also return the SurfaceReconstruction object with no post-processing applied mesh_aabb_clamp_vertices: bool Flag to clamp the vertices of the mesh to the AABB @@ -616,29 +638,58 @@ def reconstruction_pipeline( mesh_aabb_max: np.ndarray Largest corner of the axis-aligned bounding box for the mesh + + generate_quads: bool + Enable trying to convert triangles to quads if they meet quality criteria + + quad_max_edge_diag_ratio: float + Maximum allowed ratio of quad edge lengths to its diagonals to merge two triangles to a quad (inverse is used for minimum) + + quad_max_normal_angle: float + Maximum allowed angle (in degrees) between triangle normals to merge them to a quad + + quad_max_interior_angle: float + Maximum allowed vertex interior angle (in degrees) inside a quad to merge two triangles to a quad Returns ------- - tuple[TriMeshWithDataF32 | TriMeshWithDataF64, SurfaceReconstructionF32 | SurfaceReconstructionF64] + tuple[TriMeshWithDataF32 | TriMeshWithDataF64 | MixedTriQuadMeshWithDataF32 | MixedTriQuadMeshWithDataF64, Optional[SurfaceReconstructionF32] | Optional[SurfaceReconstructionF64]] Mesh with data object and SurfaceReconstruction object containing the reconstructed mesh and used grid """ if particles.dtype == 'float32': - return reconstruction_pipeline_f32(particles, attributes_to_interpolate=attributes_to_interpolate, particle_radius=particle_radius, rest_density=rest_density, - smoothing_length=smoothing_length, cube_size=cube_size, iso_surface_threshold=iso_surface_threshold, - aabb_min=aabb_min, aabb_max=aabb_max, enable_multi_threading=enable_multi_threading, - use_custom_grid_decomposition=subdomain_grid, subdomain_num_cubes_per_dim=subdomain_num_cubes_per_dim, - global_neighborhood_list=False, mesh_cleanup=mesh_cleanup, decimate_barnacles=decimate_barnacles, - keep_vertices=keep_vertices, compute_normals=compute_normals, sph_normals=sph_normals, normals_smoothing_iters=normals_smoothing_iters, - mesh_smoothing_iters=mesh_smoothing_iters, mesh_smoothing_weights=mesh_smoothing_weights, mesh_smoothing_weights_normalization=mesh_smoothing_weights_normalization, - output_mesh_smoothing_weights=output_mesh_smoothing_weights, output_raw_normals=output_raw_normals, mesh_aabb_min=mesh_aabb_min, mesh_aabb_max=mesh_aabb_max, mesh_aabb_clamp_vertices=mesh_aabb_clamp_vertices) + tri_mesh, tri_quad_mesh, reconstruction = reconstruction_pipeline_f32(particles, attributes_to_interpolate=attributes_to_interpolate, particle_radius=particle_radius, rest_density=rest_density, + smoothing_length=smoothing_length, cube_size=cube_size, iso_surface_threshold=iso_surface_threshold, + aabb_min=aabb_min, aabb_max=aabb_max, enable_multi_threading=enable_multi_threading, + use_custom_grid_decomposition=subdomain_grid, subdomain_num_cubes_per_dim=subdomain_num_cubes_per_dim, + check_mesh_closed=check_mesh_closed, check_mesh_manifold=check_mesh_manifold, check_mesh_orientation=check_mesh_orientation, check_mesh_debug=check_mesh_debug, + mesh_cleanup=mesh_cleanup, max_rel_snap_dist=mesh_cleanup_snap_dist, decimate_barnacles=decimate_barnacles, + keep_vertices=keep_vertices, compute_normals=compute_normals, sph_normals=sph_normals, normals_smoothing_iters=normals_smoothing_iters, + mesh_smoothing_iters=mesh_smoothing_iters, mesh_smoothing_weights=mesh_smoothing_weights, mesh_smoothing_weights_normalization=mesh_smoothing_weights_normalization, + output_mesh_smoothing_weights=output_mesh_smoothing_weights, output_raw_normals=output_raw_normals, output_raw_mesh=output_raw_mesh, + mesh_aabb_min=mesh_aabb_min, mesh_aabb_max=mesh_aabb_max, mesh_aabb_clamp_vertices=mesh_aabb_clamp_vertices, + generate_quads=generate_quads, quad_max_edge_diag_ratio=quad_max_edge_diag_ratio, quad_max_normal_angle=quad_max_normal_angle, quad_max_interior_angle=quad_max_interior_angle) + + if tri_mesh == None: + return (tri_quad_mesh, reconstruction) + else: + return (tri_mesh, reconstruction) + elif particles.dtype == 'float64': - return reconstruction_pipeline_f64(particles, attributes_to_interpolate=attributes_to_interpolate, particle_radius=particle_radius, rest_density=rest_density, - smoothing_length=smoothing_length, cube_size=cube_size, iso_surface_threshold=iso_surface_threshold, - aabb_min=aabb_min, aabb_max=aabb_max, enable_multi_threading=enable_multi_threading, - use_custom_grid_decomposition=subdomain_grid, subdomain_num_cubes_per_dim=subdomain_num_cubes_per_dim, - global_neighborhood_list=False, mesh_cleanup=mesh_cleanup, decimate_barnacles=decimate_barnacles, - keep_vertices=keep_vertices, compute_normals=compute_normals, sph_normals=sph_normals, normals_smoothing_iters=normals_smoothing_iters, - mesh_smoothing_iters=mesh_smoothing_iters, mesh_smoothing_weights=mesh_smoothing_weights, mesh_smoothing_weights_normalization=mesh_smoothing_weights_normalization, - output_mesh_smoothing_weights=output_mesh_smoothing_weights, output_raw_normals=output_raw_normals, mesh_aabb_min=mesh_aabb_min, mesh_aabb_max=mesh_aabb_max, mesh_aabb_clamp_vertices=mesh_aabb_clamp_vertices) + tri_mesh, tri_quad_mesh, reconstruction = reconstruction_pipeline_f64(particles, attributes_to_interpolate=attributes_to_interpolate, particle_radius=particle_radius, rest_density=rest_density, + smoothing_length=smoothing_length, cube_size=cube_size, iso_surface_threshold=iso_surface_threshold, + aabb_min=aabb_min, aabb_max=aabb_max, enable_multi_threading=enable_multi_threading, + use_custom_grid_decomposition=subdomain_grid, subdomain_num_cubes_per_dim=subdomain_num_cubes_per_dim, + check_mesh_closed=check_mesh_closed, check_mesh_manifold=check_mesh_manifold, check_mesh_orientation=check_mesh_orientation, check_mesh_debug=check_mesh_debug, + mesh_cleanup=mesh_cleanup, max_rel_snap_dist=mesh_cleanup_snap_dist, decimate_barnacles=decimate_barnacles, + keep_vertices=keep_vertices, compute_normals=compute_normals, sph_normals=sph_normals, normals_smoothing_iters=normals_smoothing_iters, + mesh_smoothing_iters=mesh_smoothing_iters, mesh_smoothing_weights=mesh_smoothing_weights, mesh_smoothing_weights_normalization=mesh_smoothing_weights_normalization, + output_mesh_smoothing_weights=output_mesh_smoothing_weights, output_raw_normals=output_raw_normals, output_raw_mesh=output_raw_mesh, + mesh_aabb_min=mesh_aabb_min, mesh_aabb_max=mesh_aabb_max, mesh_aabb_clamp_vertices=mesh_aabb_clamp_vertices, + generate_quads=generate_quads, quad_max_edge_diag_ratio=quad_max_edge_diag_ratio, quad_max_normal_angle=quad_max_normal_angle, quad_max_interior_angle=quad_max_interior_angle) + + if tri_mesh == None: + return (tri_quad_mesh, reconstruction) + else: + return (tri_mesh, reconstruction) else: raise ValueError("Invalid data type (only float32 and float64 are supported, consider explicitly specifying the dtype for particles)") diff --git a/pysplashsurf/pysplashsurf/docs/source/classes.rst b/pysplashsurf/pysplashsurf/docs/source/classes.rst index afb33aa8..8060dcfb 100644 --- a/pysplashsurf/pysplashsurf/docs/source/classes.rst +++ b/pysplashsurf/pysplashsurf/docs/source/classes.rst @@ -1,24 +1,42 @@ Classes ======= -Additionally, there exists a F64 version for every class which is otherwise identical to the F32 version. +Additionally to the classes on this page, there exists a F64 version for every class which is otherwise identical to the F32 version. + +For more information on the classes, refer to the `Rust documentation `_ of splashsurf_lib. .. currentmodule:: pysplashsurf .. autoclass:: Aabb3dF32 + See `Aabb3d `_ for more information. + .. autoclass:: MixedTriQuadMesh3dF32 + See `MixedTriQuadMesh3d `_ for more information. + .. autoclass:: MixedTriQuadMeshWithDataF32 :exclude-members: push_point_attribute_scalar_u64, push_point_attribute_scalar_real, push_point_attribute_vector_real, push_cell_attribute_scalar_real, push_cell_attribute_scalar_u64, push_cell_attribute_vector_real + See `MeshWithData `_ for more information. + .. autoclass:: SphInterpolatorF32 + See `SphInterpolator `_ for more information. + .. autoclass:: SurfaceReconstructionF32 + See `SurfaceReconstruction `_ for more information. + .. autoclass:: TriMesh3dF32 + See `TriMesh3d `_ for more information. + .. autoclass:: TriMeshWithDataF32 :exclude-members: push_point_attribute_scalar_u64, push_point_attribute_scalar_real, push_point_attribute_vector_real, push_cell_attribute_scalar_real, push_cell_attribute_scalar_u64, push_cell_attribute_vector_real -.. autoclass:: UniformGridF32 \ No newline at end of file + See `MeshWithData `_ for more information. + +.. autoclass:: UniformGridF32 + + See `UniformGrid `_ for more information. \ No newline at end of file diff --git a/pysplashsurf/src/pipeline.rs b/pysplashsurf/src/pipeline.rs index 5b451738..f10a7ba0 100644 --- a/pysplashsurf/src/pipeline.rs +++ b/pysplashsurf/src/pipeline.rs @@ -1,27 +1,23 @@ -use anyhow::anyhow; -use log::info; +use crate::{ + mesh::{ + MixedTriQuadMeshWithDataF32, MixedTriQuadMeshWithDataF64, TriMeshWithDataF32, + TriMeshWithDataF64, + }, + reconstruction::{SurfaceReconstructionF32, SurfaceReconstructionF64}, +}; use numpy::{Element, PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2}; use pyo3::{ prelude::*, types::{PyDict, PyString}, }; -use rayon::iter::{IndexedParallelIterator, IntoParallelRefIterator, ParallelIterator}; use splashsurf_lib::{ - Aabb3d, Index, Real, SurfaceReconstruction, - mesh::{AttributeData, Mesh3d, MeshAttribute, MeshWithData, TriMesh3d}, - nalgebra::{Unit, Vector3}, - profile, - sph_interpolation::SphInterpolator, -}; -use std::borrow::Cow; - -use crate::{ - mesh::{TriMeshWithDataF32, TriMeshWithDataF64}, - reconstruction::{SurfaceReconstructionF32, SurfaceReconstructionF64, reconstruct_surface_py}, + Aabb3d, GridDecompositionParameters, Index, Real, SpatialDecomposition, + mesh::{AttributeData, MeshAttribute}, + nalgebra::Vector3, }; fn reconstruction_pipeline_generic( - particles: &[Vector3], + particle_positions: &[Vector3], attributes: Vec>, particle_radius: R, rest_density: R, @@ -33,13 +29,12 @@ fn reconstruction_pipeline_generic( enable_multi_threading: bool, use_custom_grid_decomposition: bool, subdomain_num_cubes_per_dim: u32, - global_neighborhood_list: bool, - // check_mesh_closed: bool, - // check_mesh_manifold: bool, - // check_mesh_orientation: bool, - // check_mesh_debug: bool, + check_mesh_closed: bool, + check_mesh_manifold: bool, + check_mesh_orientation: bool, + check_mesh_debug: bool, mesh_cleanup: bool, - max_rel_snap_dist: Option, + max_rel_snap_dist: Option, decimate_barnacles: bool, keep_vertices: bool, compute_normals: bool, @@ -48,378 +43,93 @@ fn reconstruction_pipeline_generic( mesh_smoothing_iters: Option, mesh_smoothing_weights: bool, mesh_smoothing_weights_normalization: f64, - // generate_quads: bool, - // quad_max_edge_diag_ratio: f64, - // quad_max_normal_angle: f64, - // quad_max_interior_angle: f64, + generate_quads: bool, + quad_max_edge_diag_ratio: f64, + quad_max_normal_angle: f64, + quad_max_interior_angle: f64, output_mesh_smoothing_weights: bool, output_raw_normals: bool, - mesh_aabb_min: Option<[R; 3]>, - mesh_aabb_max: Option<[R; 3]>, + output_raw_mesh: bool, + mesh_aabb_min: Option<[f64; 3]>, + mesh_aabb_max: Option<[f64; 3]>, mesh_aabb_clamp_vertices: bool, -) -> Result<(MeshWithData>, SurfaceReconstruction), anyhow::Error> { - profile!("surface reconstruction"); +) -> Result, anyhow::Error> { + let aabb = if let (Some(aabb_min), Some(aabb_max)) = (aabb_min, aabb_max) { + // Convert the min and max arrays to Vector3 + Some(Aabb3d::new( + Vector3::from(aabb_min), + Vector3::from(aabb_max), + )) + } else { + None + }; - let compact_support_radius = R::from_f64(2.0).unwrap() * smoothing_length * particle_radius; + let spatial_decomposition = use_custom_grid_decomposition.then(|| { + let mut grid_params = GridDecompositionParameters::default(); + grid_params.subdomain_num_cubes_per_dim = subdomain_num_cubes_per_dim; + SpatialDecomposition::UniformGrid(grid_params) + }); - // Perform the surface reconstruction - let reconstruction = reconstruct_surface_py::( - particles, + let params = splashsurf_lib::Parameters { particle_radius, rest_density, - smoothing_length, - cube_size, + compact_support_radius: R::from_f64(2.0).unwrap() * smoothing_length * particle_radius, + cube_size: cube_size * particle_radius, iso_surface_threshold, + particle_aabb: aabb, enable_multi_threading, - global_neighborhood_list, - use_custom_grid_decomposition, - subdomain_num_cubes_per_dim, - aabb_min, - aabb_max, - ); - - // let grid = reconstruction.grid(); - let mut mesh_with_data: MeshWithData> = - MeshWithData::new(reconstruction.mesh().clone()); - - // Perform post-processing - { - profile!("postprocessing"); - let mut vertex_connectivity = None; - - if mesh_cleanup { - info!("Post-processing: Performing mesh cleanup"); - let tris_before = mesh_with_data.mesh.triangles.len(); - let verts_before = mesh_with_data.mesh.vertices.len(); - vertex_connectivity = Some(splashsurf_lib::postprocessing::marching_cubes_cleanup( - &mut mesh_with_data.mesh, - reconstruction.grid(), - max_rel_snap_dist, - 5, - keep_vertices, - )); - let tris_after = mesh_with_data.mesh.triangles.len(); - let verts_after = mesh_with_data.mesh.vertices.len(); - info!( - "Post-processing: Cleanup reduced number of vertices to {:.2}% and number of triangles to {:.2}% of original mesh.", - (verts_after as f64 / verts_before as f64) * 100.0, - (tris_after as f64 / tris_before as f64) * 100.0 - ) - } - - // Decimate mesh if requested - if decimate_barnacles { - info!("Post-processing: Performing decimation"); - vertex_connectivity = Some(splashsurf_lib::postprocessing::decimation( - &mut mesh_with_data.mesh, - keep_vertices, - )); - } - - // Initialize SPH interpolator if required later - let interpolator_required = mesh_smoothing_weights || sph_normals; - - let interpolator = if interpolator_required { - profile!("initialize interpolator"); - info!("Post-processing: Initializing interpolator..."); - - info!( - "Constructing global acceleration structure for SPH interpolation to {} vertices...", - mesh_with_data.vertices().len() - ); - - let particle_rest_density = rest_density; - let particle_rest_volume = - R::from_float(4.0) * R::frac_pi_3() * particle_radius.powi(3); - let particle_rest_mass = particle_rest_volume * particle_rest_density; - - let particle_densities = reconstruction - .particle_densities() - .ok_or_else(|| anyhow::anyhow!("Particle densities were not returned by surface reconstruction but are required for SPH normal computation"))? - .as_slice(); - assert_eq!( - particles.len(), - particle_densities.len(), - "There has to be one density value per particle" - ); + spatial_decomposition, + global_neighborhood_list: mesh_smoothing_weights, + }; - Some(SphInterpolator::new( - &particles, - particle_densities, - particle_rest_mass, - compact_support_radius, + let mesh_aabb = + if let (Some(mesh_aabb_min), Some(mesh_aabb_max)) = (mesh_aabb_min, mesh_aabb_max) { + // Convert the min and max arrays to Vector3 + Some(Aabb3d::new( + Vector3::from(mesh_aabb_min), + Vector3::from(mesh_aabb_max), )) } else { None }; - // Compute mesh vertex-vertex connectivity map if required later - let vertex_connectivity_required = - normals_smoothing_iters.is_some() || mesh_smoothing_iters.is_some(); - if vertex_connectivity.is_none() && vertex_connectivity_required { - vertex_connectivity = Some(mesh_with_data.mesh.vertex_vertex_connectivity()); - } - - // Compute smoothing weights if requested - let smoothing_weights = if mesh_smoothing_weights { - profile!("compute smoothing weights"); - info!("Post-processing: Computing smoothing weights..."); - - // TODO: Switch between parallel/single threaded - // TODO: Re-use data from reconstruction? - - // Global neighborhood search - let nl = reconstruction - .particle_neighbors() - .map(Cow::Borrowed) - .unwrap_or_else(|| - { - let search_radius = compact_support_radius; - - let mut domain = Aabb3d::from_points(particles); - domain.grow_uniformly(search_radius); - - let mut nl = Vec::new(); - splashsurf_lib::neighborhood_search::neighborhood_search_spatial_hashing_parallel::( - &domain, - particles, - search_radius, - &mut nl, - ); - assert_eq!(nl.len(), particles.len()); - Cow::Owned(nl) - } - ); - - // Compute weighted neighbor count - let squared_r = compact_support_radius * compact_support_radius; - let weighted_ncounts = nl - .par_iter() - .enumerate() - .map(|(i, nl)| { - nl.iter() - .copied() - .map(|j| { - let dist = (particles[i] - particles[j]).norm_squared(); - - R::one() - (dist / squared_r).clamp(R::zero(), R::one()) - }) - .fold(R::zero(), R::add) - }) - .collect::>(); - - let vertex_weighted_num_neighbors = { - profile!("interpolate weighted neighbor counts"); - interpolator - .as_ref() - .expect("interpolator is required") - .interpolate_scalar_quantity( - weighted_ncounts.as_slice(), - mesh_with_data.vertices(), - true, - ) - }; - - let smoothing_weights = { - let offset = R::zero(); - let normalization = R::from_f64(mesh_smoothing_weights_normalization).expect( - "smoothing weight normalization value cannot be represented as Real type", - ) - offset; - - // Normalize number of neighbors - let smoothing_weights = vertex_weighted_num_neighbors - .par_iter() - .copied() - .map(|n| (n - offset).max(R::zero())) - .map(|n| (n / normalization).min(R::one())) - // Smooth-Step function - .map(|x| x.powi(5).times(6) - x.powi(4).times(15) + x.powi(3).times(10)) - .collect::>(); - - if output_mesh_smoothing_weights { - // Raw distance-weighted number of neighbors value per vertex (can be used to determine normalization value) - mesh_with_data.point_attributes.push(MeshAttribute::new( - "wnn".to_string(), - AttributeData::ScalarReal(vertex_weighted_num_neighbors), - )); - // Final smoothing weights per vertex - mesh_with_data.point_attributes.push(MeshAttribute::new( - "sw".to_string(), - AttributeData::ScalarReal(smoothing_weights.clone()), - )); - } - - smoothing_weights - }; - - Some(smoothing_weights) - } else { - None - }; - - // Perform smoothing if requested - if let Some(mesh_smoothing_iters) = mesh_smoothing_iters { - profile!("mesh smoothing"); - info!("Post-processing: Smoothing mesh..."); - - // TODO: Switch between parallel/single threaded - - let smoothing_weights = smoothing_weights - .unwrap_or_else(|| vec![R::one(); mesh_with_data.vertices().len()]); - - splashsurf_lib::postprocessing::par_laplacian_smoothing_inplace( - &mut mesh_with_data.mesh, - vertex_connectivity - .as_ref() - .expect("vertex connectivity is required"), - mesh_smoothing_iters, - R::one(), - &smoothing_weights, - ); - } - - // Add normals to mesh if requested - if compute_normals { - profile!("compute normals"); - info!("Post-processing: Computing surface normals..."); - - // Compute normals - let normals = if sph_normals { - info!("Using SPH interpolation to compute surface normals"); - - let sph_normals = interpolator - .as_ref() - .expect("interpolator is required") - .interpolate_normals(mesh_with_data.vertices()); - bytemuck::allocation::cast_vec::>, Vector3>(sph_normals) - } else { - info!("Using area weighted triangle normals for surface normals"); - profile!("mesh.par_vertex_normals"); - let tri_normals = mesh_with_data.mesh.par_vertex_normals(); - - // Convert unit vectors to plain vectors - bytemuck::allocation::cast_vec::>, Vector3>(tri_normals) - }; - - // Smooth normals - if let Some(smoothing_iters) = normals_smoothing_iters { - info!("Post-processing: Smoothing normals..."); - - let mut smoothed_normals = normals.clone(); - splashsurf_lib::postprocessing::par_laplacian_smoothing_normals_inplace( - &mut smoothed_normals, - vertex_connectivity - .as_ref() - .expect("vertex connectivity is required"), - smoothing_iters, - ); - - mesh_with_data.point_attributes.push(MeshAttribute::new( - "normals".to_string(), - AttributeData::Vector3Real(smoothed_normals), - )); - if output_raw_normals { - mesh_with_data.point_attributes.push(MeshAttribute::new( - "raw_normals".to_string(), - AttributeData::Vector3Real(normals), - )); - } - } else { - mesh_with_data.point_attributes.push(MeshAttribute::new( - "normals".to_string(), - AttributeData::Vector3Real(normals), - )); - } - } - - // Interpolate attributes if requested - if !attributes.is_empty() { - profile!("interpolate attributes"); - info!("Post-processing: Interpolating attributes..."); - let interpolator = interpolator.as_ref().expect("interpolator is required"); - - for attribute in attributes.into_iter() { - info!("Interpolating attribute \"{}\"...", attribute.name); - - match attribute.data { - AttributeData::ScalarReal(values) => { - let interpolated_values = interpolator.interpolate_scalar_quantity( - values.as_slice(), - mesh_with_data.vertices(), - true, - ); - mesh_with_data.point_attributes.push(MeshAttribute::new( - attribute.name, - AttributeData::ScalarReal(interpolated_values), - )); - } - AttributeData::Vector3Real(values) => { - let interpolated_values = interpolator.interpolate_vector_quantity( - values.as_slice(), - mesh_with_data.vertices(), - true, - ); - mesh_with_data.point_attributes.push(MeshAttribute::new( - attribute.name, - AttributeData::Vector3Real(interpolated_values), - )); - } - _ => unimplemented!("Interpolation of this attribute type not implemented"), - } - } - } - } - - // Remove and clamp cells outside of AABB - let mesh_aabb = if aabb_min != None && aabb_max != None { - Some(Aabb3d::new( - Vector3::from(mesh_aabb_min.unwrap()), - Vector3::from(mesh_aabb_max.unwrap()), - )) - } else { - None - }; - - let mesh_with_data = if let Some(mesh_aabb) = &mesh_aabb { - profile!("clamp mesh to aabb"); - info!("Post-processing: Clamping mesh to AABB..."); - - mesh_with_data.par_clamp_with_aabb( - &mesh_aabb - .try_convert() - .ok_or_else(|| anyhow!("Failed to convert mesh AABB"))?, - mesh_aabb_clamp_vertices, - keep_vertices, - ) - } else { - mesh_with_data + let postprocessing_args = splashsurf::reconstruct::ReconstructionPostprocessingParameters { + check_mesh_closed, + check_mesh_manifold, + check_mesh_orientation, + check_mesh_debug, + mesh_cleanup, + mesh_cleanup_snap_dist: max_rel_snap_dist, + decimate_barnacles, + keep_vertices, + compute_normals, + sph_normals, + normals_smoothing_iters, + interpolate_attributes: Vec::new(), + mesh_smoothing_iters, + mesh_smoothing_weights, + mesh_smoothing_weights_normalization, + generate_quads, + quad_max_edge_diag_ratio, + quad_max_normal_angle, + quad_max_interior_angle, + output_mesh_smoothing_weights, + output_raw_normals, + output_raw_mesh, + mesh_aabb, + mesh_aabb_clamp_vertices, }; - // Convert triangles to quads - // let (tri_mesh, tri_quad_mesh) = if generate_quads { - // info!("Post-processing: Convert triangles to quads..."); - // let non_squareness_limit = R::from_f64(quad_max_edge_diag_ratio).unwrap(); - // let normal_angle_limit_rad = - // R::from_f64(quad_max_normal_angle.to_radians()).unwrap(); - // let max_interior_angle = - // R::from_f64(quad_max_interior_angle.to_radians()).unwrap(); - - // let tri_quad_mesh = splashsurf_lib::postprocessing::convert_tris_to_quads( - // &mesh_with_data.mesh, - // non_squareness_limit, - // normal_angle_limit_rad, - // max_interior_angle, - // ); - - // (None, Some(mesh_with_data.with_mesh(tri_quad_mesh))) - // } else { - // (Some(mesh_with_data), None) - // }; - Ok((mesh_with_data, reconstruction)) + splashsurf::reconstruct::reconstruction_pipeline( + particle_positions, + attributes, + ¶ms, + &postprocessing_args, + ) } -fn attrs_conversion<'py, R: Real + Element>( - attributes_to_interpolate: Bound<'py, PyDict>, +fn attrs_conversion( + attributes_to_interpolate: Bound, ) -> Vec> { let mut attrs: Vec> = Vec::new(); for (key, value) in attributes_to_interpolate.iter() { @@ -468,11 +178,13 @@ fn attrs_conversion<'py, R: Real + Element>( #[pyo3(signature = (particles, *, attributes_to_interpolate, particle_radius, rest_density, smoothing_length, cube_size, iso_surface_threshold, aabb_min = None, aabb_max = None, enable_multi_threading = false, - use_custom_grid_decomposition = false, subdomain_num_cubes_per_dim = 64, global_neighborhood_list = false, + use_custom_grid_decomposition = false, subdomain_num_cubes_per_dim = 64, + check_mesh_closed = false, check_mesh_manifold = false, check_mesh_orientation = false, check_mesh_debug = false, mesh_cleanup, max_rel_snap_dist = None, decimate_barnacles, keep_vertices, compute_normals, sph_normals, - normals_smoothing_iters, mesh_smoothing_iters, mesh_smoothing_weights, - mesh_smoothing_weights_normalization, output_mesh_smoothing_weights, - output_raw_normals, mesh_aabb_min, mesh_aabb_max, mesh_aabb_clamp_vertices + normals_smoothing_iters, mesh_smoothing_iters, mesh_smoothing_weights, mesh_smoothing_weights_normalization, + generate_quads = false, quad_max_edge_diag_ratio = 1.75, quad_max_normal_angle = 10.0, quad_max_interior_angle = 135.0, + output_mesh_smoothing_weights, output_raw_normals, output_raw_mesh=false, + mesh_aabb_min, mesh_aabb_max, mesh_aabb_clamp_vertices ))] pub fn reconstruction_pipeline_py_f32<'py>( particles: &Bound<'py, PyArray2>, @@ -487,9 +199,12 @@ pub fn reconstruction_pipeline_py_f32<'py>( enable_multi_threading: bool, use_custom_grid_decomposition: bool, subdomain_num_cubes_per_dim: u32, - global_neighborhood_list: bool, + check_mesh_closed: bool, + check_mesh_manifold: bool, + check_mesh_orientation: bool, + check_mesh_debug: bool, mesh_cleanup: bool, - max_rel_snap_dist: Option, + max_rel_snap_dist: Option, decimate_barnacles: bool, keep_vertices: bool, compute_normals: bool, @@ -498,20 +213,33 @@ pub fn reconstruction_pipeline_py_f32<'py>( mesh_smoothing_iters: Option, mesh_smoothing_weights: bool, mesh_smoothing_weights_normalization: f64, + generate_quads: bool, + quad_max_edge_diag_ratio: f64, + quad_max_normal_angle: f64, + quad_max_interior_angle: f64, output_mesh_smoothing_weights: bool, output_raw_normals: bool, - mesh_aabb_min: Option<[f32; 3]>, - mesh_aabb_max: Option<[f32; 3]>, + output_raw_mesh: bool, + mesh_aabb_min: Option<[f64; 3]>, + mesh_aabb_max: Option<[f64; 3]>, mesh_aabb_clamp_vertices: bool, -) -> (TriMeshWithDataF32, SurfaceReconstructionF32) { - let particles: PyReadonlyArray2 = particles.extract().unwrap(); - - let particle_positions = particles.as_slice().unwrap(); +) -> PyResult<( + Option, + Option, + Option, +)> { + let particles: PyReadonlyArray2 = particles.extract()?; + + let particle_positions = particles.as_slice()?; let particle_positions: &[Vector3] = bytemuck::cast_slice(particle_positions); let attrs = attrs_conversion(attributes_to_interpolate); - let (mesh, reconstruction) = reconstruction_pipeline_generic::( + let splashsurf::reconstruct::ReconstructionResult { + tri_mesh, + tri_quad_mesh, + raw_reconstruction: reconstruction, + } = reconstruction_pipeline_generic::( particle_positions, attrs, particle_radius, @@ -524,7 +252,10 @@ pub fn reconstruction_pipeline_py_f32<'py>( enable_multi_threading, use_custom_grid_decomposition, subdomain_num_cubes_per_dim, - global_neighborhood_list, + check_mesh_closed, + check_mesh_manifold, + check_mesh_orientation, + check_mesh_debug, mesh_cleanup, max_rel_snap_dist, decimate_barnacles, @@ -535,18 +266,24 @@ pub fn reconstruction_pipeline_py_f32<'py>( mesh_smoothing_iters, mesh_smoothing_weights, mesh_smoothing_weights_normalization, + generate_quads, + quad_max_edge_diag_ratio, + quad_max_normal_angle, + quad_max_interior_angle, output_mesh_smoothing_weights, output_raw_normals, + output_raw_mesh, mesh_aabb_min, mesh_aabb_max, mesh_aabb_clamp_vertices, ) .unwrap(); - ( - TriMeshWithDataF32::new(mesh), - SurfaceReconstructionF32::new(reconstruction), - ) + Ok(( + tri_mesh.map(TriMeshWithDataF32::new), + tri_quad_mesh.map(MixedTriQuadMeshWithDataF32::new), + reconstruction.map(SurfaceReconstructionF32::new), + )) } #[pyfunction] @@ -554,11 +291,13 @@ pub fn reconstruction_pipeline_py_f32<'py>( #[pyo3(signature = (particles, *, attributes_to_interpolate, particle_radius, rest_density, smoothing_length, cube_size, iso_surface_threshold, aabb_min = None, aabb_max = None, enable_multi_threading = false, - use_custom_grid_decomposition = false, subdomain_num_cubes_per_dim = 64, global_neighborhood_list = false, + use_custom_grid_decomposition = false, subdomain_num_cubes_per_dim = 64, + check_mesh_closed = false, check_mesh_manifold = false, check_mesh_orientation = false, check_mesh_debug = false, mesh_cleanup, max_rel_snap_dist = None, decimate_barnacles, keep_vertices, compute_normals, sph_normals, - normals_smoothing_iters, mesh_smoothing_iters, mesh_smoothing_weights, - mesh_smoothing_weights_normalization, output_mesh_smoothing_weights, - output_raw_normals, mesh_aabb_min, mesh_aabb_max, mesh_aabb_clamp_vertices + normals_smoothing_iters, mesh_smoothing_iters, mesh_smoothing_weights, mesh_smoothing_weights_normalization, + generate_quads = false, quad_max_edge_diag_ratio = 1.75, quad_max_normal_angle = 10.0, quad_max_interior_angle = 135.0, + output_mesh_smoothing_weights, output_raw_normals, output_raw_mesh=false, + mesh_aabb_min, mesh_aabb_max, mesh_aabb_clamp_vertices ))] pub fn reconstruction_pipeline_py_f64<'py>( particles: &Bound<'py, PyArray2>, @@ -573,7 +312,10 @@ pub fn reconstruction_pipeline_py_f64<'py>( enable_multi_threading: bool, use_custom_grid_decomposition: bool, subdomain_num_cubes_per_dim: u32, - global_neighborhood_list: bool, + check_mesh_closed: bool, + check_mesh_manifold: bool, + check_mesh_orientation: bool, + check_mesh_debug: bool, mesh_cleanup: bool, max_rel_snap_dist: Option, decimate_barnacles: bool, @@ -584,20 +326,33 @@ pub fn reconstruction_pipeline_py_f64<'py>( mesh_smoothing_iters: Option, mesh_smoothing_weights: bool, mesh_smoothing_weights_normalization: f64, + generate_quads: bool, + quad_max_edge_diag_ratio: f64, + quad_max_normal_angle: f64, + quad_max_interior_angle: f64, output_mesh_smoothing_weights: bool, output_raw_normals: bool, + output_raw_mesh: bool, mesh_aabb_min: Option<[f64; 3]>, mesh_aabb_max: Option<[f64; 3]>, mesh_aabb_clamp_vertices: bool, -) -> (TriMeshWithDataF64, SurfaceReconstructionF64) { - let particles: PyReadonlyArray2 = particles.extract().unwrap(); - - let particle_positions = particles.as_slice().unwrap(); +) -> PyResult<( + Option, + Option, + Option, +)> { + let particles: PyReadonlyArray2 = particles.extract()?; + + let particle_positions = particles.as_slice()?; let particle_positions: &[Vector3] = bytemuck::cast_slice(particle_positions); let attrs = attrs_conversion(attributes_to_interpolate); - let (mesh, reconstruction) = reconstruction_pipeline_generic::( + let splashsurf::reconstruct::ReconstructionResult { + tri_mesh, + tri_quad_mesh, + raw_reconstruction: reconstruction, + } = reconstruction_pipeline_generic::( particle_positions, attrs, particle_radius, @@ -610,7 +365,10 @@ pub fn reconstruction_pipeline_py_f64<'py>( enable_multi_threading, use_custom_grid_decomposition, subdomain_num_cubes_per_dim, - global_neighborhood_list, + check_mesh_closed, + check_mesh_manifold, + check_mesh_orientation, + check_mesh_debug, mesh_cleanup, max_rel_snap_dist, decimate_barnacles, @@ -621,16 +379,22 @@ pub fn reconstruction_pipeline_py_f64<'py>( mesh_smoothing_iters, mesh_smoothing_weights, mesh_smoothing_weights_normalization, + generate_quads, + quad_max_edge_diag_ratio, + quad_max_normal_angle, + quad_max_interior_angle, output_mesh_smoothing_weights, output_raw_normals, + output_raw_mesh, mesh_aabb_min, mesh_aabb_max, mesh_aabb_clamp_vertices, ) .unwrap(); - ( - TriMeshWithDataF64::new(mesh), - SurfaceReconstructionF64::new(reconstruction), - ) + Ok(( + tri_mesh.map(TriMeshWithDataF64::new), + tri_quad_mesh.map(MixedTriQuadMeshWithDataF64::new), + reconstruction.map(SurfaceReconstructionF64::new), + )) } diff --git a/pysplashsurf/tests/test_calling.py b/pysplashsurf/tests/test_calling.py index eee399d8..07442f7c 100644 --- a/pysplashsurf/tests/test_calling.py +++ b/pysplashsurf/tests/test_calling.py @@ -84,8 +84,8 @@ def reconstruction_pipeline(input_file, output_file, *, attributes_to_interpolat iso_surface_threshold=0.6, mesh_smoothing_weights=False, output_mesh_smoothing_weights=False, sph_normals=False, mesh_smoothing_weights_normalization=13.0, mesh_smoothing_iters=5, normals_smoothing_iters=5, mesh_aabb_min=None, mesh_aabb_max=None, mesh_cleanup=False, decimate_barnacles=False, keep_vertices=False, - compute_normals=False, output_raw_normals=False, mesh_aabb_clamp_vertices=False, - check_mesh_closed=False, check_mesh_manifold=False, check_mesh_debug=False, + compute_normals=False, output_raw_normals=False, output_raw_mesh=False, mesh_aabb_clamp_vertices=False, + check_mesh_closed=False, check_mesh_manifold=False, check_mesh_orientation=False, check_mesh_debug=False, generate_quads=False, quad_max_edge_diag_ratio=1.75, quad_max_normal_angle=10.0, quad_max_interior_angle=135.0, subdomain_grid=False, subdomain_num_cubes_per_dim=64): @@ -106,25 +106,12 @@ def reconstruction_pipeline(input_file, output_file, *, attributes_to_interpolat mesh_smoothing_weights_normalization=mesh_smoothing_weights_normalization, mesh_smoothing_iters=mesh_smoothing_iters, normals_smoothing_iters=normals_smoothing_iters, mesh_aabb_min=mesh_aabb_min, mesh_aabb_max=mesh_aabb_max, mesh_cleanup=mesh_cleanup, decimate_barnacles=decimate_barnacles, - keep_vertices=keep_vertices, compute_normals=compute_normals, output_raw_normals=output_raw_normals, - mesh_aabb_clamp_vertices=mesh_aabb_clamp_vertices, subdomain_grid=subdomain_grid, subdomain_num_cubes_per_dim=subdomain_num_cubes_per_dim, output_mesh_smoothing_weights=output_mesh_smoothing_weights) - - # Convert triangles to quads - if generate_quads: - mesh_with_data = pysplashsurf.convert_tris_to_quads(mesh_with_data, non_squareness_limit=quad_max_edge_diag_ratio, normal_angle_limit_rad=math.radians(quad_max_normal_angle), max_interior_angle=math.radians(quad_max_interior_angle)) + keep_vertices=keep_vertices, compute_normals=compute_normals, output_raw_normals=output_raw_normals, output_raw_mesh=output_raw_mesh, + mesh_aabb_clamp_vertices=mesh_aabb_clamp_vertices, subdomain_grid=subdomain_grid, subdomain_num_cubes_per_dim=subdomain_num_cubes_per_dim, output_mesh_smoothing_weights=output_mesh_smoothing_weights, + check_mesh_closed=check_mesh_closed, check_mesh_manifold=check_mesh_manifold, check_mesh_orientation=check_mesh_orientation, check_mesh_debug=check_mesh_debug, + generate_quads=generate_quads, quad_max_edge_diag_ratio=quad_max_edge_diag_ratio, quad_max_normal_angle=quad_max_normal_angle, quad_max_interior_angle=quad_max_interior_angle) - pysplashsurf.write_to_file(mesh_with_data, output_file, consume_object=False) - - mesh = mesh_with_data.take_mesh() - - if type(mesh) is pysplashsurf.TriMesh3dF64 or type(mesh) is pysplashsurf.TriMesh3dF32: - # Mesh checks - if check_mesh_closed or check_mesh_manifold: - pysplashsurf.check_mesh_consistency(reconstruction.grid, mesh, check_closed=check_mesh_closed, check_manifold=check_mesh_manifold, debug=check_mesh_debug) - - - # Left out: Mesh orientation check - + pysplashsurf.write_to_file(mesh_with_data, output_file, consume_object=True) def test_no_post_processing(): @@ -134,9 +121,8 @@ def test_no_post_processing(): start = time.time() reconstruction_pipeline(VTK_PATH, DIR.joinpath("test.vtk"), particle_radius=np.float64(0.025), smoothing_length=np.float64(2.0), - cube_size=np.float64(0.5), iso_surface_threshold=np.float64(0.6), mesh_smoothing_weights=True, - mesh_smoothing_weights_normalization=np.float64(13.0), mesh_smoothing_iters=0, normals_smoothing_iters=0, - generate_quads=False, mesh_cleanup=False, compute_normals=False, subdomain_grid=True) + cube_size=np.float64(0.5), iso_surface_threshold=np.float64(0.6), mesh_smoothing_weights=False, + mesh_smoothing_iters=0, normals_smoothing_iters=0, mesh_cleanup=False, compute_normals=False, subdomain_grid=True) print("Python done in", time.time() - start) binary_mesh = meshio.read(DIR.joinpath("test_bin.vtk")) diff --git a/splashsurf/src/cli.rs b/splashsurf/src/cli.rs index f0404bc6..05c8a136 100644 --- a/splashsurf/src/cli.rs +++ b/splashsurf/src/cli.rs @@ -4,7 +4,7 @@ //! The reconstruction procedure and other internals of the CLI are provided by the [`splashsurf_lib`] crate. use crate::allocator::GetPeakAllocatedMemory; -use crate::{convert, logging, reconstruction}; +use crate::{convert, logging, reconstruct}; use anyhow::Context; use clap::Parser; use log::info; @@ -44,7 +44,7 @@ struct CommandlineArgs { enum Subcommand { /// Reconstruct a surface from particle data #[command(help_template = HELP_TEMPLATE)] - Reconstruct(reconstruction::ReconstructSubcommandArgs), + Reconstruct(reconstruct::ReconstructSubcommandArgs), /// Convert particle or mesh files between different file formats #[command(help_template = HELP_TEMPLATE)] Convert(convert::ConvertSubcommandArgs), @@ -82,7 +82,9 @@ impl Switch { /// Runs the splashsurf CLI with the provided command line arguments. /// -/// This function behaves like the binary `splashsurf` command line tool including output to stdout and stderr. +/// This function behaves like the binary `splashsurf` command line tool including output to stdout +/// and stderr. It will also exit the process depending on the command line arguments, so it should +/// not be used in typical library contexts. /// Note that the first argument is always ignored - this is typically the binary name when called using /// `std::env::args()` from the terminal: /// ``` @@ -115,7 +117,7 @@ where // Delegate to subcommands let result = match &cmd_args.subcommand { - Subcommand::Reconstruct(cmd_args) => reconstruction::reconstruct_subcommand(cmd_args), + Subcommand::Reconstruct(cmd_args) => reconstruct::reconstruct_subcommand(cmd_args), Subcommand::Convert(cmd_args) => convert::convert_subcommand(cmd_args), }; @@ -189,7 +191,7 @@ mod cli_args_tests { #[test] fn verify_reconstruct_cli() { use clap::CommandFactory; - crate::reconstruction::ReconstructSubcommandArgs::command().debug_assert() + crate::reconstruct::ReconstructSubcommandArgs::command().debug_assert() } #[test] diff --git a/splashsurf/src/lib.rs b/splashsurf/src/lib.rs index 05a5b34f..c20f824d 100644 --- a/splashsurf/src/lib.rs +++ b/splashsurf/src/lib.rs @@ -1,7 +1,24 @@ +//! Library target of the `splashsurf` CLI. +//! +//! To use the CLI you can install it via Cargo: +//! ```bash +//! cargo install splashsurf +//! ``` +//! For documentation of the CLI see the [README](https://github.com/InteractiveComputerGraphics/splashsurf) in the project repository. +//! +//! This library target exposes some high-level functionality such as [running the CLI](cli) with a +//! list of command line arguments or specifically running the [reconstruction pipeline](reconstruct) +//! of the CLI (including a set of postprocessing steps) as a library function. +//! This functionality is mainly used by the `pySplashsurf` Pythong bindings and the CLI binary +//! target itself. +//! +//! If you only want to use a subset of this functionality (e.g. only the reconstruction itself, +//! without any postprocessing) in your crates, please refer to the [`splashsurf_lib`] crate instead. + pub mod cli; mod convert; mod io; -mod reconstruction; +pub mod reconstruct; #[macro_use] mod allocator; mod logging; diff --git a/splashsurf/src/main.rs b/splashsurf/src/main.rs index 66b5152d..b16ca167 100644 --- a/splashsurf/src/main.rs +++ b/splashsurf/src/main.rs @@ -1,7 +1,7 @@ pub mod cli; mod convert; mod io; -mod reconstruction; +mod reconstruct; #[macro_use] mod allocator; mod logging; diff --git a/splashsurf/src/reconstruction.rs b/splashsurf/src/reconstruct.rs similarity index 86% rename from splashsurf/src/reconstruction.rs rename to splashsurf/src/reconstruct.rs index c48b306d..a505a798 100644 --- a/splashsurf/src/reconstruction.rs +++ b/splashsurf/src/reconstruct.rs @@ -1,15 +1,19 @@ +//! Implementation of the `reconstruct` subcommand of the splashsurf CLI. + use crate::cli::Switch; -use crate::reconstruction::arguments::*; +use crate::reconstruct::arguments::*; use crate::{io, logging}; use anyhow::{Context, anyhow}; use clap::value_parser; use indicatif::{ProgressBar, ProgressStyle}; use log::{error, info, warn}; use rayon::prelude::*; -use splashsurf_lib::mesh::{AttributeData, Mesh3d, MeshAttribute, MeshWithData}; +use splashsurf_lib::mesh::{ + AttributeData, Mesh3d, MeshAttribute, MeshWithData, MixedTriQuadMesh3d, TriMesh3d, +}; use splashsurf_lib::nalgebra::{Unit, Vector3}; use splashsurf_lib::sph_interpolation::SphInterpolator; -use splashsurf_lib::{Aabb3d, Index, Real, profile}; +use splashsurf_lib::{Aabb3d, Index, Real, SurfaceReconstruction, profile}; use std::borrow::Cow; use std::collections::HashMap; use std::convert::TryFrom; @@ -32,7 +36,7 @@ static ARGS_OTHER: &str = "Remaining options"; #[derive(Clone, Debug, clap::Parser)] #[clap(group = clap::ArgGroup::new("input").required(true))] #[command(next_help_heading = ARGS_OTHER)] -pub struct ReconstructSubcommandArgs { +pub(crate) struct ReconstructSubcommandArgs { /// Path to the input file where the particle positions are stored (supported formats: VTK 4.2, VTU, binary f32 XYZ, PLY, BGEO), use "{}" in the filename to indicate a placeholder for a sequence. #[arg(help_heading = ARGS_IO, group = "input", value_parser = value_parser!(PathBuf))] pub input_file_or_sequence: PathBuf, @@ -355,7 +359,9 @@ pub struct ReconstructSubcommandArgs { } /// Executes the `reconstruct` subcommand -pub fn reconstruct_subcommand(cmd_args: &ReconstructSubcommandArgs) -> Result<(), anyhow::Error> { +pub(crate) fn reconstruct_subcommand( + cmd_args: &ReconstructSubcommandArgs, +) -> Result<(), anyhow::Error> { profile!("reconstruct subcommand"); let paths = ReconstructionRunnerPathCollection::try_from(cmd_args) @@ -377,7 +383,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_pipeline(path, &args) + reconstruction_pipeline_from_args(path, &args) .with_context(|| { format!( "Error while processing input file \"{}\" from a file sequence", @@ -396,7 +402,7 @@ pub fn reconstruct_subcommand(cmd_args: &ReconstructSubcommandArgs) -> Result<() }) } else { paths.iter().try_for_each(|path| { - reconstruction_pipeline(path, &args).map(|_| { + reconstruction_pipeline_from_args(path, &args).map(|_| { if let Some(pb) = logging::get_progress_bar() { pb.inc(1) } @@ -418,9 +424,73 @@ pub fn reconstruct_subcommand(cmd_args: &ReconstructSubcommandArgs) -> Result<() result } +/// Struct to hold the result of the reconstruction pipeline +pub struct ReconstructionResult { + /// Holds the reconstructed triangle mesh with optional data (if `generate_quads` was not enabled) + pub tri_mesh: Option>>, + /// Holds the reconstructed quad mesh with optional data (only if `generate_quads` was enabled) + pub tri_quad_mesh: Option>>, + /// Holds the surface reconstruction with no post-processing applied if `output_raw_mesh` was enabled + pub raw_reconstruction: Option>, +} + +/// Parameters for the post-processing steps in the reconstruction pipeline +#[derive(Clone, Debug)] +pub struct ReconstructionPostprocessingParameters { + /// Enable checking the final mesh for holes + pub check_mesh_closed: bool, + /// Enable checking the final mesh for non-manifold edges and vertices + pub check_mesh_manifold: bool, + /// Enable checking the final mesh for inverted triangles (compares angle between vertex normals and adjacent face normals) + pub check_mesh_orientation: bool, + /// Enable additional debug output for the check-mesh operations (has no effect if no other check-mesh option is enabled) + pub check_mesh_debug: bool, + /// Enable MC specific mesh decimation/simplification which removes bad quality triangles typically generated by MC by snapping (enabled by default if smoothing is enabled) + pub mesh_cleanup: bool, + /// If MC mesh cleanup is enabled, vertex snapping can be limited to this distance relative to the MC edge length (should be in the interval `[0.0,0.5]`) + pub mesh_cleanup_snap_dist: Option, + /// Enable decimation of some typical bad marching cubes triangle configurations (resulting in "barnacles" after Laplacian smoothing) + pub decimate_barnacles: bool, + /// Enable preserving vertices without connectivity during decimation instead of filtering them out (faster and helps with debugging) + pub keep_vertices: bool, + /// Enable computing surface normals at the mesh vertices and write them to the output object + pub compute_normals: bool, + /// Enable computing the normals using SPH interpolation instead of using the area weighted triangle normals + pub sph_normals: bool, + /// Number of smoothing iterations to apply to normals if normal interpolation is enabled + pub normals_smoothing_iters: Option, + /// Interpolate point attributes with the given name from the input attributes to the reconstructed surface + pub interpolate_attributes: Vec, + /// Number of smoothing iterations to run on the reconstructed mesh + pub mesh_smoothing_iters: Option, + /// Enable feature weights for mesh smoothing if mesh smoothing enabled. Preserves isolated particles even under strong smoothing. + pub mesh_smoothing_weights: bool, + /// Override a manual normalization value from weighted number of neighbors to mesh smoothing weights + pub mesh_smoothing_weights_normalization: f64, + /// Enable conversion of triangles to quads if they meet quality criteria + pub generate_quads: bool, + /// Maximum allowed ratio of quad edge lengths to its diagonals to merge two triangles to a quad (inverse is used for minimum) + pub quad_max_edge_diag_ratio: f64, + /// Maximum allowed angle (in degrees) between triangle normals to merge them to a quad + pub quad_max_normal_angle: f64, + /// Maximum allowed vertex interior angle (in degrees) inside a quad to merge two triangles to a quad + pub quad_max_interior_angle: f64, + /// Enable writing the smoothing weights as a vertex attribute to the output mesh file + pub output_mesh_smoothing_weights: bool, + /// Enable writing raw normals without smoothing to the output mesh if normal smoothing is enabled + pub output_raw_normals: bool, + // For reconstruction_pipeline_from_args/path: Enable writing the raw reconstructed mesh before applying any post-processing steps (like smoothing or decimation) + /// When true, also return the SurfaceReconstruction object with no post-processing applied + pub output_raw_mesh: bool, + /// Bounding-box for the surface mesh, triangles completely outside are removed + pub mesh_aabb: Option>, + /// Enable clamping of vertices outside the specified mesh AABB to the AABB (only has an effect if mesh-aabb is specified) + pub mesh_aabb_clamp_vertices: bool, +} + /// Conversion and validation of command line arguments -mod arguments { - use super::ReconstructSubcommandArgs; +pub(crate) mod arguments { + use super::{ReconstructSubcommandArgs, ReconstructionPostprocessingParameters}; use crate::io; use anyhow::{Context, anyhow}; use log::info; @@ -433,40 +503,13 @@ mod arguments { use std::str::FromStr; use walkdir::WalkDir; - pub struct ReconstructionRunnerPostprocessingArgs { - pub check_mesh_closed: bool, - pub check_mesh_manifold: bool, - pub check_mesh_orientation: bool, - pub check_mesh_debug: bool, - pub mesh_cleanup: bool, - pub mesh_cleanup_snap_dist: Option, - pub decimate_barnacles: bool, - pub keep_vertices: bool, - pub compute_normals: bool, - pub sph_normals: bool, - pub normals_smoothing_iters: Option, - pub interpolate_attributes: Vec, - pub mesh_smoothing_iters: Option, - pub mesh_smoothing_weights: bool, - pub mesh_smoothing_weights_normalization: f64, - pub generate_quads: bool, - pub quad_max_edge_diag_ratio: f64, - pub quad_max_normal_angle: f64, - pub quad_max_interior_angle: f64, - pub output_mesh_smoothing_weights: bool, - pub output_raw_normals: bool, - pub output_raw_mesh: bool, - pub mesh_aabb: Option>, - pub mesh_aabb_clamp_vertices: bool, - } - /// All arguments that can be supplied to the surface reconstruction tool converted to useful types pub struct ReconstructionRunnerArgs { /// Parameters passed directly to the surface reconstruction pub params: splashsurf_lib::Parameters, pub use_double_precision: bool, pub io_params: io::FormatParameters, - pub postprocessing: ReconstructionRunnerPostprocessingArgs, + pub postprocessing: ReconstructionPostprocessingParameters, } fn try_aabb_from_min_max( @@ -562,7 +605,7 @@ mod arguments { splashsurf_lib::initialize_thread_pool(num_threads)?; } - let postprocessing = ReconstructionRunnerPostprocessingArgs { + let postprocessing = ReconstructionPostprocessingParameters { check_mesh_closed: args.check_mesh.into_bool() || args.check_mesh_closed.into_bool(), check_mesh_manifold: args.check_mesh.into_bool() @@ -883,13 +926,13 @@ mod arguments { } /// Calls the reconstruction pipeline for single or double precision depending on the runtime parameters -pub(crate) fn reconstruction_pipeline( +pub(crate) fn reconstruction_pipeline_from_args( paths: &ReconstructionRunnerPaths, args: &ReconstructionRunnerArgs, ) -> Result<(), anyhow::Error> { if args.use_double_precision { info!("Using double precision (f64) for surface reconstruction."); - reconstruction_pipeline_generic::( + reconstruction_pipeline_from_path::( paths, &args.params, &args.io_params, @@ -897,7 +940,7 @@ pub(crate) fn reconstruction_pipeline( )?; } else { info!("Using single precision (f32) for surface reconstruction."); - reconstruction_pipeline_generic::( + reconstruction_pipeline_from_path::( paths, &args.params.try_convert().ok_or(anyhow!( "Unable to convert surface reconstruction parameters from f64 to f32." @@ -910,63 +953,33 @@ pub(crate) fn reconstruction_pipeline( Ok(()) } -/// Wrapper for the reconstruction pipeline: loads input file, runs reconstructions, stores output files -pub(crate) fn reconstruction_pipeline_generic( - paths: &ReconstructionRunnerPaths, +/// Performs a surface reconstruction including optional post-processing steps +/// +/// This function implements the surface reconstruction pipeline used by the `reconstruct` subcommand +/// of the `splashsurf` CLI. +/// Inputs are the particle positions, a (possibly empty) list of attributes defined on the particles, +/// [`Parameters`](splashsurf_lib::Parameters) for the surface reconstruction itself, and a set of parameters for optional +/// post-processing steps. +/// Please note that, unlike the CLI, the parameters for the surface reconstruction are not relative +/// to the particle radius but absolute values. +pub fn reconstruction_pipeline( + particle_positions: &[Vector3], + attributes: Vec>, params: &splashsurf_lib::Parameters, - io_params: &io::FormatParameters, - postprocessing: &ReconstructionRunnerPostprocessingArgs, -) -> Result<(), anyhow::Error> { - profile!("surface reconstruction"); - - // Load particle positions and attributes to interpolate - let (particle_positions, attributes) = io::read_particle_positions_with_attributes( - &paths.input_file, - &postprocessing.interpolate_attributes, - &io_params.input, - ) - .with_context(|| { - format!( - "Failed to load particle positions from file \"{}\"", - paths.input_file.display() - ) - })?; - + postprocessing: &ReconstructionPostprocessingParameters, +) -> Result, anyhow::Error> { // Perform the surface reconstruction - let reconstruction = - splashsurf_lib::reconstruct_surface::(particle_positions.as_slice(), params)?; + let reconstruction = splashsurf_lib::reconstruct_surface::(particle_positions, params)?; + + let reconstruction_output = if postprocessing.output_raw_mesh { + Some(reconstruction.clone()) + } else { + None + }; let grid = reconstruction.grid(); let mut mesh_with_data = MeshWithData::new(Cow::Borrowed(reconstruction.mesh())); - if postprocessing.output_raw_mesh { - profile!("write surface mesh to file"); - - let output_path = paths - .output_file - .parent() - // Add a trailing separator if the parent is non-empty - .map(|p| p.join("")) - .unwrap_or_default(); - let output_filename = format!( - "raw_{}", - paths.output_file.file_name().unwrap().to_string_lossy() - ); - let raw_output_file = output_path.join(output_filename); - - info!( - "Writing unprocessed surface mesh to \"{}\"...", - raw_output_file.display() - ); - - io::write_mesh(&mesh_with_data, raw_output_file, &io_params.output).with_context(|| { - anyhow!( - "Failed to write raw output mesh to file \"{}\"", - paths.output_file.display() - ) - })?; - } - // Perform post-processing { profile!("postprocessing"); @@ -1064,13 +1077,13 @@ pub(crate) fn reconstruction_pipeline_generic( { let search_radius = params.compact_support_radius; - let mut domain = Aabb3d::from_points(particle_positions.as_slice()); + let mut domain = Aabb3d::from_points(particle_positions); domain.grow_uniformly(search_radius); let mut nl = Vec::new(); splashsurf_lib::neighborhood_search::neighborhood_search_spatial_hashing_parallel::( &domain, - particle_positions.as_slice(), + particle_positions, search_radius, &mut nl, ); @@ -1277,7 +1290,7 @@ pub(crate) fn reconstruction_pipeline_generic( }; // Convert triangles to quads - let (tri_mesh, tri_quad_mesh) = if postprocessing.generate_quads { + let (mut tri_mesh, tri_quad_mesh) = if postprocessing.generate_quads { info!("Post-processing: Convert triangles to quads..."); let non_squareness_limit = R::from_f64(postprocessing.quad_max_edge_diag_ratio).unwrap(); let normal_angle_limit_rad = @@ -1312,28 +1325,6 @@ pub(crate) fn reconstruction_pipeline_generic( (Some(mesh_with_data), None) }; - // Store the surface mesh - { - profile!("write surface mesh to file"); - - match (&tri_mesh, &tri_quad_mesh) { - (Some(mesh), None) => { - io::write_mesh(mesh, paths.output_file.clone(), &io_params.output) - } - (None, Some(mesh)) => { - io::write_mesh(mesh, paths.output_file.clone(), &io_params.output) - } - - _ => unreachable!(), - } - .with_context(|| { - anyhow!( - "Failed to write output mesh to file \"{}\"", - paths.output_file.display() - ) - })?; - } - // TODO: Option to continue processing sequences even if checks fail. Maybe return special error type? if postprocessing.check_mesh_closed || postprocessing.check_mesh_manifold { @@ -1347,7 +1338,11 @@ pub(crate) fn reconstruction_pipeline_generic( ), (None, Some(_mesh)) => { info!("Checking for mesh consistency not implemented for quad mesh at the moment."); - return Ok(()); + return Ok(ReconstructionResult { + tri_mesh: None, + tri_quad_mesh: Some(_mesh.to_owned()), + raw_reconstruction: reconstruction_output, + }); } _ => unreachable!(), } { @@ -1358,7 +1353,7 @@ pub(crate) fn reconstruction_pipeline_generic( error!("{}", err); return Err(anyhow!("{}", err)) .context(format!("Checked mesh for problems (holes: {}, non-manifold edges/vertices: {}), problems were found!", postprocessing.check_mesh_closed, postprocessing.check_mesh_manifold)) - .context(format!("Problem found with mesh file \"{}\"", paths.output_file.display())); + .context("Problem found with mesh"); } else { info!( "Checked mesh for problems (holes: {}, non-manifold edges/vertices: {}), no problems were found.", @@ -1416,7 +1411,7 @@ pub(crate) fn reconstruction_pipeline_generic( info!( "Checking for normal orientation not implemented for quad mesh at the moment." ); - return Ok(()); + Ok(()) } _ => unreachable!(), } { @@ -1424,14 +1419,121 @@ pub(crate) fn reconstruction_pipeline_generic( error!("{}", err); return Err(anyhow!("{}", err)) .context("Checked mesh orientation (flipped normals), problems were found!") - .context(format!( - "Problem found with mesh file \"{}\"", - paths.output_file.display() - )); + .context("Problem found with mesh"); } else { info!("Checked mesh orientation (flipped normals), no problems were found."); } } + match (&mut tri_mesh, &tri_quad_mesh) { + (Some(mesh), None) => { + let mut res: MeshWithData> = + MeshWithData::new(mesh.to_owned().mesh.into_owned()); + res.point_attributes = std::mem::take(&mut mesh.point_attributes); + res.cell_attributes = std::mem::take(&mut mesh.cell_attributes); + + Ok(ReconstructionResult { + tri_mesh: Some(res), + tri_quad_mesh: None, + raw_reconstruction: reconstruction_output, + }) + } + (None, Some(_mesh)) => Ok(ReconstructionResult { + tri_mesh: None, + tri_quad_mesh: Some(_mesh.to_owned()), + raw_reconstruction: reconstruction_output, + }), + _ => unreachable!(), + } +} + +/// Wrapper for the reconstruction pipeline: loads input file, runs reconstructions, stores output files +pub(crate) fn reconstruction_pipeline_from_path( + paths: &ReconstructionRunnerPaths, + params: &splashsurf_lib::Parameters, + io_params: &io::FormatParameters, + postprocessing: &ReconstructionPostprocessingParameters, +) -> Result<(), anyhow::Error> { + profile!("surface reconstruction"); + + // Load particle positions and attributes to interpolate + let (particle_positions, attributes) = io::read_particle_positions_with_attributes( + &paths.input_file, + &postprocessing.interpolate_attributes, + &io_params.input, + ) + .with_context(|| { + format!( + "Failed to load particle positions from file \"{}\"", + paths.input_file.display() + ) + })?; + + let ReconstructionResult { + tri_mesh, + tri_quad_mesh, + raw_reconstruction: reconstruction, + } = reconstruction_pipeline::(&particle_positions, attributes, params, postprocessing)?; + + if postprocessing.output_raw_mesh { + profile!("write surface mesh to file"); + + let reconstruction = reconstruction.expect( + "reconstruction_pipeline_from_data did not return a SurfaceReconstruction object", + ); + let mesh = reconstruction.mesh(); + + let output_path = paths + .output_file + .parent() + // Add a trailing separator if the parent is non-empty + .map(|p| p.join("")) + .unwrap_or_default(); + let output_filename = format!( + "raw_{}", + paths.output_file.file_name().unwrap().to_string_lossy() + ); + let raw_output_file = output_path.join(output_filename); + + info!( + "Writing unprocessed surface mesh to \"{}\"...", + raw_output_file.display() + ); + + io::write_mesh( + &MeshWithData::new(mesh.to_owned()), + raw_output_file, + &io_params.output, + ) + .with_context(|| { + anyhow!( + "Failed to write raw output mesh to file \"{}\"", + paths.output_file.display() + ) + })?; + } + + // Store the surface mesh + { + profile!("write surface mesh to file"); + + match (&tri_mesh, &tri_quad_mesh) { + (Some(mesh), None) => { + io::write_mesh(mesh, paths.output_file.clone(), &io_params.output) + } + (None, Some(mesh)) => { + io::write_mesh(mesh, paths.output_file.clone(), &io_params.output) + } + + _ => unreachable!(), + } + .with_context(|| { + anyhow!( + "Failed to write output mesh to file \"{}\"", + paths.output_file.display() + ) + })?; + } + Ok(()) }