use alloc::vec::Vec;
use core::fmt;
use crate::{MeshBuffer, Real};
mod accuracy;
mod determinism;
mod field_bound;
mod isotopy;
mod self_intersection;
mod tri_grid;
pub use accuracy::{AccuracyConfig, AccuracyReport, DistanceStats, accuracy};
pub use determinism::{DeterminismReport, Divergence, RunPair, check_determinism};
pub use field_bound::{EIKONAL_TOLERANCE, FieldBoundReport, field_bound_report};
pub use isotopy::{IsotopyReport, cell_is_certified, isotopy_report};
pub use self_intersection::{SelfIntersectionReport, self_intersections};
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct ValidateConfig {
cell_size: f64,
weld_epsilon: f64,
area_epsilon_rel: f64,
}
impl ValidateConfig {
pub const WELD_EPSILON_REL: f64 = 1e-4;
pub const AREA_EPSILON_REL: f64 = 1e-6;
pub fn from_cell_size(cell_size: f64) -> crate::Result<Self> {
if !cell_size.is_finite() || cell_size <= 0.0 {
return Err(crate::Error::InvalidCellSize { value: cell_size });
}
Ok(Self {
cell_size,
weld_epsilon: cell_size * Self::WELD_EPSILON_REL,
area_epsilon_rel: Self::AREA_EPSILON_REL,
})
}
#[must_use]
pub fn cell_size(&self) -> f64 {
self.cell_size
}
#[must_use]
pub fn weld_epsilon(&self) -> f64 {
self.weld_epsilon
}
#[must_use]
pub fn area_epsilon_rel(&self) -> f64 {
self.area_epsilon_rel
}
}
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct MeshReport {
pub vertices: u64,
pub referenced_vertices: u64,
pub edges: u64,
pub faces: u64,
pub faces_skipped: u64,
pub euler_characteristic: i64,
pub components: u64,
pub boundary_loops: u64,
pub genus: Option<i64>,
pub non_manifold_edges: u64,
pub non_manifold_vertices: u64,
pub boundary_edges: u64,
pub inconsistently_oriented_edges: u64,
pub degenerate_triangles: u64,
pub repeated_index_triangles: u64,
pub weld_buckets: u64,
pub duplicate_vertices: u64,
pub unreferenced_vertices: u64,
pub out_of_range_indices: u64,
pub trailing_indices: u64,
pub non_finite_positions: u64,
pub normal_count_mismatch: bool,
pub config: ValidateConfig,
}
impl MeshReport {
#[must_use]
pub fn has_structural_errors(&self) -> bool {
self.out_of_range_indices > 0
|| self.trailing_indices > 0
|| self.repeated_index_triangles > 0
|| self.non_finite_positions > 0
|| self.normal_count_mismatch
}
#[must_use]
pub fn is_manifold(&self) -> bool {
!self.has_structural_errors()
&& self.non_manifold_edges == 0
&& self.non_manifold_vertices == 0
&& self.inconsistently_oriented_edges == 0
}
#[must_use]
pub fn is_closed(&self) -> bool {
self.is_manifold() && self.boundary_edges == 0 && self.euler_characteristic % 2 == 0
}
#[must_use]
pub fn satisfies(&self, gate: SurfaceGate) -> bool {
if self.has_structural_errors() {
return false;
}
match gate {
SurfaceGate::Closed => self.is_closed(),
SurfaceGate::Manifold => self.is_manifold(),
SurfaceGate::ClosedAllowingUnresolvedTopology => {
self.boundary_edges == 0 && self.inconsistently_oriented_edges == 0
}
}
}
#[must_use]
pub fn violations(&self) -> u64 {
self.non_manifold_edges
+ self.non_manifold_vertices
+ self.inconsistently_oriented_edges
+ self.repeated_index_triangles
+ self.out_of_range_indices
+ self.trailing_indices
+ self.non_finite_positions
+ u64::from(self.normal_count_mismatch)
}
pub fn panic_if_invalid(&self, expect_closed: bool) {
let ok = if expect_closed {
self.is_closed()
} else {
self.is_manifold()
};
assert!(ok, "{self}");
}
}
impl fmt::Display for MeshReport {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
let partial = self.faces_skipped > 0 || self.trailing_indices > 0;
let rule = " ------------------------------------------------------------";
writeln!(f, "mesh report")?;
writeln!(
f,
" vertices {:8} ({} referenced, {} unreferenced)",
self.vertices, self.referenced_vertices, self.unreferenced_vertices
)?;
writeln!(f, " edges {:8}", self.edges)?;
if partial {
writeln!(
f,
" faces {:8} ({} considered, {} skipped)",
self.faces + self.faces_skipped,
self.faces,
self.faces_skipped
)?;
writeln!(
f,
" euler characteristic {:8} (over the valid subset only)",
self.euler_characteristic
)?;
} else {
writeln!(f, " faces {:8}", self.faces)?;
writeln!(
f,
" euler characteristic {:8}",
self.euler_characteristic
)?;
}
writeln!(f, " components {:8}", self.components)?;
writeln!(f, " boundary loops {:8}", self.boundary_loops)?;
match self.genus {
Some(g) => writeln!(f, " genus {g:8}")?,
None => writeln!(
f,
" genus - (not a single oriented manifold component)"
)?,
}
writeln!(f, "{rule}")?;
writeln!(
f,
" non-manifold edges {:8}",
self.non_manifold_edges
)?;
writeln!(
f,
" non-manifold vertices {:8}",
self.non_manifold_vertices
)?;
writeln!(f, " boundary edges {:8}", self.boundary_edges)?;
writeln!(
f,
" inconsistently oriented {:8}",
self.inconsistently_oriented_edges
)?;
writeln!(
f,
" degenerate triangles {:8} (area <= {:e} * h^2, h = {})",
self.degenerate_triangles, self.config.area_epsilon_rel, self.config.cell_size
)?;
writeln!(
f,
" repeated-index triangles {:8}",
self.repeated_index_triangles
)?;
writeln!(
f,
" duplicate vertices {:8} (within {:e})",
self.duplicate_vertices, self.config.weld_epsilon
)?;
writeln!(
f,
" out-of-range indices {:8}",
self.out_of_range_indices
)?;
writeln!(f, " trailing indices {:8}", self.trailing_indices)?;
writeln!(
f,
" non-finite positions {:8}",
self.non_finite_positions
)?;
if self.normal_count_mismatch {
writeln!(f, " normal count mismatch yes")?;
}
writeln!(f, "{rule}")?;
if self.has_structural_errors() {
writeln!(
f,
" !! STRUCTURAL ERRORS - derived metrics cover the valid subset only"
)?;
}
if self.is_closed() {
write!(f, " MANIFOLD, CLOSED")
} else if self.is_manifold() {
write!(f, " MANIFOLD, WITH BOUNDARY")
} else {
write!(f, " INVALID: {} violations", self.violations())
}
}
}
struct Dsu {
parent: Vec<u32>,
size: Vec<u32>,
}
impl Dsu {
fn new(n: usize) -> Self {
Self {
parent: (0..n as u32).collect(),
size: alloc::vec![1; n],
}
}
fn reset(&mut self, n: usize) {
self.parent.clear();
self.parent.extend(0..n as u32);
self.size.clear();
self.size.resize(n, 1);
}
fn find(&mut self, mut x: u32) -> u32 {
while self.parent[x as usize] != x {
let grand = self.parent[self.parent[x as usize] as usize];
self.parent[x as usize] = grand;
x = grand;
}
x
}
fn union(&mut self, a: u32, b: u32) {
let (mut ra, mut rb) = (self.find(a), self.find(b));
if ra == rb {
return;
}
if self.size[ra as usize] < self.size[rb as usize] {
core::mem::swap(&mut ra, &mut rb);
}
self.parent[rb as usize] = ra;
self.size[ra as usize] += self.size[rb as usize];
}
}
#[derive(Clone, Debug, Default, PartialEq, Eq)]
pub struct NonManifoldFeatures {
pub edges: Vec<[u32; 2]>,
pub vertices: Vec<u32>,
pub boundary_edges: Vec<[u32; 2]>,
pub inconsistently_oriented_edges: Vec<[u32; 2]>,
}
#[must_use]
pub fn validate_indexed<R: Real>(
positions: &[[R; 3]],
indices: &[u32],
cfg: &ValidateConfig,
) -> MeshReport {
validate_features(positions, indices, cfg).0
}
pub fn validate_features<R: Real>(
positions: &[[R; 3]],
indices: &[u32],
cfg: &ValidateConfig,
) -> (MeshReport, NonManifoldFeatures) {
let mut features = NonManifoldFeatures::default();
let vertex_count = positions.len();
let mut report = MeshReport {
vertices: vertex_count as u64,
referenced_vertices: 0,
edges: 0,
faces: 0,
faces_skipped: 0,
euler_characteristic: 0,
components: 0,
boundary_loops: 0,
genus: None,
non_manifold_edges: 0,
non_manifold_vertices: 0,
boundary_edges: 0,
inconsistently_oriented_edges: 0,
degenerate_triangles: 0,
repeated_index_triangles: 0,
duplicate_vertices: 0,
weld_buckets: 0,
unreferenced_vertices: 0,
out_of_range_indices: 0,
trailing_indices: (indices.len() % 3) as u64,
non_finite_positions: 0,
normal_count_mismatch: false,
config: *cfg,
};
for p in positions {
if !(p[0].is_finite() && p[1].is_finite() && p[2].is_finite()) {
report.non_finite_positions += 1;
}
}
let whole = indices.len() - indices.len() % 3;
let mut faces: Vec<[u32; 3]> = Vec::with_capacity(whole / 3);
for tri in indices[..whole].chunks_exact(3) {
let mut out_of_range = 0u64;
for &i in tri {
if i as usize >= vertex_count {
out_of_range += 1;
}
}
let repeated = tri[0] == tri[1] || tri[1] == tri[2] || tri[0] == tri[2];
if repeated {
report.repeated_index_triangles += 1;
}
if out_of_range > 0 || repeated {
report.out_of_range_indices += out_of_range;
report.faces_skipped += 1;
continue;
}
faces.push([tri[0], tri[1], tri[2]]);
}
report.faces = faces.len() as u64;
let mut referenced = alloc::vec![false; vertex_count];
for f in &faces {
for &v in f {
referenced[v as usize] = true;
}
}
report.referenced_vertices = referenced.iter().filter(|r| **r).count() as u64;
report.unreferenced_vertices = report.vertices - report.referenced_vertices;
let mut half_edges: Vec<(u32, u32, u32, u8)> = Vec::with_capacity(faces.len() * 3);
for (fi, f) in faces.iter().enumerate() {
for k in 0..3 {
let (a, b) = (f[k], f[(k + 1) % 3]);
let (lo, hi, forward) = if a < b { (a, b, 1) } else { (b, a, 0) };
half_edges.push((lo, hi, fi as u32, forward));
}
}
half_edges.sort_unstable();
let mut face_dsu = Dsu::new(faces.len());
let mut boundary_dsu = Dsu::new(vertex_count);
let mut on_boundary = alloc::vec![false; vertex_count];
let mut start = 0usize;
while start < half_edges.len() {
let (lo, hi, _, _) = half_edges[start];
let mut end = start + 1;
while end < half_edges.len() && half_edges[end].0 == lo && half_edges[end].1 == hi {
end += 1;
}
let run = &half_edges[start..end];
report.edges += 1;
match run.len() {
1 => {
report.boundary_edges += 1;
features.boundary_edges.push([lo, hi]);
on_boundary[lo as usize] = true;
on_boundary[hi as usize] = true;
boundary_dsu.union(lo, hi);
}
2 => {
if run[0].3 == run[1].3 {
report.inconsistently_oriented_edges += 1;
features.inconsistently_oriented_edges.push([lo, hi]);
}
}
_ => {
report.non_manifold_edges += 1;
features.edges.push([lo, hi]);
}
}
for w in run.windows(2) {
face_dsu.union(w[0].2, w[1].2);
}
start = end;
}
{
let mut seen: Vec<u32> = (0..faces.len() as u32).map(|f| face_dsu.find(f)).collect();
seen.sort_unstable();
seen.dedup();
report.components = seen.len() as u64;
let mut loops: Vec<u32> = (0..vertex_count as u32)
.filter(|v| on_boundary[*v as usize])
.map(|v| boundary_dsu.find(v))
.collect();
loops.sort_unstable();
loops.dedup();
report.boundary_loops = loops.len() as u64;
}
{
let mut vertex_faces: Vec<(u32, u32)> = Vec::with_capacity(faces.len() * 3);
for (fi, f) in faces.iter().enumerate() {
for &v in f {
vertex_faces.push((v, fi as u32));
}
}
vertex_faces.sort_unstable();
let mut wings: Vec<(u32, u32)> = Vec::new();
let mut link = Dsu::new(0);
let mut start = 0usize;
while start < vertex_faces.len() {
let v = vertex_faces[start].0;
let mut end = start + 1;
while end < vertex_faces.len() && vertex_faces[end].0 == v {
end += 1;
}
let degree = end - start;
wings.clear();
for (local, &(_, fi)) in vertex_faces[start..end].iter().enumerate() {
for &w in &faces[fi as usize] {
if w != v {
wings.push((w, local as u32));
}
}
}
wings.sort_unstable();
link.reset(degree);
let mut branching = false;
let mut i = 0usize;
while i < wings.len() {
let mut j = i + 1;
while j < wings.len() && wings[j].0 == wings[i].0 {
j += 1;
}
if j - i >= 3 {
branching = true;
}
for k in i + 1..j {
link.union(wings[i].1, wings[k].1);
}
i = j;
}
let mut roots: Vec<u32> = (0..degree as u32).map(|d| link.find(d)).collect();
roots.sort_unstable();
roots.dedup();
if branching || roots.len() > 1 {
report.non_manifold_vertices += 1;
features.vertices.push(v);
}
start = end;
}
}
{
let two_area_limit =
R::from_f64(2.0 * cfg.area_epsilon_rel * cfg.cell_size * cfg.cell_size);
let limit_sq = two_area_limit * two_area_limit;
for f in &faces {
let a = positions[f[0] as usize];
let b = positions[f[1] as usize];
let c = positions[f[2] as usize];
let u = [b[0] - a[0], b[1] - a[1], b[2] - a[2]];
let v = [c[0] - a[0], c[1] - a[1], c[2] - a[2]];
let cross = [
u[1] * v[2] - u[2] * v[1],
u[2] * v[0] - u[0] * v[2],
u[0] * v[1] - u[1] * v[0],
];
let len_sq = cross[0] * cross[0] + cross[1] * cross[1] + cross[2] * cross[2];
if len_sq <= limit_sq {
report.degenerate_triangles += 1;
}
}
}
{
let eps = R::from_f64(cfg.weld_epsilon);
let eps_sq = eps * eps;
let lattice = crate::weld::Lattice::new(positions, eps);
let key_of = |p: [R; 3]| lattice.key_of(p);
let mut cells: Vec<([i64; 3], u32)> = (0..vertex_count)
.map(|i| (key_of(positions[i]), i as u32))
.collect();
cells.sort_unstable();
{
let mut distinct = 0u64;
let mut previous: Option<[i64; 3]> = None;
for (k, _) in &cells {
if previous != Some(*k) {
distinct += 1;
previous = Some(*k);
}
}
report.weld_buckets = distinct;
}
for v in 0..vertex_count as u32 {
let p = positions[v as usize];
let base = key_of(p);
let mut found = false;
'probe: for dz in -1..=1i64 {
for dy in -1..=1i64 {
for dx in -1..=1i64 {
let key = [base[0] + dx, base[1] + dy, base[2] + dz];
let lo = cells.partition_point(|(k, _)| *k < key);
for &(k, u) in &cells[lo..] {
if k != key {
break;
}
if u >= v {
continue;
}
let q = positions[u as usize];
let d = [p[0] - q[0], p[1] - q[1], p[2] - q[2]];
if d[0] * d[0] + d[1] * d[1] + d[2] * d[2] <= eps_sq {
found = true;
break 'probe;
}
}
}
}
}
if found {
report.duplicate_vertices += 1;
}
}
}
report.euler_characteristic =
report.referenced_vertices as i64 - report.edges as i64 + report.faces as i64;
if report.components == 1
&& report.faces_skipped == 0
&& report.non_manifold_edges == 0
&& report.non_manifold_vertices == 0
&& report.inconsistently_oriented_edges == 0
{
let doubled = 2 - report.euler_characteristic - report.boundary_loops as i64;
if doubled % 2 == 0 {
report.genus = Some(doubled / 2);
}
}
(report, features)
}
#[must_use]
pub fn validate<R: Real>(mesh: &MeshBuffer<R>, cfg: &ValidateConfig) -> MeshReport {
let mut report = validate_indexed(&mesh.positions, &mesh.indices, cfg);
report.normal_count_mismatch = mesh.normals.len() != mesh.positions.len();
report
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub enum SurfaceGate {
Closed,
Manifold,
ClosedAllowingUnresolvedTopology,
}
#[cfg(test)]
mod tests;