linked_mesh.rs

use std::borrow::Cow;
use std::{f32::consts::PI, fmt::Debug, hash::Hash, marker::PhantomData, num::NonZeroU32};

use arrayvec::ArrayVec;
use bitvec::{bitarr, slice::BitSlice};
use fxhash::{FxBuildHasher, FxHashMap, FxHashSet};
use nalgebra::{Matrix3, Matrix4, Vector3};
use slotmap::{new_key_type, Key, SecondaryMap, SlotMap};
use smallvec::SmallVec;

use crate::attrs;
use crate::bvh::{Aabb, TriBvh};
#[cfg(not(target_arch = "wasm32"))]
use crate::models::{
  stanford_bunny::{STANFORD_BUNNY_INDICES, STANFORD_BUNNY_VERTICES},
  utah_teapot::{UTAH_TEAPOT_INDICES, UTAH_TEAPOT_VERTICES},
};
use crate::{
  slotmap_utils::{
    build_slotmap_from_iter, build_slotmap_from_iter_with_key, ekey_ix, fkey_ix,
    slotmap_insert_dense, slotmap_insert_dense_with_key, vkey, vkey_ix, vkey_ix_mut, vkey_version,
    LocalKeyData, LocalSlotMap,
  },
  OwnedIndexedMesh, OwnedIndexedMeshBuilder, OwnedMesh, Triangle,
};

pub type Vec3 = Vector3<f32>;
pub type Mat4 = Matrix4<f32>;

new_key_type! {
  pub struct VertexKey;
  pub struct FaceKey;
  pub struct EdgeKey;
}

#[derive(Clone, Debug, Default)]
pub struct Vertex {
  pub position: Vec3,
  pub edges: SmallVec<[EdgeKey; 9]>,
}

impl Vertex {
  pub fn new(position: Vec3) -> Self {
    Vertex {
      position,
      edges: SmallVec::new(),
    }
  }
}

#[derive(Clone, Debug)]
pub struct Face<T> {
  /// Counter-clockwise winding
  pub vertices: [VertexKey; 3],
  /// Unordered
  pub edges: [EdgeKey; 3],
  pub data: T,
}

static mut GRAPHVIZ_PRINT_INNER: fn(&str) -> () = |s| {
  println!("{}", s);
};

fn graphviz_print(s: &str) {
  unsafe {
    GRAPHVIZ_PRINT_INNER(s);
  }
}

pub fn set_graphviz_print(f: fn(&str) -> ()) {
  unsafe {
    GRAPHVIZ_PRINT_INNER = f;
  }
}

static mut DEBUG_PRINT_INNER: fn(&str) -> () = |s| {
  println!("{s}");
};

#[allow(dead_code)]
fn debug_print(s: &str) {
  unsafe {
    DEBUG_PRINT_INNER(s);
  }
}

pub fn set_debug_print(f: fn(&str) -> ()) {
  unsafe {
    DEBUG_PRINT_INNER = f;
  }
}

// /// Same as Vec::retain but it doesn't preserve order of elements and can possibly work a bit
// faster /// if only removing a few elements.
// fn swap_retain<T>(vec: &mut Vec<T>, mut f: impl FnMut(&mut T) -> bool) {
//   let mut i = 0;
//   while i < vec.len() {
//     if !f(&mut vec[i]) {
//       vec.swap_remove(i);
//     } else {
//       i += 1;
//     }
//   }
// }

/// Same as SmallVec::retain but it doesn't preserve order of elements and can possibly work a bit
/// faster if only removing a few elements.
fn swap_retain_sv<T: smallvec::Array>(
  vec: &mut SmallVec<T>,
  mut f: impl FnMut(&mut <T as smallvec::Array>::Item) -> bool,
) {
  let mut i = 0;
  while i < vec.len() {
    if !f(&mut vec[i]) {
      vec.swap_remove(i);
    } else {
      i += 1;
    }
  }
}

impl<T> Face<T> {
  pub fn vertex_positions(&self, verts: &SlotMap<VertexKey, Vertex>) -> [Vec3; 3] {
    [
      verts[self.vertices[0]].position,
      verts[self.vertices[1]].position,
      verts[self.vertices[2]].position,
    ]
  }

  pub fn iter_vtx_positions<'a>(
    &'a self,
    verts: &'a SlotMap<VertexKey, Vertex>,
  ) -> impl Iterator<Item = &'a Vec3> + 'a {
    self
      .vertices
      .iter()
      .map(move |&vtx_key| &verts[vtx_key].position)
  }

  pub fn normal(&self, verts: &SlotMap<VertexKey, Vertex>) -> Vector3<f32> {
    let [a, b, c] = self.vertex_positions(verts);
    (b - a).cross(&(c - a)).normalize()
  }

  pub fn area(&self, verts: &SlotMap<VertexKey, Vertex>) -> f32 {
    let [a, b, c] = self.vertex_positions(verts);
    0.5 * (b - a).cross(&(c - a)).magnitude()
  }

  pub fn center(&self, verts: &SlotMap<VertexKey, Vertex>) -> Vector3<f32> {
    let [a, b, c] = self.vertex_positions(verts);
    (a + b + c) / 3.
  }

  pub fn compute_angle_at_vertex(&self, vtx_ix: usize, verts: &SlotMap<VertexKey, Vertex>) -> f32 {
    let (target_vtx_key, b, c) = match vtx_ix {
      0 => (self.vertices[0], self.vertices[1], self.vertices[2]),
      1 => (self.vertices[1], self.vertices[2], self.vertices[0]),
      2 => (self.vertices[2], self.vertices[0], self.vertices[1]),
      _ => panic!("Vertex index {vtx_ix} out of bounds for triangle with 3 vertices"),
    };
    let [target_vtx_pos, b_vtx_pos, c_vtx_pos] = [
      verts[target_vtx_key].position,
      verts[b].position,
      verts[c].position,
    ];

    let edge_0 = b_vtx_pos - target_vtx_pos;
    let edge_1 = c_vtx_pos - target_vtx_pos;
    edge_0.angle(&edge_1)
  }

  pub fn compute_angle_at_vertex_key(
    &self,
    vtx_key: VertexKey,
    verts: &SlotMap<VertexKey, Vertex>,
  ) -> f32 {
    let vtx_ix = self
      .vertices
      .iter()
      .position(|&v| v == vtx_key)
      .unwrap_or_else(|| panic!("Vertex key {vtx_key:?} not found in face"));
    self.compute_angle_at_vertex(vtx_ix, verts)
  }

  pub fn to_triangle(&self, verts: &SlotMap<VertexKey, Vertex>) -> Triangle {
    Triangle {
      a: verts[self.vertices[0]].position,
      b: verts[self.vertices[1]].position,
      c: verts[self.vertices[2]].position,
    }
  }

  pub fn is_degenerate(&self, verts: &SlotMap<VertexKey, Vertex>) -> bool {
    self.to_triangle(verts).is_degenerate()
  }
}

#[derive(Clone, Debug)]
pub struct Edge {
  // Ordered such that the first vertex key is always less than the second
  pub vertices: [VertexKey; 2],
  pub faces: SmallVec<[FaceKey; 2]>,
  pub sharp: bool,
}

impl Edge {
  pub fn length(&self, verts: &SlotMap<VertexKey, Vertex>) -> f32 {
    let a = verts[self.vertices[0]].position;
    let b = verts[self.vertices[1]].position;
    (a - b).magnitude()
  }

  pub fn other_vtx(&self, vtx_key: VertexKey) -> VertexKey {
    if self.vertices[0] == vtx_key {
      self.vertices[1]
    } else if self.vertices[1] == vtx_key {
      self.vertices[0]
    } else {
      panic!("Vertex key {vtx_key:?} not found in edge");
    }
  }
}

fn distance(p0: Vec3, p1: Vec3) -> f32 {
  (p0 - p1).magnitude()
}

/// Fast `acos` approximation (Abramowitz & Stegun 4.4.45), max abs error ~6.7e-5 rad. Much cheaper
/// than libm `acos`, which has no hardware support on wasm.
#[cfg_attr(target_arch = "wasm32", allow(dead_code))]
#[inline]
fn fast_acos(x: f32) -> f32 {
  let neg = if x < 0. { 1.0f32 } else { 0.0 };
  let x = x.abs();
  let mut ret = -0.0187293;
  ret = ret * x + 0.074261;
  ret = ret * x - 0.2121144;
  ret = ret * x + 1.5707288;
  ret *= (1.0 - x).sqrt();
  ret -= 2. * neg * ret;
  neg * PI + ret
}

/// SIMD 4-lane `fast_acos`.
#[cfg(target_arch = "wasm32")]
#[inline]
fn fast_acos_x4(x: core::arch::wasm32::v128) -> core::arch::wasm32::v128 {
  use core::arch::wasm32::*;
  let neg = v128_and(f32x4_lt(x, f32x4_splat(0.)), f32x4_splat(1.0));
  let x = f32x4_abs(x);
  let mut ret = f32x4_splat(-0.0187293);
  ret = f32x4_add(f32x4_mul(ret, x), f32x4_splat(0.074261));
  ret = f32x4_add(f32x4_mul(ret, x), f32x4_splat(-0.2121144));
  ret = f32x4_add(f32x4_mul(ret, x), f32x4_splat(1.5707288));
  ret = f32x4_mul(ret, f32x4_sqrt(f32x4_sub(f32x4_splat(1.0), x)));
  ret = f32x4_sub(ret, f32x4_mul(f32x4_mul(f32x4_splat(2.0), neg), ret));
  f32x4_add(f32x4_mul(neg, f32x4_splat(PI)), ret)
}

#[inline]
fn batch_acos(cos0: &[f32; 4], cos1: &[f32; 4]) -> ([f32; 4], [f32; 4]) {
  #[cfg(target_arch = "wasm32")]
  {
    use core::arch::wasm32::*;
    let (mut a0, mut a1) = ([0f32; 4], [0f32; 4]);
    unsafe {
      v128_store(
        a0.as_mut_ptr().cast(),
        fast_acos_x4(v128_load(cos0.as_ptr().cast())),
      );
      v128_store(
        a1.as_mut_ptr().cast(),
        fast_acos_x4(v128_load(cos1.as_ptr().cast())),
      );
    }
    (a0, a1)
  }
  #[cfg(not(target_arch = "wasm32"))]
  {
    (cos0.map(fast_acos), cos1.map(fast_acos))
  }
}

/// Batches 4 faces' corner cosines so each `acos` resolves via one SIMD `fast_acos_x4`. The third
/// corner angle is derived from `a0 + a1 + a2 == PI`, saving a third `acos`.
#[derive(Default)]
struct AcosBatch {
  n: usize,
  key: [FaceKey; 4],
  normal: [Vec3; 4],
  cos0: [f32; 4],
  cos1: [f32; 4],
}

impl AcosBatch {
  #[inline]
  fn push(&mut self, key: FaceKey, normal: Vec3, cos0: f32, cos1: f32, data: &mut [FaceShading]) {
    let n = self.n;
    self.key[n] = key;
    self.normal[n] = normal;
    self.cos0[n] = cos0;
    self.cos1[n] = cos1;
    self.n = n + 1;
    if self.n == 4 {
      self.flush(data);
    }
  }

  #[inline]
  fn flush(&mut self, data: &mut [FaceShading]) {
    let (a0s, a1s) = batch_acos(&self.cos0, &self.cos1);
    for i in 0..self.n {
      let a2 = (PI - a0s[i] - a1s[i]).max(0.);
      data[fkey_ix(&self.key[i]) as usize] = FaceShading {
        normal: self.normal[i],
        angles: [a0s[i], a1s[i], a2],
      };
    }
    self.n = 0;
  }
}

/// Mesh representation that maintains topological information between vertices,
/// edges, and faces.
///
/// Similar to a half-edge data structure, but with a some differences such as
/// all faces being triangles and no intrinsic directionality to edges.
/// Instead, faces store their vertices directly in-order.
#[derive(Clone, Default)]
pub struct LinkedMesh<FaceData = ()> {
  pub vertices: SlotMap<VertexKey, Vertex>,
  pub faces: SlotMap<FaceKey, Face<FaceData>>,
  pub edges: SlotMap<EdgeKey, Edge>,
  pub transform: Option<Mat4>,
  /// Per-vertex normal used for shading/lighting. Sparse: absence = no normal.
  pub shading_normals: SecondaryMap<VertexKey, Vec3>,
  /// Per-vertex normal used for displacement mapping.
  pub displacement_normals: SecondaryMap<VertexKey, Vec3>,
  /// Per-edge normal used for displacement mapping.
  pub edge_displacement_normals: SecondaryMap<EdgeKey, Vec3>,
  /// Generic passive per-vertex attribute channels (uv, tangent, color, user-defined), keyed by
  /// channel name.
  pub vertex_channels: FxHashMap<String, Channel<VertexKey>>,
  /// Generic passive per-edge attribute channels.
  pub edge_channels: FxHashMap<String, Channel<EdgeKey>>,
  /// `mesh_flags` bitset controlling the render finalize pass. Default 0 = full default treatment.
  pub flags: u32,
}

/// Opt-in per-mesh flags. Whether to recompute shading normals is decided from ground truth (are
/// the `shading_normals` complete), but welding coincident vertices can't be inferred that way — a
/// mesh with attribute seams (UV cuts, duplicated boundary rings) has position-coincident verts
/// that must stay distinct. Ops that author such seams set `NO_WELD` to opt out of the
/// distance-weld.
pub mod mesh_flags {
  /// Skip `merge_vertices_by_distance` at render finalize.
  pub const NO_WELD: u32 = 1 << 0;
}

impl<FaceData> LinkedMesh<FaceData> {
  #[inline]
  pub fn has_flag(&self, flag: u32) -> bool {
    self.flags & flag != 0
  }

  pub fn shading_normal(&self, key: VertexKey) -> Option<Vec3> {
    self.shading_normals.get(key).copied()
  }

  pub fn displacement_normal(&self, key: VertexKey) -> Option<Vec3> {
    self.displacement_normals.get(key).copied()
  }

  pub fn edge_displacement_normal(&self, key: EdgeKey) -> Option<Vec3> {
    self.edge_displacement_normals.get(key).copied()
  }

  pub fn set_shading_normal(&mut self, key: VertexKey, normal: Option<Vec3>) {
    set_channel(&mut self.shading_normals, key, normal);
  }

  pub fn set_displacement_normal(&mut self, key: VertexKey, normal: Option<Vec3>) {
    set_channel(&mut self.displacement_normals, key, normal);
  }

  pub fn set_edge_displacement_normal(&mut self, key: EdgeKey, normal: Option<Vec3>) {
    set_channel(&mut self.edge_displacement_normals, key, normal);
  }

  /// Removes a vertex along with any per-vertex channel data keyed to it. All vertex removals
  /// should go through here so the attribute channels don't accumulate stale entries as the mesh
  /// is edited.
  pub fn remove_vertex(&mut self, key: VertexKey) -> Option<Vertex> {
    self.shading_normals.remove(key);
    self.displacement_normals.remove(key);
    for ch in self.vertex_channels.values_mut() {
      ch.store.remove(key);
    }
    self.vertices.remove(key)
  }

  /// Removes an edge along with any per-edge channel data keyed to it.
  pub fn remove_edge(&mut self, key: EdgeKey) -> Option<Edge> {
    self.edge_displacement_normals.remove(key);
    for ch in self.edge_channels.values_mut() {
      ch.store.remove(key);
    }
    self.edges.remove(key)
  }

  /// Inserts a fresh vertex, cloning every passive channel value from `src` (the seam-duplicate
  /// case). Named normals are the caller's responsibility — they're recomputed/assigned bespoke.
  pub fn add_vertex_cloned_from(&mut self, src: VertexKey, position: Vec3) -> VertexKey {
    let key = self.vertices.insert(Vertex::new(position));
    for ch in self.vertex_channels.values_mut() {
      ch.store.clone_slot(key, src);
    }
    key
  }

  /// Inserts a fresh vertex, blending every passive channel from the weighted `srcs` per the
  /// channel's interp rule (the split/cut case). Named normals stay the caller's responsibility.
  pub fn add_vertex_interpolated(
    &mut self,
    srcs: &[(VertexKey, f32)],
    position: Vec3,
  ) -> VertexKey {
    let key = self.vertices.insert(Vertex::new(position));
    for ch in self.vertex_channels.values_mut() {
      ch.store.interp_slot(key, srcs, ch.interp);
    }
    key
  }

  /// Clones `src` into a fresh vertex inheriting its position, normals, and passive channels, then
  /// repoints `faces` onto the clone. Introduces a UV/material seam at a shared vertex *without*
  /// recomputing or creasing the smooth normal — the caller typically then overrides a channel
  /// (e.g. the UV `V` of a closed-profile wrap) on the returned clone. `faces` must reference
  /// `src`.
  pub fn split_off_faces(&mut self, src: VertexKey, faces: &[FaceKey]) -> VertexKey
  where
    FaceData: Default,
  {
    let position = self.vertices[src].position;
    let clone = self.add_vertex_cloned_from(src, position);
    let shading_normal = self.shading_normal(src);
    self.set_shading_normal(clone, shading_normal);
    let displacement_normal = self.displacement_normal(src);
    self.set_displacement_normal(clone, displacement_normal);
    for &face_key in faces {
      self.replace_vertex_in_face(face_key, src, clone);
    }
    clone
  }
}

fn set_channel<K: Key>(channel: &mut SecondaryMap<K, Vec3>, key: K, value: Option<Vec3>) {
  match value {
    Some(v) => {
      channel.insert(key, v);
    }
    None => {
      channel.remove(key);
    }
  }
}

#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum Arity {
  Scalar,
  Vec2,
  Vec3,
  Vec4,
}

/// How a passive channel's value at a newly-created vertex/edge is derived from its sources.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum Interp {
  /// Weighted average of the sources (uv, color, …).
  Lerp,
  /// Weighted average, then renormalized — for unit-length data like tangents.
  LerpNormalize,
}

/// How a passive channel transforms when the surface is flipped/inverted (winding reversed).
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum FlipXform {
  Identity,
  Negate,
}

/// How a passive channel responds to a spatial (affine) transform baked into the mesh.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum SpatialXform {
  /// Not a spatial quantity (uv, color, …); left unchanged.
  Identity,
  /// A direction lying in the surface (tangent); its `xyz` lanes are multiplied by the transform's
  /// linear part.
  Direction,
}

/// Tight per-arity SoA storage for one passive attribute channel.
#[derive(Clone, Debug)]
pub enum ChannelStore<K: Key> {
  Scalar(SecondaryMap<K, f32>),
  Vec2(SecondaryMap<K, [f32; 2]>),
  Vec3(SecondaryMap<K, [f32; 3]>),
  Vec4(SecondaryMap<K, [f32; 4]>),
}

impl<K: Key> ChannelStore<K> {
  pub fn new(arity: Arity) -> Self {
    match arity {
      Arity::Scalar => ChannelStore::Scalar(SecondaryMap::new()),
      Arity::Vec2 => ChannelStore::Vec2(SecondaryMap::new()),
      Arity::Vec3 => ChannelStore::Vec3(SecondaryMap::new()),
      Arity::Vec4 => ChannelStore::Vec4(SecondaryMap::new()),
    }
  }

  pub fn arity(&self) -> usize {
    match self {
      ChannelStore::Scalar(_) => 1,
      ChannelStore::Vec2(_) => 2,
      ChannelStore::Vec3(_) => 3,
      ChannelStore::Vec4(_) => 4,
    }
  }

  /// Reads a slot as a fixed-width float vector (unused lanes zeroed) so the propagation math can
  /// stay arity-agnostic; `set` writes back only the lanes this arity actually stores.
  fn get(&self, key: K) -> Option<[f32; 4]> {
    match self {
      ChannelStore::Scalar(m) => m.get(key).map(|&v| [v, 0., 0., 0.]),
      ChannelStore::Vec2(m) => m.get(key).map(|&[a, b]| [a, b, 0., 0.]),
      ChannelStore::Vec3(m) => m.get(key).map(|&[a, b, c]| [a, b, c, 0.]),
      ChannelStore::Vec4(m) => m.get(key).copied(),
    }
  }

  fn set(&mut self, key: K, v: [f32; 4]) {
    match self {
      ChannelStore::Scalar(m) => {
        m.insert(key, v[0]);
      }
      ChannelStore::Vec2(m) => {
        m.insert(key, [v[0], v[1]]);
      }
      ChannelStore::Vec3(m) => {
        m.insert(key, [v[0], v[1], v[2]]);
      }
      ChannelStore::Vec4(m) => {
        m.insert(key, v);
      }
    }
  }

  fn remove(&mut self, key: K) {
    match self {
      ChannelStore::Scalar(m) => {
        m.remove(key);
      }
      ChannelStore::Vec2(m) => {
        m.remove(key);
      }
      ChannelStore::Vec3(m) => {
        m.remove(key);
      }
      ChannelStore::Vec4(m) => {
        m.remove(key);
      }
    }
  }

  fn clone_slot(&mut self, dst: K, src: K) {
    if let Some(v) = self.get(src) {
      self.set(dst, v);
    }
  }

  /// Weighted blend of `srcs`; `None` if any source lacks a value (matches the named-normal `zip`
  /// semantics).
  fn blend(&self, srcs: &[(K, f32)], interp: Interp) -> Option<[f32; 4]> {
    let arity = self.arity();
    let mut acc = [0f32; 4];
    for &(key, weight) in srcs {
      let v = self.get(key)?;
      for i in 0..4 {
        acc[i] += v[i] * weight;
      }
    }
    if interp == Interp::LerpNormalize && arity >= 2 {
      let norm = (0..arity).map(|i| acc[i] * acc[i]).sum::<f32>().sqrt();
      if norm > 1e-12 {
        for x in acc.iter_mut().take(arity) {
          *x /= norm;
        }
      }
    }
    Some(acc)
  }

  fn interp_slot(&mut self, dst: K, srcs: &[(K, f32)], interp: Interp) {
    if let Some(v) = self.blend(srcs, interp) {
      self.set(dst, v);
    }
  }

  fn flip_slot(&mut self, key: K, flip: FlipXform) {
    if flip == FlipXform::Negate {
      if let Some(mut v) = self.get(key) {
        for x in &mut v {
          *x = -*x;
        }
        self.set(key, v);
      }
    }
  }

  fn flip_all(&mut self, flip: FlipXform) {
    if flip != FlipXform::Negate {
      return;
    }
    match self {
      ChannelStore::Scalar(m) => m.values_mut().for_each(|v| *v = -*v),
      ChannelStore::Vec2(m) => m
        .values_mut()
        .for_each(|v| v.iter_mut().for_each(|x| *x = -*x)),
      ChannelStore::Vec3(m) => m
        .values_mut()
        .for_each(|v| v.iter_mut().for_each(|x| *x = -*x)),
      ChannelStore::Vec4(m) => m
        .values_mut()
        .for_each(|v| v.iter_mut().for_each(|x| *x = -*x)),
    }
  }

  /// Multiplies the `xyz` lanes of every value by `m`, renormalizing to preserve unit length. Used
  /// to carry direction channels (tangents) through a spatial transform; no-op for scalar/vec2.
  fn transform_directions(&mut self, m: &Matrix3<f32>) {
    let apply = |xyz: Vec3| {
      let out = m * xyz;
      let n = out.norm();
      if n > 1e-12 {
        out / n
      } else {
        out
      }
    };
    let write = |v: &mut [f32]| {
      let t = apply(Vec3::new(v[0], v[1], v[2]));
      v[0] = t.x;
      v[1] = t.y;
      v[2] = t.z;
    };
    match self {
      ChannelStore::Scalar(_) | ChannelStore::Vec2(_) => {}
      ChannelStore::Vec3(map) => map.values_mut().for_each(|v| write(v)),
      ChannelStore::Vec4(map) => map.values_mut().for_each(|v| write(v)),
    }
  }
}

/// A passive (non-normal) attribute channel: arity-tagged SoA storage plus the rules for how it
/// propagates through vertex/edge creation and surface flips.
#[derive(Clone, Debug)]
pub struct Channel<K: Key> {
  pub interp: Interp,
  pub flip: FlipXform,
  pub spatial: SpatialXform,
  pub store: ChannelStore<K>,
}

impl<K: Key> Channel<K> {
  pub fn new(arity: Arity, interp: Interp, flip: FlipXform, spatial: SpatialXform) -> Self {
    Channel {
      interp,
      flip,
      spatial,
      store: ChannelStore::new(arity),
    }
  }

  /// Carries a `Direction` channel through a spatial transform via its linear part; no-op
  /// otherwise.
  pub fn transform_directions(&mut self, m: &Matrix3<f32>) {
    if self.spatial == SpatialXform::Direction {
      self.store.transform_directions(m);
    }
  }

  /// Weighted blend of `srcs` per this channel's interp rule; `None` if any source lacks a value.
  pub fn blend(&self, srcs: &[(K, f32)]) -> Option<[f32; 4]> {
    self.store.blend(srcs, self.interp)
  }

  pub fn transform_direction_slot(&mut self, key: K, m: &Matrix3<f32>) {
    if self.spatial != SpatialXform::Direction {
      return;
    }
    if let Some(mut v) = self.store.get(key) {
      let out = m * Vec3::new(v[0], v[1], v[2]);
      let n = out.norm();
      let out = if n > 1e-12 { out / n } else { out };
      v[0] = out.x;
      v[1] = out.y;
      v[2] = out.z;
      self.store.set(key, v);
    }
  }

  /// Applies this channel's flip rule to one slot (surface winding reversed there).
  pub fn flip_slot(&mut self, key: K) {
    self.store.flip_slot(key, self.flip);
  }

  /// Writes `value` for `key`; lanes above this channel's arity are ignored, so callers pass
  /// `[x, y, 0., 0.]` for a `Vec2` channel, etc.
  pub fn set(&mut self, key: K, value: [f32; 4]) {
    self.store.set(key, value);
  }

  /// Reads `key` as a fixed-width vector; lanes above this channel's arity are zero.
  pub fn get(&self, key: K) -> Option<[f32; 4]> {
    self.store.get(key)
  }

  /// An empty channel sharing this one's arity, interp, and flip rules.
  pub fn empty_like(&self) -> Self {
    let store = match &self.store {
      ChannelStore::Scalar(_) => ChannelStore::Scalar(SecondaryMap::new()),
      ChannelStore::Vec2(_) => ChannelStore::Vec2(SecondaryMap::new()),
      ChannelStore::Vec3(_) => ChannelStore::Vec3(SecondaryMap::new()),
      ChannelStore::Vec4(_) => ChannelStore::Vec4(SecondaryMap::new()),
    };
    Channel {
      interp: self.interp,
      flip: self.flip,
      spatial: self.spatial,
      store,
    }
  }
}

impl<T> Debug for LinkedMesh<T> {
  fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
    write!(
      f,
      "LinkedMesh {{ vertices: {}, faces: {}, edges: {} }}",
      self.vertices.len(),
      self.faces.len(),
      self.edges.len()
    )
  }
}

fn sort_edge(v0: VertexKey, v1: VertexKey) -> [VertexKey; 2] {
  let v0_ix = vkey_ix(&v0);
  let v1_ix = vkey_ix(&v1);
  if v0_ix < v1_ix {
    [v0, v1]
  } else {
    [v1, v0]
  }
}

fn sort_edge_inplace<const DENSE: bool>(vtxs: &mut [VertexKey; 2]) {
  if DENSE {
    let [vtx0, vtx1] = vtxs;
    debug_assert_eq!(vkey_version(vtx0), 1);
    debug_assert_eq!(vkey_version(vtx1), 1);
    let v0_ix = vkey_ix_mut(vtx0);
    let v1_ix = vkey_ix_mut(vtx1);

    if *v0_ix > *v1_ix {
      std::mem::swap(v0_ix, v1_ix);
    }
  } else {
    let v0_ix = vkey_ix(&vtxs[0]);
    let v1_ix = vkey_ix(&vtxs[1]);
    if v0_ix > v1_ix {
      vtxs.swap(0, 1);
    }
  }
}

/// Per-face shading data precomputed once before fan-walking. Indexed by `fkey_ix`. A degenerate
/// face is flagged by an all-zero `normal` (a valid face normal is unit).
#[derive(Clone, Default)]
struct FaceShading {
  normal: Vec3,
  /// Corner angle at the owning face's i-th vertex; the normal's accumulation weight.
  angles: [f32; 3],
}

impl FaceShading {
  /// Angle-weighted normal at `center`; zero for degenerate faces (a no-op when accumulated).
  #[inline]
  fn weighted_at(&self, verts: &[VertexKey; 3], center: VertexKey) -> Vec3 {
    let corner = if verts[0] == center {
      0
    } else if verts[1] == center {
      1
    } else {
      2
    };
    self.normal * self.angles[corner]
  }
}

#[derive(Debug)]
struct NormalAcc {
  pub accumulated_normal: Vec3,
}

impl NormalAcc {
  pub fn get(&self) -> Option<Vec3> {
    if self.accumulated_normal == Vec3::zeros() {
      None
    } else {
      let normalized = self.accumulated_normal.normalize();
      if normalized.x.is_nan() || normalized.y.is_nan() || normalized.z.is_nan() {
        panic!("Non-zero but still NaN normalized: {self:?}");
      }
      Some(normalized)
    }
  }

  fn new() -> Self {
    Self {
      accumulated_normal: Vec3::zeros(),
    }
  }
}

struct SmoothFan {
  pub old_key: VertexKey,
  pub face_keys: SmallVec<[FaceKey; 8]>,
  pub normal: Vec3,
}

/// Determines the method used to compute displacement normals for new vertices
/// produced when splitting edges
#[repr(u32)]
#[derive(Clone, Copy)]
pub enum DisplacementNormalMethod {
  /// Displacement normals are generated by mixing the normals of the vertices
  /// of edges when splitting them
  Interpolate = 0,
  /// Displacement normals are assigned from the normal of the edge when
  /// splitting them
  EdgeNormal = 1,
}

impl TryFrom<u32> for DisplacementNormalMethod {
  type Error = ();

  fn try_from(value: u32) -> Result<Self, Self::Error> {
    match value {
      0 => Ok(DisplacementNormalMethod::Interpolate),
      1 => Ok(DisplacementNormalMethod::EdgeNormal),
      _ => Err(()),
    }
  }
}

#[derive(Debug, Clone)]
pub enum NonManifoldError {
  /// No faces in the mesh
  EmptyMesh,
  /// Edge is not connected to any faces
  LooseEdge { edge_key: EdgeKey },
  /// Vertex is either not connected to any edges or not used in any faces
  LooseVertex { vtx_key: VertexKey },
  /// Edge is connected to more than two faces
  NonManifoldEdge {
    edge_key: EdgeKey,
    face_count: usize,
  },
  /// Vertex has more than one fan of faces around it
  MultipleFans {
    vtx_key: VertexKey,
    incident_face_count: usize,
    visited_face_count: usize,
  },
  /// The fan of faces around a vertex is not closed
  NonClosedFan { vtx_key: VertexKey },
}

pub struct EdgeSplitPos {
  /// Number from 0. to 1. representing the position along the edge from the
  /// starting vertex to split at
  pub pos: f32,
  pub start_vtx_key: VertexKey,
}

impl EdgeSplitPos {
  fn get(&self, v0_key: VertexKey) -> f32 {
    if v0_key == self.start_vtx_key {
      self.pos
    } else {
      1. - self.pos
    }
  }

  pub fn middle() -> Self {
    Self {
      pos: 0.5,
      start_vtx_key: VertexKey::null(),
    }
  }
}

const EPSILON: f32 = 1e-5;

#[derive(Clone, Debug)]
pub struct Plane {
  pub normal: Vec3,
  pub w: f32,
}

#[derive(Debug, Clone, Copy, PartialEq)]
enum PolygonClass {
  Coplanar = 0,
  Front = 1,
  Back = 2,
  Spanning = 3,
}

impl From<u8> for PolygonClass {
  fn from(val: u8) -> Self {
    match val {
      0 => PolygonClass::Coplanar,
      1 => PolygonClass::Front,
      2 => PolygonClass::Back,
      3 => PolygonClass::Spanning,
      _ => panic!("Invalid PolygonClass value"),
    }
  }
}

impl std::ops::BitOr for PolygonClass {
  type Output = Self;

  fn bitor(self, rhs: Self) -> Self::Output {
    let out = self as u8 | rhs as u8;
    out.into()
  }
}

/// Generates vertices arranged in rings at different heights.
///
/// Creates `height_segments + 1` rings, each with `radial_segments` vertices.
/// - `radius_fn(h)`: returns the radius at normalized height `h` (0.0 to 1.0)
/// - `y_fn(h)`: returns the y-coordinate at normalized height `h` (0.0 to 1.0)
fn generate_radial_vertices(
  vertices: &mut Vec<Vec3>,
  radial_segments: usize,
  height_segments: usize,
  radius_fn: impl Fn(f32) -> f32,
  y_fn: impl Fn(f32) -> f32,
) {
  for h in 0..=height_segments {
    let h_norm = h as f32 / height_segments as f32;
    let y = y_fn(h_norm);
    let radius = radius_fn(h_norm);
    for r in 0..radial_segments {
      let theta = (r as f32 / radial_segments as f32) * 2. * PI;
      let x = radius * theta.cos();
      let z = radius * theta.sin();
      vertices.push(Vec3::new(x, y, z));
    }
  }
}

/// Generates quad (two triangle) indices connecting adjacent height rings.
///
/// Connects rings of vertices arranged in the pattern created by `generate_radial_vertices`.
fn generate_radial_side_indices(
  indices: &mut Vec<u32>,
  radial_segments: usize,
  height_segments: usize,
) {
  for h in 0..height_segments {
    for r in 0..radial_segments {
      let v0 = h * radial_segments + r;
      let v1 = v0 + radial_segments;
      let next_r = (r + 1) % radial_segments;
      let v0_next = h * radial_segments + next_r;
      let v1_next = v0_next + radial_segments;

      indices.push(v0 as u32);
      indices.push(v1 as u32);
      indices.push(v0_next as u32);
      indices.push(v1 as u32);
      indices.push(v1_next as u32);
      indices.push(v0_next as u32);
    }
  }
}

/// Generates cap (fan) indices from a center vertex to a ring of vertices.
///
/// - `center_ix`: index of the center vertex
/// - `radial_segments`: number of vertices in the ring
/// - `ring_offset`: starting index of the ring vertices
/// - `reverse_winding`: if true, reverses the winding order (for top caps)
fn generate_cap_indices(
  indices: &mut Vec<u32>,
  center_ix: u32,
  radial_segments: usize,
  ring_offset: usize,
  reverse_winding: bool,
) {
  for i in 0..radial_segments {
    let curr = (ring_offset + i) as u32;
    let next = (ring_offset + ((i + 1) % radial_segments)) as u32;
    if reverse_winding {
      indices.push(next);
      indices.push(curr);
      indices.push(center_ix);
    } else {
      indices.push(center_ix);
      indices.push(curr);
      indices.push(next);
    }
  }
}

/// Regular icosahedron with unit circumradius: three orthogonal golden rectangles aligned to the
/// coordinate planes (2-fold symmetric about each axis).  Shared by the icosphere, icosahedron, and
/// dodecahedron (built as its dual).
const ICOSAHEDRON_VERTICES: &[Vec3] = &[
  Vec3::new(-0.5257311, 0.8506508, 0.),
  Vec3::new(0.5257311, 0.8506508, 0.),
  Vec3::new(-0.5257311, -0.8506508, 0.),
  Vec3::new(0.5257311, -0.8506508, 0.),
  Vec3::new(0., -0.5257311, 0.8506508),
  Vec3::new(0., 0.5257311, 0.8506508),
  Vec3::new(0., -0.5257311, -0.8506508),
  Vec3::new(0., 0.5257311, -0.8506508),
  Vec3::new(0.8506508, 0., -0.5257311),
  Vec3::new(0.8506508, 0., 0.5257311),
  Vec3::new(-0.8506508, 0., -0.5257311),
  Vec3::new(-0.8506508, 0., 0.5257311),
];

const ICOSAHEDRON_FACES: &[[usize; 3]] = &[
  [0, 11, 5],
  [0, 5, 1],
  [0, 1, 7],
  [0, 7, 10],
  [0, 10, 11],
  [1, 5, 9],
  [5, 11, 4],
  [11, 10, 2],
  [10, 7, 6],
  [7, 1, 8],
  [3, 9, 4],
  [3, 4, 2],
  [3, 2, 6],
  [3, 6, 8],
  [3, 8, 9],
  [4, 9, 5],
  [2, 4, 11],
  [6, 2, 10],
  [8, 6, 7],
  [9, 8, 1],
];

/// Triangle indices wound so the face normal points away from the origin, for origin-centered
/// convex solids (where "outward" means the normal agrees with the face centroid direction).
fn outward_tri(verts: &[Vec3], a: usize, b: usize, c: usize) -> [u32; 3] {
  let normal = (verts[b] - verts[a]).cross(&(verts[c] - verts[a]));
  let centroid = verts[a] + verts[b] + verts[c];
  if normal.dot(&centroid) >= 0. {
    [a as u32, b as u32, c as u32]
  } else {
    [a as u32, c as u32, b as u32]
  }
}

impl<FaceData: Default> LinkedMesh<FaceData> {
  pub fn new(vertex_count: usize, face_count: usize, transform: Option<Mat4>) -> Self {
    Self {
      vertices: SlotMap::with_capacity_and_key(vertex_count),
      faces: SlotMap::with_capacity_and_key(face_count),
      edges: SlotMap::with_key(),
      transform,
      shading_normals: SecondaryMap::new(),
      displacement_normals: SecondaryMap::new(),
      edge_displacement_normals: SecondaryMap::new(),
      vertex_channels: FxHashMap::default(),
      edge_channels: FxHashMap::default(),
      flags: 0,
    }
  }

  pub fn iter_faces(&self) -> impl Iterator<Item = (FaceKey, &Face<FaceData>)> {
    self.faces.iter()
  }

  pub fn iter_edges(&self) -> impl Iterator<Item = (EdgeKey, &Edge)> {
    self.edges.iter()
  }

  pub fn iter_vertices(&self) -> impl Iterator<Item = (VertexKey, &Vertex)> {
    self.vertices.iter()
  }

  pub fn from_indexed_vertices(
    vertices: &[Vec3],
    indices: &[u32],
    normals: Option<&[Vec3]>,
    transform: Option<Mat4>,
  ) -> Self {
    if indices.len() % 3 != 0 {
      panic!(
        "Indices length must be a multiple of 3; got {}",
        indices.len()
      );
    }

    let vertices: SlotMap<VertexKey, Vertex> =
      build_slotmap_from_iter(vertices.iter().map(|&vtx_pos| Vertex::new(vtx_pos)));

    // `build_slotmap_from_iter` assigns dense keys `{ ix: i + 1, version: 1 }`.
    let shading_normals = match normals {
      Some(normals) => normals
        .iter()
        .enumerate()
        .map(|(i, &n)| (vkey(i as u32 + 1, 1), n))
        .collect(),
      None => SecondaryMap::new(),
    };

    let temp_faces: LocalSlotMap<FaceKey, Face<()>> = LocalSlotMap {
      slots: Vec::new(),
      free_head: 1,
      num_elems: 0,
      _k: PhantomData,
    };
    let mut mesh = Self {
      vertices,
      faces: unsafe { std::mem::transmute(temp_faces) },
      edges: SlotMap::with_key(),
      transform,
      shading_normals,
      displacement_normals: SecondaryMap::new(),
      edge_displacement_normals: SecondaryMap::new(),
      vertex_channels: FxHashMap::default(),
      edge_channels: FxHashMap::default(),
      flags: 0,
    };

    mesh.faces = build_slotmap_from_iter_with_key(
      indices.as_chunks::<3>().0.iter(),
      |[a_ix, b_ix, c_ix], face_key| {
        // we've inserted all the vertices in a fixed order without any removals, so the keys are
        // deterministic
        let a = vkey(a_ix + 1, 1);
        let b = vkey(b_ix + 1, 1);
        let c = vkey(c_ix + 1, 1);
        let vertices = [a, b, c];

        let edges = [(0, 1), (1, 2), (2, 0)];
        let face = Face {
          vertices,
          edges: edges.map(|(a, b)| {
            let edge_key = mesh.get_or_create_edge::<true>([vertices[a], vertices[b]]);

            let edge = unsafe { mesh.edges.get_unchecked_mut(edge_key) };
            edge.faces.push(face_key);

            unsafe {
              std::mem::transmute(LocalKeyData {
                idx: ekey_ix(&edge_key),
                version: NonZeroU32::new_unchecked(1),
              })
            }
          }),
          data: Default::default(),
        };

        face
      },
    );

    mesh
  }

  pub fn from_triangles(triangles: impl IntoIterator<Item = Triangle>) -> Self {
    let triangles = triangles.into_iter();
    let (min_size, max_size) = triangles.size_hint();
    let size = max_size.unwrap_or(min_size);
    let mut mesh = Self::new(size * 3, size, None);

    // TODO: optimize using fastpath slotmap methods
    for tri in triangles {
      // This might break mesh topology in some cases, but it saves us from dealing
      // with NaNs
      if tri.is_degenerate() {
        log::warn!("Skipping degenerate triangle: {tri:?}");
        continue;
      }

      let [a_key, b_key, c_key] = [
        mesh.vertices.insert(Vertex::new(tri.a)),
        mesh.vertices.insert(Vertex::new(tri.b)),
        mesh.vertices.insert(Vertex::new(tri.c)),
      ];

      mesh.add_face::<true>([a_key, b_key, c_key], Default::default());
    }

    mesh
  }

  /// Removes `v1` and updates all references to it to point to `v0` instead.
  pub fn merge_vertices(&mut self, v0_key: VertexKey, v1_key: VertexKey) {
    let removed_vtx = self.remove_vertex(v1_key).unwrap_or_else(|| {
      panic!(
        "Tried to merge vertex that doesn't exist; key={v1_key:?}. Was referenced by removed \
         vertex with key={v0_key:?}",
      )
    });

    for &edge_key in &removed_vtx.edges {
      let (pair_vtx_key, new_edge_vertices) = {
        let edge = &mut self.edges[edge_key];
        for &face_key in &edge.faces {
          let face = &mut self.faces[face_key];
          if face.vertices.contains(&v0_key) && face.vertices.contains(&v1_key) {
            let v0 = &self.vertices[v0_key];
            panic!(
              "Triangle {face_key:?} contains both vertices to merge: {v0_key:?}, {v1_key:?} with \
               positions {} and {}",
              v0.position, removed_vtx.position
            );
          }
          for vert_key in &mut face.vertices {
            if *vert_key == v1_key {
              *vert_key = v0_key;
            }
          }
        }

        let pair_vtx_key = if edge.vertices[0] == v1_key {
          edge.vertices[0] = v0_key;
          edge.vertices[1]
        } else if edge.vertices[1] == v1_key {
          edge.vertices[1] = v0_key;
          edge.vertices[0]
        } else {
          panic!("Edge {edge_key:?} doesn't contain vertex {v1_key:?} to merge",)
        };

        let v0 = &mut self.vertices[v0_key];
        v0.edges.push(edge_key);

        sort_edge_inplace::<false>(&mut edge.vertices);
        (pair_vtx_key, edge.vertices)
      };

      // We've updated the edge between (v1 and v_pair) to be between (v0 and v_pair)
      //
      // However, it's possible that there's already an edge between v0 and v_pair.
      // If that's the case, we merge this edge with that one and drop it from the
      // graph.
      let pair_vtx = &self.vertices[pair_vtx_key];
      let mut edge_key_to_merge_into = None;
      for &pair_edge_key in &pair_vtx.edges {
        if pair_edge_key == edge_key {
          continue;
        }

        let pair_edge = &self.edges[pair_edge_key];
        if pair_edge.vertices == new_edge_vertices {
          if edge_key_to_merge_into.is_some() {
            panic!(
              "Multiple edges found to merge into; \
               edge_key_to_merge_into={edge_key_to_merge_into:?}; pair_edge_key={pair_edge_key:?}",
            );
          }
          edge_key_to_merge_into = Some(pair_edge_key);
        }
      }

      if let Some(edge_key_to_merge_into) = edge_key_to_merge_into {
        let dropped_edge = self.remove_edge(edge_key).unwrap();
        let merged_into_edge = &mut self.edges[edge_key_to_merge_into];

        for &face_key in &dropped_edge.faces {
          let face = &mut self.faces[face_key];
          for o_edge_key in &mut face.edges {
            if *o_edge_key == edge_key {
              *o_edge_key = edge_key_to_merge_into;
            }
          }

          if !merged_into_edge.faces.contains(&face_key) {
            merged_into_edge.faces.push(face_key);
          }
        }

        for &vert_key in &dropped_edge.vertices {
          let vert = &mut self.vertices[vert_key];
          swap_retain_sv(&mut vert.edges, |&mut e| e != edge_key);
        }
      }
    }

    if cfg!(debug_assertions) {
      let mut uniq_edges = FxHashSet::default();
      for (_edge_key, edge) in &self.edges {
        if edge.vertices.contains(&v1_key) {
          panic!("Edge still contains removed vertex; edge={edge:?}; v0={v0_key:?}; v1={v1_key:?}",);
        }

        if !uniq_edges.insert(edge.vertices) {
          let dupe_edges = self
            .edges
            .iter()
            .filter(|(_, e)| e.vertices == edge.vertices)
            .collect::<Vec<_>>();
          let removed_vtx_edges_after = removed_vtx
            .edges
            .iter()
            .map(|&k| &self.edges[k])
            .collect::<Vec<_>>();
          graphviz_print(&self.debug());
          panic!(
            "Duplicate edge found after merging vertices removed_edge={edge:?};\n \
             dupe_edges={dupe_edges:?};\nv0={v0_key:?};\n v1={v1_key:?};\n \
             removed_vtx_edges_after={removed_vtx_edges_after:?}"
          );
        }
      }
    }
  }

  /// Whether two vertices agree on every passive channel (uv, tangent, etc.)
  fn vertex_channels_match(&self, a: VertexKey, b: VertexKey) -> bool {
    self
      .vertex_channels
      .values()
      .all(|ch| match (ch.get(a), ch.get(b)) {
        (Some(x), Some(y)) => x.iter().zip(&y).all(|(p, q)| (p - q).abs() < 1e-6),
        (None, None) => true,
        _ => false,
      })
  }

  // TODO: The spatial bucketing used here could maybe be replaced with the `rstar` crate
  // https://docs.rs/rstar/latest/rstar/
  pub fn merge_vertices_by_distance_cb(
    &mut self,
    max_distance: f32,
    mut cb: impl FnMut(FaceKey, FaceData),
  ) -> usize {
    // Spatial hashing. Raising `buckets_per_dim` would speed up clustered meshes (the per-bucket
    // scan is O(k^2) in occupancy) but changes which near-coincident verts greedily weld, altering
    // the merge result — a faster grouping-stable merge needs a union-find rewrite.
    let buckets_per_dim = 32usize;
    let mut buckets: FxHashMap<usize, Vec<_>> = FxHashMap::default();

    let (mut mins, mut maxs) = self.vertices.values().fold(
      (
        Vec3::new(f32::MAX, f32::MAX, f32::MAX),
        Vec3::new(f32::MIN, f32::MIN, f32::MIN),
      ),
      |(mins, maxs), vtx| (mins.inf(&vtx.position), maxs.sup(&vtx.position)),
    );
    // expand bounds and offset center slightly to allow more chances to avoid
    // searching neighboring buckets
    mins -= Vec3::new(0.3828593, 0.2859821, 0.18533925);
    maxs += Vec3::new(0.3828593, 0.2859821, 0.18533925);

    let x_range = maxs.x - mins.x;
    let y_range = maxs.y - mins.y;
    let z_range = maxs.z - mins.z;

    let get_bucket_ix = |x_ix: usize, y_ix: usize, z_ix: usize| -> usize {
      x_ix + y_ix * buckets_per_dim + z_ix * buckets_per_dim * buckets_per_dim
    };

    let get_bucket_coords = |pos: Vec3| -> ([usize; 3], bool) {
      let x_ix = ((pos.x - mins.x) / x_range * buckets_per_dim as f32).floor() as usize;
      let y_ix = ((pos.y - mins.y) / y_range * buckets_per_dim as f32).floor() as usize;
      let z_ix = ((pos.z - mins.z) / z_range * buckets_per_dim as f32).floor() as usize;

      // we can completely skip searching neighboring buckets if the distance between
      // this point and all of the bucket boundaries is greater than the max_distance
      let needs_neighbor_search = 'b: {
        let [bucket_mins_x, bucket_maxs_x] = [
          mins.x + x_ix as f32 / buckets_per_dim as f32 * x_range,
          mins.x + (x_ix + 1) as f32 / buckets_per_dim as f32 * x_range,
        ];
        if pos.x - bucket_mins_x < max_distance || bucket_maxs_x - pos.x < max_distance {
          break 'b true;
        }

        let [bucket_mins_y, bucket_maxs_y] = [
          mins.y + y_ix as f32 / buckets_per_dim as f32 * y_range,
          mins.y + (y_ix + 1) as f32 / buckets_per_dim as f32 * y_range,
        ];
        if pos.y - bucket_mins_y < max_distance || bucket_maxs_y - pos.y < max_distance {
          break 'b true;
        }

        let [bucket_mins_z, bucket_maxs_z] = [
          mins.z + z_ix as f32 / buckets_per_dim as f32 * z_range,
          mins.z + (z_ix + 1) as f32 / buckets_per_dim as f32 * z_range,
        ];
        if pos.z - bucket_mins_z < max_distance || bucket_maxs_z - pos.z < max_distance {
          break 'b true;
        }

        false
      };

      ([x_ix, y_ix, z_ix], needs_neighbor_search)
    };

    for (vtx_key, vtx) in &self.vertices {
      let (bucket_coords, _) = get_bucket_coords(vtx.position);
      let bucket_ix = get_bucket_ix(bucket_coords[0], bucket_coords[1], bucket_coords[2]);
      buckets
        .entry(bucket_ix)
        .or_insert_with(Vec::new)
        .push(vtx_key);
    }

    let mut removed_vert_keys = FxHashSet::default();
    let all_vert_keys = self.vertices.keys().collect::<Vec<_>>();

    let mut vertices_to_merge = Vec::new();
    for vtx_key in all_vert_keys {
      if removed_vert_keys.contains(&vtx_key) {
        continue;
      }

      let vtx = &self.vertices[vtx_key];
      let (bucket_coords, needs_neighbor_search) = get_bucket_coords(vtx.position);
      let bucket_ix = get_bucket_ix(bucket_coords[0], bucket_coords[1], bucket_coords[2]);

      let bucket = buckets.get_mut(&bucket_ix).unwrap();
      let mut o_vtx_ix = 0usize;
      while o_vtx_ix < bucket.len() {
        let o_vtx_key = bucket[o_vtx_ix];
        if o_vtx_key == vtx_key {
          o_vtx_ix += 1;
          continue;
        }

        let o_vtx = &self.vertices[o_vtx_key];
        // Don't weld a UV/attribute seam (coincident verts disagreeing on a passive channel).
        if distance(vtx.position, o_vtx.position) < max_distance
          && self.vertex_channels_match(vtx_key, o_vtx_key)
        {
          vertices_to_merge.push(o_vtx_key);
          bucket.swap_remove(o_vtx_ix);
        } else {
          o_vtx_ix += 1;
        }
      }

      if needs_neighbor_search {
        // Straddles a bucket boundary: scan the 26 neighbor cells (home cell already done above).
        for dx in -1isize..=1 {
          let x_ix = bucket_coords[0] as isize + dx;
          if x_ix < 0 || x_ix >= buckets_per_dim as isize {
            continue;
          }
          for dy in -1isize..=1 {
            let y_ix = bucket_coords[1] as isize + dy;
            if y_ix < 0 || y_ix >= buckets_per_dim as isize {
              continue;
            }
            for dz in -1isize..=1 {
              if dx == 0 && dy == 0 && dz == 0 {
                continue;
              }
              let z_ix = bucket_coords[2] as isize + dz;
              if z_ix < 0 || z_ix >= buckets_per_dim as isize {
                continue;
              }

              let neighbor_bucket_ix = get_bucket_ix(x_ix as usize, y_ix as usize, z_ix as usize);
              let Some(bucket) = buckets.get_mut(&neighbor_bucket_ix) else {
                continue;
              };

              let mut i = 0;
              while i < bucket.len() {
                let o_vtx_key = bucket[i];
                if o_vtx_key == vtx_key {
                  i += 1;
                  continue;
                }

                let o_vtx = &self.vertices[o_vtx_key];
                if distance(vtx.position, o_vtx.position) < max_distance
                  && self.vertex_channels_match(vtx_key, o_vtx_key)
                {
                  vertices_to_merge.push(o_vtx_key);
                  bucket.swap_remove(i);
                } else {
                  i += 1;
                }
              }
            }
          }
        }
      }

      for o_vtx_key in vertices_to_merge.drain(..) {
        // for every face that contains both vertices to merge, we remove it
        let mut face_keys_to_remove = Vec::new();
        for &edge_key in &self.vertices[o_vtx_key].edges {
          let edge = &self.edges[edge_key];
          for &face_key in &edge.faces {
            let face = &self.faces[face_key];
            if face.vertices.contains(&vtx_key) {
              face_keys_to_remove.push(face_key);
            }
          }
        }
        if !face_keys_to_remove.is_empty() {
          for face_key in face_keys_to_remove {
            // face might have been removed already
            if self.faces.contains_key(face_key) {
              let face_data = self.remove_face(face_key);
              cb(face_key, face_data);
            }
          }
        }

        self.merge_vertices(vtx_key, o_vtx_key);
        removed_vert_keys.insert(o_vtx_key);
      }
    }

    removed_vert_keys.len()
  }

  pub fn merge_vertices_by_distance(&mut self, max_distance: f32) -> usize {
    self.merge_vertices_by_distance_cb(max_distance, |_, _| ())
  }

  pub fn debug(&self) -> String {
    fn format_vtx(key: VertexKey, _vtx: &Vertex) -> String {
      format!("{key:?}",)
    }

    fn format_edge(key: EdgeKey, edge: &Edge) -> String {
      format!("{key:?} {:?} -> {:?}", edge.vertices[0], edge.vertices[1])
    }

    fn format_face<FaceData>(key: FaceKey, face: &Face<FaceData>) -> String {
      format!(
        "{key:?} {:?} -> {:?} -> {:?}",
        face.vertices[0], face.vertices[1], face.vertices[2]
      )
    }

    let mut connections = Vec::new();
    for (vtx_key, vtx) in self.iter_vertices() {
      let vtx_name = format_vtx(vtx_key, vtx);

      for &edge_key in &vtx.edges {
        let edge = &self.edges.get(edge_key).unwrap_or_else(|| {
          panic!(
            "Tried to get edge that doesn't exist; key={edge_key:?}. Was referenced by vertex \
             with key={vtx_key:?}",
          )
        });
        let edge = format_edge(edge_key, edge);
        connections.push((vtx_name.clone(), edge));
      }
    }

    for (edge_key, edge) in self.iter_edges() {
      let edge_name = format_edge(edge_key, edge);

      for &face_key in &edge.faces {
        let face = &self.faces[face_key];
        let face_name = format_face(face_key, face);
        connections.push((edge_name.clone(), face_name));
      }

      for &vtx_key in &edge.vertices {
        let vtx_name = self
          .vertices
          .get(vtx_key)
          .map(|vtx| format_vtx(vtx_key, vtx))
          .unwrap_or_else(|| {
            panic!(
              "Tried to get vertex that doesn't exist; key={vtx_key:?}. Was referenced by edge \
               with key={edge_key:?}"
            )
          });
        connections.push((edge_name.clone(), vtx_name));
      }
    }

    for (face_key, face) in self.iter_faces() {
      let face_name = format_face(face_key, face);

      for &edge_key in &face.edges {
        let edge = &self.edges[edge_key];
        let edge_name = format_edge(edge_key, edge);
        connections.push((face_name.clone(), edge_name));
      }

      for &vtx_key in &face.vertices {
        let vtx = &self.vertices[vtx_key];
        let vtx_name = format_vtx(vtx_key, vtx);
        connections.push((face_name.clone(), vtx_name));
      }
    }

    connections
      .into_iter()
      .map(|(a, b)| format!("{a}::{b}"))
      .collect::<Vec<_>>()
      .join("\n")
  }

  pub fn from_raw_indexed(
    vertices: &[f32],
    indices: &[u32],
    normals: Option<&[f32]>,
    transform: Option<Mat4>,
  ) -> Self {
    assert!(
      vertices.len() % 3 == 0,
      "Vertices length must be a multiple of 3; got {}",
      vertices.len()
    );
    let vertices =
      unsafe { std::slice::from_raw_parts(vertices.as_ptr() as *const Vec3, vertices.len() / 3) };
    let normals = if let Some(normals) = normals {
      if normals.is_empty() {
        None
      } else {
        assert!(
          normals.len() % 3 == 0,
          "Normals length must be a multiple of 3; got {}",
          normals.len()
        );
        Some(unsafe {
          std::slice::from_raw_parts(normals.as_ptr() as *const Vec3, normals.len() / 3)
        })
      }
    } else {
      None
    };

    Self::from_indexed_vertices(vertices, indices, normals, transform)
  }

  fn replace_vertex_in_face(
    &mut self,
    face_key: FaceKey,
    old_vtx_key: VertexKey,
    new_vtx_key: VertexKey,
  ) {
    let (
      edge_indices_to_alter,
      old_edge_keys,
      pair_vtx_keys,
      edge_displacement_normals,
      edge_sharpness,
    ) = {
      let face = &mut self.faces[face_key];
      let mut edge_indices_to_alter: [usize; 2] = [usize::MAX, usize::MAX];
      let mut pair_vtx_keys: [VertexKey; 2] = [VertexKey::null(), VertexKey::null()];
      let mut edge_displacement_normals: [Option<Vec3>; 2] = [None, None];
      let mut edge_sharpness: [bool; 2] = [false, false];
      let mut old_edge_keys: [EdgeKey; 2] = [EdgeKey::null(), EdgeKey::null()];

      for (edge_ix, &edge_key) in face.edges.iter().enumerate() {
        let edge = &self.edges[edge_key];
        if edge.vertices.contains(&old_vtx_key) {
          if edge_indices_to_alter[0] == usize::MAX {
            edge_indices_to_alter[0] = edge_ix;
            old_edge_keys[0] = edge_key;
            pair_vtx_keys[0] = if edge.vertices[0] == old_vtx_key {
              edge.vertices[1]
            } else {
              edge.vertices[0]
            };
            edge_displacement_normals[0] = self.edge_displacement_normals.get(edge_key).copied();
            edge_sharpness[0] = edge.sharp;
          } else {
            edge_indices_to_alter[1] = edge_ix;
            old_edge_keys[1] = edge_key;
            pair_vtx_keys[1] = if edge.vertices[0] == old_vtx_key {
              edge.vertices[1]
            } else {
              edge.vertices[0]
            };
            edge_displacement_normals[1] = self.edge_displacement_normals.get(edge_key).copied();
            edge_sharpness[1] = edge.sharp;
            break;
          }
        }
      }

      (
        edge_indices_to_alter,
        old_edge_keys,
        pair_vtx_keys,
        edge_displacement_normals,
        edge_sharpness,
      )
    };

    // remove this face from the old edges and delete them if they no longer
    // reference any faces
    for old_edge_key in old_edge_keys {
      let edge = &mut self.edges[old_edge_key];
      swap_retain_sv(&mut edge.faces, |&mut f| f != face_key);

      if edge.faces.is_empty() {
        for vtx_key in edge.vertices {
          let vtx = &mut self.vertices[vtx_key];
          swap_retain_sv(&mut vtx.edges, |&mut e| e != old_edge_key);
        }
        self.remove_edge(old_edge_key);
      }
    }

    let new_edge_key_0 = self.get_or_create_edge::<false>([new_vtx_key, pair_vtx_keys[0]]);
    let new_edge_key_1 = self.get_or_create_edge::<false>([new_vtx_key, pair_vtx_keys[1]]);
    self.edges[new_edge_key_0].faces.push(face_key);
    self.edges[new_edge_key_0].sharp = edge_sharpness[0];
    self.set_edge_displacement_normal(new_edge_key_0, edge_displacement_normals[0]);
    self.edges[new_edge_key_1].faces.push(face_key);
    self.edges[new_edge_key_1].sharp = edge_sharpness[1];
    self.set_edge_displacement_normal(new_edge_key_1, edge_displacement_normals[1]);

    let face = &mut self.faces[face_key];

    // update the edge keys
    face.edges[edge_indices_to_alter[0]] = new_edge_key_0;
    face.edges[edge_indices_to_alter[1]] = new_edge_key_1;

    // actually update the vertex key in the face
    let vtx_key_ix_to_alter = face
      .vertices
      .iter()
      .position(|&v| v == old_vtx_key)
      .unwrap();
    face.vertices[vtx_key_ix_to_alter] = new_vtx_key;
  }

  /// This check assumes that the mesh consists of a single connected component.  If there are
  /// islands or disconnected parts, this may produce incorrect results.
  ///
  /// If `TWO_MANIFOLD` is `true`, additionally enforces that every edge is shared by *exactly* two
  /// faces, meaning that the mesh is watertight and forms a continuous surface.
  ///
  /// The test is entirely topological; positions, normals, triangle winding order, etc. are not
  /// checked.
  pub fn check_is_manifold<const TWO_MANIFOLD: bool>(&self) -> Result<(), NonManifoldError> {
    for (edge_key, edge) in self.edges.iter() {
      let face_count = edge.faces.len();
      if TWO_MANIFOLD {
        if face_count != 2 {
          return Err(NonManifoldError::NonManifoldEdge {
            edge_key,
            face_count,
          });
        }
      } else {
        if face_count == 0 {
          log::error!("Found edge with no faces: {edge:?}");
          return Err(NonManifoldError::LooseEdge { edge_key });
        } else if face_count > 2 {
          return Err(NonManifoldError::NonManifoldEdge {
            edge_key,
            face_count,
          });
        }
      }
    }

    if self.faces.is_empty() {
      return Err(NonManifoldError::EmptyMesh);
    }

    // For each vertex, walk the fan of incident faces
    let mut incident_faces = Vec::new();
    let mut visited_faces = FxHashSet::default();
    for (vtx_key, vtx) in &self.vertices {
      if vtx.edges.is_empty() {
        return Err(NonManifoldError::LooseVertex { vtx_key });
      }

      incident_faces.clear();
      for &edge_key in &vtx.edges {
        incident_faces.extend(self.edges[edge_key].faces.iter().copied());
      }
      incident_faces.sort_unstable();
      incident_faces.dedup();
      let Some(&start_face_key) = incident_faces.first() else {
        return Err(NonManifoldError::LooseVertex { vtx_key });
      };
      let start_face = &self.faces[start_face_key];

      visited_faces.clear();

      // The two edges of the start face incident to this vertex. Each triangle touches a vertex
      // via exactly two edges; walking outward through one sends the traversal in one direction of
      // the fan. A boundary vertex's fan is an open path whose lowest-index (start) face may sit in
      // its middle, so we must try both directions to cover it.
      let start_vtx_edges: Vec<_> = start_face
        .edges
        .iter()
        .copied()
        .filter(|&e| self.edges[e].vertices.contains(&vtx_key))
        .collect();

      // Walk the fan from `start_face_key`, entering via `entry_edge` (so we leave through the
      // *other* vertex-edge), until a border edge or a full loop back to a visited face. Returns
      // whether the fan closed on itself; faces are unioned into `visited_faces`.
      let walk_fan = |entry_edge, visited_faces: &mut FxHashSet<_>| -> bool {
        let mut cur_face_key = start_face_key;
        let mut cur_edge_key = entry_edge;
        loop {
          visited_faces.insert(cur_face_key);
          let face = &self.faces[cur_face_key];
          let Some(next_edge_key) = face
            .edges
            .iter()
            .copied()
            .find(|&e| self.edges[e].vertices.contains(&vtx_key) && e != cur_edge_key)
          else {
            return false; // reached a border edge
          };
          let Some(&next_face_key) = self.edges[next_edge_key].faces.iter().find(|&&face_key| {
            face_key != cur_face_key && self.faces[face_key].vertices.contains(&vtx_key)
          }) else {
            return false; // reached a border edge
          };
          if visited_faces.contains(&next_face_key) {
            return true; // closed fan
          }
          cur_face_key = next_face_key;
          cur_edge_key = next_edge_key;
        }
      };

      let mut full_fan = walk_fan(start_vtx_edges[0], &mut visited_faces);
      if !full_fan {
        if let Some(&other_edge) = start_vtx_edges.get(1) {
          full_fan = walk_fan(other_edge, &mut visited_faces) || full_fan;
        }
      }

      if visited_faces.len() != incident_faces.len() {
        return Err(NonManifoldError::MultipleFans {
          vtx_key,
          incident_face_count: incident_faces.len(),
          visited_face_count: visited_faces.len(),
        });
      }

      if TWO_MANIFOLD {
        if !full_fan {
          return Err(NonManifoldError::NonClosedFan { vtx_key });
        }
      }
    }

    Ok(())
  }

  /// Returns `true` when the mesh is a single connected component.  This check assumes that the
  /// mesh consists of a single connected component.  If there are islands or disconnected parts,
  /// this may produce incorrect results.
  ///
  /// If `TWO_MANIFOLD` is `true`, additionally enforces that every edge is shared by *exactly* two
  /// faces, meaning that the mesh is watertight and forms a continuous surface.
  ///
  /// The test is entirely topological; positions, normals, triangle winding order, etc. are not
  /// checked.
  pub fn check_is_manifold_dynamic(&self, two_manifold: bool) -> Result<(), NonManifoldError> {
    for (edge_key, edge) in self.edges.iter() {
      let face_count = edge.faces.len();
      if two_manifold {
        if face_count != 2 {
          return Err(NonManifoldError::NonManifoldEdge {
            edge_key,
            face_count,
          });
        }
      } else {
        if face_count == 0 {
          log::error!("Found edge with no faces: {edge:?}");
          return Err(NonManifoldError::LooseEdge { edge_key });
        } else if face_count > 2 {
          return Err(NonManifoldError::NonManifoldEdge {
            edge_key,
            face_count,
          });
        }
      }
    }

    if self.faces.is_empty() {
      return Err(NonManifoldError::EmptyMesh);
    }

    // For each vertex, walk the fan of incident faces
    let mut incident_faces = Vec::new();
    let mut visited_faces = FxHashSet::default();
    let mut visited_edges = FxHashSet::default();
    for (vtx_key, vtx) in &self.vertices {
      if vtx.edges.is_empty() {
        return Err(NonManifoldError::LooseVertex { vtx_key });
      }

      incident_faces.clear();
      for &edge_key in &vtx.edges {
        incident_faces.extend(self.edges[edge_key].faces.iter().copied());
      }
      incident_faces.sort_unstable();
      incident_faces.dedup();
      let Some(&start_face_key) = incident_faces.first() else {
        return Err(NonManifoldError::LooseVertex { vtx_key });
      };
      let start_face = &self.faces[start_face_key];

      visited_faces.clear();
      visited_edges.clear();

      // Find an edge of the face that contains this vertex
      let start_edge_key = start_face
        .edges
        .iter()
        .find(|&&e| self.edges[e].vertices.contains(&vtx_key))
        .copied()
        .unwrap();
      let mut cur_face_key = start_face_key;
      let mut cur_edge_key = start_edge_key;
      let mut full_fan = false;
      loop {
        visited_faces.insert(cur_face_key);
        visited_edges.insert(cur_edge_key);
        // choose the next edge in the current face that contains this vertex and is not the
        // previous edge
        let face = &self.faces[cur_face_key];
        let next_edge_key = face
          .edges
          .iter()
          .find(|&&edge_key| {
            self.edges[edge_key].vertices.contains(&vtx_key) && edge_key != cur_edge_key
          })
          .copied();
        if next_edge_key.is_none() {
          break;
        }

        let next_edge_key = next_edge_key.unwrap();
        // Find the other face (besides cur_face_key) that shares this edge and contains the vertex
        let edge = &self.edges[next_edge_key];
        let Some(&next_face_key) = edge.faces.iter().find(|&&face_key| {
          face_key != cur_face_key && self.faces[face_key].vertices.contains(&vtx_key)
        }) else {
          // reached a border edge
          break;
        };

        if visited_faces.contains(&next_face_key) {
          full_fan = true;
          break;
        }
        cur_face_key = next_face_key;
        cur_edge_key = next_edge_key;
      }

      if visited_faces.len() != incident_faces.len() {
        return Err(NonManifoldError::MultipleFans {
          vtx_key,
          incident_face_count: incident_faces.len(),
          visited_face_count: visited_faces.len(),
        });
      }

      if two_manifold {
        if !full_fan {
          return Err(NonManifoldError::NonClosedFan { vtx_key });
        }
      }
    }

    Ok(())
  }

  pub fn compute_edge_displacement_normals(&mut self) {
    let edge_keys: Vec<EdgeKey> = self.edges.keys().collect();
    for edge_key in edge_keys {
      let edge_displacement_normal = self.edges[edge_key]
        .faces
        .iter()
        .map(|&face_key| {
          let face = &self.faces[face_key];
          face.normal(&self.vertices)
        })
        .sum::<Vec3>()
        .normalize();
      self
        .edge_displacement_normals
        .insert(edge_key, edge_displacement_normal);
    }
  }

  pub fn mark_edge_sharpness(&mut self, sharp_edge_threshold_rads: f32) {
    for edge in self.edges.values_mut() {
      // If edge has explicitly been marked as sharp, don't change it
      if edge.sharp {
        continue;
      }

      // Border edges as well as edges belonging to more than 3 faces are
      // automatically sharp
      if edge.faces.len() != 2 {
        edge.sharp = true;
        continue;
      }

      let [face0, face1] = [edge.faces[0], edge.faces[1]];
      let [normal0, normal1] = [
        self.faces[face0].normal(&self.vertices),
        self.faces[face1].normal(&self.vertices),
      ];
      let angle = normal0.angle(&normal1);
      edge.sharp = angle > sharp_edge_threshold_rads;
    }
  }

  /// Precomputes per face the normalized normal and three corner angles the fan walk needs
  fn precompute_face_shading<const COMPUTE_INCIDENT_COUNTS: bool>(
    &self,
  ) -> (Vec<FaceShading>, Vec<u32>) {
    let slot_count = self
      .faces
      .keys()
      .map(|k| fkey_ix(&k))
      .max()
      .map_or(0, |m| m as usize + 1);
    let mut data = vec![FaceShading::default(); slot_count];

    let vtx_slot_count = self
      .vertices
      .keys()
      .map(|k| vkey_ix(&k))
      .max()
      .map_or(0, |m| m as usize + 1);
    let mut incident_counts = if COMPUTE_INCIDENT_COUNTS {
      vec![0u32; vtx_slot_count]
    } else {
      Vec::new()
    };

    let mut batch = AcosBatch::default();
    for (face_key, face) in self.faces.iter() {
      if COMPUTE_INCIDENT_COUNTS {
        for &v in &face.vertices {
          unsafe { *incident_counts.get_unchecked_mut(vkey_ix(&v) as usize) += 1 };
        }
      }
      let [a, b, c] = unsafe {
        [
          self.vertices.get_unchecked(face.vertices[0]).position,
          self.vertices.get_unchecked(face.vertices[1]).position,
          self.vertices.get_unchecked(face.vertices[2]).position,
        ]
      };
      let ab = b - a;
      let ac = c - a;
      let cross = ab.cross(&ac);
      let mag = cross.magnitude();
      if !(mag > 0.) {
        // degenerate; leave the `Default` zero-normal sentinel
        continue;
      }
      let bc = c - b;
      let l_ab2 = ab.dot(&ab);
      let l_ac2 = ac.dot(&ac);
      let l_bc2 = bc.dot(&bc);
      let cos0 = (ab.dot(&ac) / (l_ab2 * l_ac2).sqrt()).clamp(-1., 1.);
      let cos1 = (-ab.dot(&bc) / (l_ab2 * l_bc2).sqrt()).clamp(-1., 1.);
      batch.push(face_key, cross / mag, cos0, cos1, &mut data);
    }
    batch.flush(&mut data);

    (data, incident_counts)
  }

  /// Walks one direction of a smooth fan from `(cur_edge, cur_face)` until a sharp/border/visited
  /// edge, accumulating each face's angle-weighted normal and recording visited edges/faces.
  /// `cur_edge_ix` is `cur_edge_key`'s index in `vtx.edges`, carried so the per-visit scan is gone.
  /// All keys come from the frozen-until-split topology, so the unchecked slotmap access is sound.
  #[inline]
  fn walk_fan_from(
    &self,
    vtx: &Vertex,
    vtx_key: VertexKey,
    face_data: &[FaceShading],
    mut cur_edge_key: EdgeKey,
    mut cur_edge_ix: usize,
    mut cur_face_key: FaceKey,
    visited_edges: &mut BitSlice,
    visited_faces: &mut SmallVec<[FaceKey; 16]>,
    smooth_fan_faces: &mut SmallVec<[FaceKey; 16]>,
    fan_normal_acc: &mut NormalAcc,
    vtx_normal_acc: &mut NormalAcc,
  ) {
    loop {
      smooth_fan_faces.push(cur_face_key);
      visited_faces.push(cur_face_key);
      let face = unsafe { self.faces.get_unchecked(cur_face_key) };
      let weighted_normal = unsafe { face_data.get_unchecked(fkey_ix(&cur_face_key) as usize) }
        .weighted_at(&face.vertices, vtx_key);
      fan_normal_acc.accumulated_normal += weighted_normal;
      vtx_normal_acc.accumulated_normal += weighted_normal;
      visited_edges.set(cur_edge_ix, true);

      // Next edge: cur_face's other spoke at `vtx`; the scan that finds it also yields its index.
      let next = face
        .edges
        .iter()
        .filter(|&&edge_key| edge_key != cur_edge_key)
        .find_map(|&edge_key| {
          vtx
            .edges
            .iter()
            .position(|&e| e == edge_key)
            .map(|ix| (edge_key, ix))
        });
      let Some((next_edge_key, next_edge_ix)) = next else {
        break;
      };
      let next_edge = unsafe { self.edges.get_unchecked(next_edge_key) };
      if visited_edges[next_edge_ix] {
        break;
      }

      let next_face_key = next_edge
        .faces
        .iter()
        .find(|&&face_key| face_key != cur_face_key && !visited_faces.contains(&face_key));
      let Some(&next_face_key) = next_face_key else {
        // Border edge
        visited_edges.set(next_edge_ix, true);
        break;
      };

      if next_edge.sharp {
        break;
      }

      cur_edge_key = next_edge_key;
      cur_edge_ix = next_edge_ix;
      cur_face_key = next_face_key;
    }
  }

  /// Accumulates faces that are part of a single smooth fan around `vtx_key`
  /// into `smooth_fan_faces`.
  ///
  /// Returns the normal of the smooth fan on `vtx_key`, weighted by the angle
  /// between the each face's edges that meet at `vtx_key`
  fn walk_one_smooth_fan(
    &mut self,
    vtx_key: VertexKey,
    face_data: &[FaceShading],
    visited_edges: &mut BitSlice,
    visited_faces: &mut SmallVec<[FaceKey; 16]>,
    smooth_fan_faces: &mut SmallVec<[FaceKey; 16]>,
    vtx_normal_acc: &mut NormalAcc,
  ) -> Option<Vec3> {
    let vtx = unsafe { self.vertices.get_unchecked(vtx_key) };

    let mut fan_normal_acc = NormalAcc::new();

    let start_edge_ix = visited_edges.iter().position(|visited| !*visited);
    let Some(start_edge_ix) = start_edge_ix else {
      // We've visited all edges
      return fan_normal_acc.get();
    };

    let start_edge_key = vtx.edges[start_edge_ix];
    let start_face_key = unsafe { self.edges.get_unchecked(start_edge_key) }
      .faces
      .iter()
      .find(|&&face_key| {
        if visited_faces.contains(&face_key) {
          return false;
        }

        let face = unsafe { self.faces.get_unchecked(face_key) };
        // Find the edges that include the vertex
        let [mut edge_key_0, mut edge_key_1] = [EdgeKey::null(), EdgeKey::null()];
        for &edge_key in &face.edges {
          if !unsafe { self.edges.get_unchecked(edge_key) }
            .vertices
            .contains(&vtx_key)
          {
            continue;
          }
          if edge_key_0.is_null() {
            edge_key_0 = edge_key
          } else {
            edge_key_1 = edge_key;
            break;
          }
        }

        // One edge must be not visited in order for us to walk from it
        [edge_key_0, edge_key_1].into_iter().any(|edge_key| {
          // TODO: why is this necessary
          let Some(edge_key_ix) = vtx.edges.iter().position(|&e| e == edge_key) else {
            return false;
          };
          !visited_edges[edge_key_ix]
        })
      });
    let Some(&start_face_key) = start_face_key else {
      return fan_normal_acc.get();
    };

    // If the starting edge is smooth, we have to walk the other way once we hit a
    // sharp edge.
    let needs_walk_the_other_way = !unsafe { self.edges.get_unchecked(start_edge_key) }.sharp;

    self.walk_fan_from(
      vtx,
      vtx_key,
      face_data,
      start_edge_key,
      start_edge_ix,
      start_face_key,
      visited_edges,
      visited_faces,
      smooth_fan_faces,
      &mut fan_normal_acc,
      vtx_normal_acc,
    );

    if !needs_walk_the_other_way {
      return fan_normal_acc.get();
    }

    let first_other_way_face_key = unsafe { self.edges.get_unchecked(start_edge_key) }
      .faces
      .iter()
      .find(|&&face_key| {
        let face = unsafe { self.faces.get_unchecked(face_key) };
        face.vertices.contains(&vtx_key) && !visited_faces.contains(&face_key)
      });
    let Some(first_other_way_face_key) = first_other_way_face_key else {
      return fan_normal_acc.get();
    };

    self.walk_fan_from(
      vtx,
      vtx_key,
      face_data,
      start_edge_key,
      start_edge_ix,
      *first_other_way_face_key,
      visited_edges,
      visited_faces,
      smooth_fan_faces,
      &mut fan_normal_acc,
      vtx_normal_acc,
    );

    fan_normal_acc.get()
  }

  /// Returns the normal of the full vertex, including all smooth fans.  That normal should be used
  /// for displacement.
  ///
  /// Can return `None` for cases like disconnected vertices or degenerate triangles.
  fn separate_and_compute_normals_for_vertex(
    &mut self,
    smooth_fans_acc: &mut Vec<SmoothFan>,
    face_data: &[FaceShading],
    incident_counts: &[u32],
    vtx_key: VertexKey,
  ) -> Option<Vec3> {
    let mut visited_edges = bitarr![0; 1024*8];
    let mut visited_faces = SmallVec::<[_; 16]>::new();
    let mut smooth_fan_faces: SmallVec<[FaceKey; 16]> = SmallVec::new();

    let mut vtx_normal_acc = NormalAcc::new();

    let edge_count = self.vertices[vtx_key].edges.len();
    // keeps track of which edges have been visited.  Indices match the indices of
    // `vtx.edges`
    if edge_count > 1024 * 8 {
      panic!("Vertex has too many edges; vtx_key={vtx_key:?}; edge_count={edge_count}",);
    }

    if edge_count < 3 {
      let mut seen_face_keys = SmallVec::<[_; 8]>::new();
      for &edge_key in &self.vertices[vtx_key].edges {
        for &face_key in &self.edges[edge_key].faces {
          if seen_face_keys.contains(&face_key) {
            continue;
          }
          seen_face_keys.push(face_key);

          let verts = &self.faces[face_key].vertices;
          vtx_normal_acc.accumulated_normal +=
            face_data[fkey_ix(&face_key) as usize].weighted_at(verts, vtx_key);
        }
      }
      let computed_normal = vtx_normal_acc.get();
      self.set_shading_normal(vtx_key, computed_normal);
      return computed_normal;
    }

    loop {
      let visited_edges_slice = &mut visited_edges[..edge_count];
      // Track progress to prevent infinite loops on loose edges.
      let visited_count_start = visited_edges_slice.count_ones();

      let computed_fan_normal = self.walk_one_smooth_fan(
        vtx_key,
        face_data,
        visited_edges_slice,
        &mut visited_faces,
        &mut smooth_fan_faces,
        &mut vtx_normal_acc,
      );

      let all_faces_visited = visited_edges_slice.all();
      let made_progress = visited_edges_slice.count_ones() > visited_count_start;
      if all_faces_visited || !made_progress {
        // All faces from this fan can retain the same vertex, and we can assign
        // it the computed normal directly
        self.set_shading_normal(
          vtx_key,
          computed_fan_normal.or_else(|| vtx_normal_acc.get()),
        );

        break;
      }

      if !smooth_fan_faces.is_empty() {
        let Some(computed_fan_normal) = computed_fan_normal else {
          continue;
        };

        smooth_fans_acc.push(SmoothFan {
          old_key: vtx_key,
          face_keys: smooth_fan_faces.iter().copied().collect(),
          normal: computed_fan_normal,
        });

        smooth_fan_faces.clear();
      }
    }

    // The edge walk can orphan a face whose two spokes here are both sharp and were each consumed
    // as a neighbor fan's start edge
    let incident_count =
      unsafe { *incident_counts.get_unchecked(vkey_ix(&vtx_key) as usize) } as usize;
    if core::hint::unlikely(visited_faces.len() < incident_count) {
      self.split_orphaned_fans(
        vtx_key,
        face_data,
        &mut visited_faces,
        &mut vtx_normal_acc,
        smooth_fans_acc,
      );
    }

    if cfg!(debug_assertions) && self.shading_normal(vtx_key).is_none() {
      panic!("Vertex {vtx_key:?} has no shading normal after walking smooth fans",);
    }

    vtx_normal_acc.get()
  }

  #[cold]
  fn split_orphaned_fans(
    &self,
    vtx_key: VertexKey,
    face_data: &[FaceShading],
    visited_faces: &mut SmallVec<[FaceKey; 16]>,
    vtx_normal_acc: &mut NormalAcc,
    smooth_fans_acc: &mut Vec<SmoothFan>,
  ) {
    for &seed_edge in &self.vertices[vtx_key].edges {
      for &seed in &self.edges[seed_edge].faces {
        if visited_faces.contains(&seed) {
          continue;
        }
        // Flood-fill the smooth-connected fan containing this orphan and split it off with its own
        // normal so it doesn't keep an adjacent fan's normal.
        let mut stack: SmallVec<[FaceKey; 16]> = SmallVec::new();
        stack.push(seed);
        let mut fan_faces: SmallVec<[FaceKey; 8]> = SmallVec::new();
        let mut fan_acc = NormalAcc::new();
        while let Some(fk) = stack.pop() {
          if visited_faces.contains(&fk) {
            continue;
          }
          visited_faces.push(fk);
          fan_faces.push(fk);
          let weighted =
            face_data[fkey_ix(&fk) as usize].weighted_at(&self.faces[fk].vertices, vtx_key);
          fan_acc.accumulated_normal += weighted;
          vtx_normal_acc.accumulated_normal += weighted;
          for &fedge in &self.faces[fk].edges {
            let edge = &self.edges[fedge];
            if edge.sharp || !edge.vertices.contains(&vtx_key) {
              continue;
            }
            for &adj in &edge.faces {
              if adj != fk && !visited_faces.contains(&adj) {
                stack.push(adj);
              }
            }
          }
        }
        if let Some(normal) = fan_acc.get() {
          smooth_fans_acc.push(SmoothFan {
            old_key: vtx_key,
            face_keys: fan_faces,
            normal,
          });
        }
      }
    }
  }

  fn compute_displacement_normal(
    &mut self,
    face_data: &[FaceShading],
    vtx_key: VertexKey,
  ) -> Option<Vec3> {
    let mut visited_edges = bitarr![0; 1024*8];
    // TODO: should re-use these buffers and pass them in from outer
    let mut visited_faces = SmallVec::<[_; 16]>::new();
    let mut smooth_fan_faces: SmallVec<[FaceKey; 16]> = SmallVec::new();
    let mut vtx_normal_acc = NormalAcc::new();
    let edge_count = self.vertices[vtx_key].edges.len();

    loop {
      let visited_edges = &mut visited_edges[..edge_count];

      self.walk_one_smooth_fan(
        vtx_key,
        face_data,
        visited_edges,
        &mut visited_faces,
        &mut smooth_fan_faces,
        &mut vtx_normal_acc,
      );

      let all_faces_visited = visited_edges.all();
      if all_faces_visited {
        break;
      }

      smooth_fan_faces.clear();
    }

    let displacement_normal = vtx_normal_acc.get();
    let Some(displacement_normal) = displacement_normal else {
      // panic!(
      //   "Vertex {vtx_key:?} has no displacement normal after walking smooth fans; \
      //    visited_faces.len()={}",
      //   visited_faces.len()
      // );
      // return None;
      return Some(Vec3::zeros());
    };

    Some(displacement_normal)
  }

  /// Computes displacement normals for all vertices, replacing any existing
  /// values.
  pub fn compute_vertex_displacement_normals(&mut self) {
    let (face_data, _) = self.precompute_face_shading::<false>();
    let all_vtx_keys = self.vertices.keys().collect::<Vec<_>>();
    for vtx_key in all_vtx_keys {
      let displacement_normal = self.compute_displacement_normal(&face_data, vtx_key);
      self.set_displacement_normal(vtx_key, displacement_normal);
    }
  }

  /// Does something akin to "shade auto-smooth" from Blender.  For each edge in the
  /// mesh, it is determined to be sharp or smooth based on the angle
  /// between it and the face that shares it.
  ///
  /// Then, each vertex is considered and potentially duplicated one or more
  /// times so that the normals of each are shared only with faces that are
  /// smooth wrt. each other.
  ///
  /// Heavily inspired by the Blender implementation:
  /// https://github.com/blender/blender/blob/a4aa5faa2008472413403600382f419280ac8b20/source/blender/bmesh/intern/bmesh_mesh_normals.cc#L1081
  pub fn separate_vertices_and_compute_normals(&mut self) {
    self.separate_vertices_and_compute_normals_impl::<true, false>();
  }

  /// Variant that also stores per-vertex displacement normals (a largely unused legacy feature the
  /// default no longer computes); prefer compute_vertex_displacement_normals` if you need them.
  #[doc(hidden)]
  pub fn separate_vertices_and_compute_normals_with_displacement(&mut self) {
    self.separate_vertices_and_compute_normals_impl::<true, true>();
  }

  /// In-place finalized split, leaving the mesh topologically inconsistent (`edges`/`Vertex::edges`
  /// no longer match `faces`). Bench-only; production should use separate_normals_and_finalize`.
  #[doc(hidden)]
  pub fn separate_vertices_and_compute_normals_finalized(&mut self) {
    self.separate_vertices_and_compute_normals_impl::<false, false>();
  }

  /// Separates fan vertices, computes shading normals, and exports to an indexed mesh in one shot,
  /// consuming `self`. Same output as `separate_vertices_and_compute_normals` + `to_raw_indexed`
  /// but skips the edge rebuild; the intermediate topological inconsistency can't be observed
  /// since `self` is dropped here. Preferred finalize-then-render path.
  pub fn separate_normals_and_finalize(
    mut self,
    include_shading_normals: bool,
    include_displacement_normals: bool,
    include_degenerate_faces: bool,
  ) -> OwnedIndexedMesh {
    self.separate_vertices_and_compute_normals_impl::<false, false>();
    self.to_raw_indexed(
      include_shading_normals,
      include_displacement_normals,
      include_degenerate_faces,
    )
  }

  fn separate_vertices_and_compute_normals_impl<
    const REBUILD_EDGES: bool,
    const STORE_DISPLACEMENT: bool,
  >(
    &mut self,
  ) {
    let all_vtx_keys = self.vertices.keys().collect::<Vec<_>>();
    // For each vertex in the mesh, we partition its edges into "smooth fans"
    // (as Blender calls them).  These are group of edges that all share the
    // vertex and are connected together by faces which each share one edge.
    //
    // They can either wrap all the way around the vertex or be bounded on
    // either side by a sharp edge.
    //
    // For each smooth fan, a duplicate vertex is created which gets a normal
    // computed by a weighted average of the face normals of all faces in the
    // fan.
    let (face_data, incident_counts) = self.precompute_face_shading::<true>();
    let mut smooth_fans = Vec::new();
    for vtx_key in all_vtx_keys {
      let computed_normal = self.separate_and_compute_normals_for_vertex(
        &mut smooth_fans,
        &face_data,
        &incident_counts,
        vtx_key,
      );
      if STORE_DISPLACEMENT {
        self.set_displacement_normal(vtx_key, computed_normal);
      }
    }

    // We have to wait until we've walked all the fans for the mesh before splitting
    // vertices because if we do it dynamically while walking, the mesh topology
    // will get torn up and make smooth fan walking incorrect.
    for SmoothFan {
      old_key,
      face_keys,
      normal,
    } in smooth_fans
    {
      // We have to create a new vertex for the faces that are part of the
      // smooth fan so that it can have a distinct normal
      let position = self.vertices[old_key].position;
      let new_vtx_key = self.add_vertex_cloned_from(old_key, position);
      self.set_shading_normal(new_vtx_key, Some(normal));
      if STORE_DISPLACEMENT {
        let displacement_normal = self.displacement_normal(old_key);
        self.set_displacement_normal(new_vtx_key, displacement_normal);
      }

      if REBUILD_EDGES {
        for face_key in face_keys {
          self.replace_vertex_in_face(face_key, old_key, new_vtx_key);
        }
      } else {
        // Finalized fast path: repoint the face's vertex without maintaining edges.
        for face_key in face_keys {
          let face = &mut self.faces[face_key];
          for v in &mut face.vertices {
            if *v == old_key {
              *v = new_vtx_key;
            }
          }
        }
      }
    }
  }

  /// Recomputes shading normals with dihedral auto-smooth creasing (`mark_edge_sharpness` +
  /// `separate_vertices_and_compute_normals`) while keeping authored attribute seams smooth. Any
  /// set of position-coincident verts present *before* the split (a UV seam cut, a duplicated
  /// boundary ring) is re-averaged per normal cluster afterward, so a seam running through a
  /// smooth region reads seamless while genuine geometric creases and cap edges keep their crisp
  /// split normals. Passive channels (uv/tangent) are cloned onto split verts by the separation
  /// pass and never touched here, so procedural UVs survive intact.
  ///
  /// Coincidence is bit-exact on position, matching the way seam ops (`rail_sweep`, `compute_uvs`)
  /// author their duplicates.
  pub fn recompute_shading_normals_preserving_seams(&mut self, sharp_edge_threshold_rads: f32) {
    let pos_key = |p: Vec3| [p.x.to_bits(), p.y.to_bits(), p.z.to_bits()];

    // Authored seams are the coincident verts that exist before the auto-smooth split; the split
    // itself only ever *adds* coincident verts, at fresh creases.
    let mut counts: FxHashMap<[u32; 3], u32> = FxHashMap::default();
    for v in self.vertices.values() {
      *counts.entry(pos_key(v.position)).or_insert(0) += 1;
    }
    let seam_positions: FxHashSet<[u32; 3]> = counts
      .into_iter()
      .filter(|&(_, n)| n > 1)
      .map(|(k, _)| k)
      .collect();

    self.mark_edge_sharpness(sharp_edge_threshold_rads);
    self.separate_vertices_and_compute_normals();

    if seam_positions.is_empty() {
      return;
    }

    // Re-gather the (possibly further crease-split) verts sitting on each seam position.
    let mut groups: FxHashMap<[u32; 3], Vec<VertexKey>> = FxHashMap::default();
    for (key, v) in self.vertices.iter() {
      let p = pos_key(v.position);
      if seam_positions.contains(&p) {
        groups.entry(p).or_default().push(key);
      }
    }

    // Within each seam, greedily cluster verts whose recomputed normals stay within the smooth
    // threshold and collapse each multi-member cluster onto its averaged normal; distinct clusters
    // (real creases, cap turns) stay sharp.
    let cos_threshold = sharp_edge_threshold_rads.cos();
    let mut clusters: Vec<(Vec3, Vec<VertexKey>)> = Vec::new();
    for verts in groups.values() {
      if verts.len() < 2 {
        continue;
      }
      clusters.clear();
      for &vk in verts {
        let Some(n) = self.shading_normal(vk) else {
          continue;
        };
        match clusters
          .iter_mut()
          .find(|(rep, _)| rep.dot(&n) >= cos_threshold)
        {
          Some((_, members)) => members.push(vk),
          None => clusters.push((n, vec![vk])),
        }
      }
      for (_, members) in &clusters {
        if members.len() < 2 {
          continue;
        }
        let sum: Vec3 = members
          .iter()
          .filter_map(|&vk| self.shading_normal(vk))
          .sum();
        if let Some(avg) = sum.try_normalize(1e-12) {
          for &vk in members {
            self.set_shading_normal(vk, Some(avg));
          }
        }
      }
    }
  }

  /// Passive channels sorted by name with their export arity, plus a scratch row for `add_vtx`.
  fn export_channels(&self) -> (Vec<(&String, &Channel<VertexKey>)>, Vec<(String, usize)>) {
    let mut channels: Vec<_> = self.vertex_channels.iter().collect();
    channels.sort_unstable_by(|a, b| a.0.cmp(b.0));
    let specs = channels
      .iter()
      .map(|(name, ch)| {
        let arity = match ch.store {
          ChannelStore::Scalar(_) => Arity::Scalar,
          ChannelStore::Vec2(_) => Arity::Vec2,
          ChannelStore::Vec3(_) => Arity::Vec3,
          ChannelStore::Vec4(_) => Arity::Vec4,
        };
        ((*name).clone(), attrs::export_arity(name, arity))
      })
      .collect();
    (channels, specs)
  }

  fn add_vtx_with_channels(
    &self,
    builder: &mut OwnedIndexedMeshBuilder,
    channels: &[(&String, &Channel<VertexKey>)],
    vals: &mut [[f32; 4]],
    vtx_key: VertexKey,
  ) {
    for (slot, (name, ch)) in vals.iter_mut().zip(channels) {
      *slot = attrs::export_value(name, ch.get(vtx_key));
    }
    builder.add_vtx(
      vtx_key,
      self.vertices[vtx_key].position,
      self.shading_normal(vtx_key),
      self.displacement_normal(vtx_key),
      vals,
    );
  }

  pub fn to_raw_indexed(
    &self,
    include_shading_normals: bool,
    include_displacement_normals: bool,
    include_degenerate_faces: bool,
  ) -> OwnedIndexedMesh {
    let (channels, specs) = self.export_channels();
    let mut builder = OwnedIndexedMeshBuilder::with_capacity(
      self.vertices.len(),
      self.faces.len(),
      include_displacement_normals,
      include_shading_normals,
      &specs,
    );
    let mut vals = vec![[0f32; 4]; channels.len()];

    for face in self.faces.values() {
      if !include_degenerate_faces && face.is_degenerate(&self.vertices) {
        continue;
      }
      for &vtx_key in &face.vertices {
        self.add_vtx_with_channels(&mut builder, &channels, &mut vals, vtx_key);
      }
    }

    builder.build(self.transform)
  }

  /// Works the same as `to_raw_indexed` but splits the triangles into multiple different meshes
  /// based on `partition_fn`.
  pub fn to_raw_indexed_multi<T: Hash + Eq>(
    &self,
    include_shading_normals: bool,
    include_displacement_normals: bool,
    partition_fn: impl Fn(&Face<FaceData>) -> T,
  ) -> Vec<OwnedIndexedMesh> {
    let (channels, specs) = self.export_channels();
    let mut vals = vec![[0f32; 4]; channels.len()];
    let mut out_meshes = FxHashMap::default();

    for face in self.faces.values() {
      if face.is_degenerate(&self.vertices) {
        continue;
      }

      let out_key = partition_fn(face);
      let builder = out_meshes.entry(out_key).or_insert_with(|| {
        OwnedIndexedMeshBuilder::new(
          include_displacement_normals,
          include_shading_normals,
          &specs,
        )
      });

      for &vtx_key in &face.vertices {
        self.add_vtx_with_channels(builder, &channels, &mut vals, vtx_key);
      }
    }

    out_meshes
      .into_values()
      .filter(|builder| !builder.is_empty())
      .map(|builder| builder.build(self.transform))
      .collect()
  }

  pub fn to_owned_mesh(&self, transform: Option<Mat4>) -> OwnedMesh {
    OwnedMesh {
      vertices: self
        .faces
        .values()
        .flat_map(|face| face.vertices.map(|vtx_key| self.vertices[vtx_key].position))
        .collect(),
      normals: Some(
        self
          .faces
          .values()
          .flat_map(|face| {
            face.vertices.map(|vtx_key| {
              self
                .displacement_normal(vtx_key)
                .unwrap_or_else(Vec3::zeros)
            })
          })
          .collect(),
      ),
      transform,
    }
  }

  /// Splits an edge in half, creating a new vertex in the middle.  All faces
  /// that include this edge are split into two faces with new edges being
  /// created for each.
  ///
  /// Returns the key of the new vertex.
  pub fn split_edge_cb(
    &mut self,
    edge_key_to_split: EdgeKey,
    split_pos: EdgeSplitPos,
    displacement_normal_method: DisplacementNormalMethod,
    mut face_split_cb: impl FnMut(&Self, FaceKey, FaceData, [FaceKey; 2]),
  ) -> VertexKey {
    let edge = self.edges.get(edge_key_to_split).unwrap_or_else(|| {
      panic!("Tried to split edge that doesn't exist; key={edge_key_to_split:?}")
    });
    let edge_id_to_split = edge.vertices;
    let [v0_key, v1_key] = edge_id_to_split;
    let [v0, v1] = [&self.vertices[v0_key], &self.vertices[v1_key]];

    let split_pos = split_pos.get(v0_key);
    let vm_position = v0.position.lerp(&v1.position, split_pos);

    let edge_displacement_normal = self
      .edge_displacement_normals
      .get(edge_key_to_split)
      .copied();
    let faces_to_split = edge.faces.clone();

    let shading_normal = self
      .shading_normals
      .get(edge_id_to_split[0])
      .copied()
      .zip(self.shading_normals.get(edge_id_to_split[1]).copied())
      .map(|(n0, n1)| n0.lerp(&n1, split_pos).normalize());

    let displacement_normal = match displacement_normal_method {
      DisplacementNormalMethod::Interpolate => {
        match self
          .displacement_normals
          .get(edge_id_to_split[0])
          .copied()
          .zip(self.displacement_normals.get(edge_id_to_split[1]).copied())
        {
          Some((n0, n1)) => {
            let merged_normal = n0.lerp(&n1, split_pos).normalize();
            if merged_normal.x.is_nan() || merged_normal.y.is_nan() || merged_normal.z.is_nan() {
              panic!("Merged normal is NaN; n0={n0:?}; n1={n1:?}; merged_normal={merged_normal:?}");
            } else {
              Some(merged_normal)
            }
          }
          None => None,
        }
      }
      // The new vertex inherits the edge's displacement normal directly. Averaging the two endpoint
      // vertices instead pulls in unrelated faces' normals, inflating straight edges into curves.
      DisplacementNormalMethod::EdgeNormal => edge_displacement_normal,
    };

    let middle_vertex_key = self.add_vertex_interpolated(
      &[
        (edge_id_to_split[0], 1. - split_pos),
        (edge_id_to_split[1], split_pos),
      ],
      vm_position,
    );
    self.set_shading_normal(middle_vertex_key, shading_normal);
    self.set_displacement_normal(middle_vertex_key, displacement_normal);

    // Split each adjacent face
    for old_face_key in faces_to_split {
      let old_face = &self.faces[old_face_key];
      let old_face_normal = old_face.normal(&self.vertices);
      let [v0_key, v1_key, v2_key] = old_face.vertices;
      let split_edge_ix = if edge_id_to_split == sort_edge(v0_key, v1_key) {
        0
      } else if edge_id_to_split == sort_edge(v1_key, v2_key) {
        1
      } else {
        2
      };

      let vertex_keys = [v0_key, v1_key, v2_key, middle_vertex_key];
      let edge_displacement_normals = [
        self
          .edge_displacement_normals
          .get(old_face.edges[0])
          .copied(),
        self
          .edge_displacement_normals
          .get(old_face.edges[1])
          .copied(),
        self
          .edge_displacement_normals
          .get(old_face.edges[2])
          .copied(),
        Some(old_face_normal),
      ];
      let edge_sharpnesses = [
        self.edges[old_face.edges[0]].sharp,
        self.edges[old_face.edges[1]].sharp,
        self.edges[old_face.edges[2]].sharp,
        false,
      ];

      let orders = match split_edge_ix {
        0 => [[0, 3, 2], [3, 1, 2]],
        1 => [[0, 1, 3], [0, 3, 2]],
        2 => [[0, 1, 3], [1, 2, 3]],
        _ => unreachable!(),
      };

      let old_face_data = self.remove_face(old_face_key);

      let mut add_face = |order_ix: usize| {
        let order = &orders[order_ix];
        let face_key = self.add_face::<false>(
          [
            vertex_keys[order[0]],
            vertex_keys[order[1]],
            vertex_keys[order[2]],
          ],
          Default::default(),
        );
        let edge_keys = &self.faces[face_key].edges;
        for edge_ix in 0..3 {
          let edge_key = edge_keys[edge_ix];
          self.edges[edge_key].sharp = edge_sharpnesses[order[edge_ix]];
          set_channel(
            &mut self.edge_displacement_normals,
            edge_key,
            edge_displacement_normals[order[edge_ix]],
          );
        }
        face_key
      };

      let new_face_keys = [add_face(0), add_face(1)];
      face_split_cb(&*self, old_face_key, old_face_data, new_face_keys);
    }

    assert!(self.remove_edge(edge_key_to_split).is_none());

    middle_vertex_key
  }

  /// Splits an edge in half, creating a new vertex in the middle.  All faces
  /// that include this edge are split into two faces with new edges being
  /// created for each.
  ///
  /// Returns the key of the new vertex.
  pub fn split_edge(
    &mut self,
    edge_key_to_split: EdgeKey,
    split_pos: EdgeSplitPos,
    displacement_normal_method: DisplacementNormalMethod,
  ) -> VertexKey {
    self.split_edge_cb(
      edge_key_to_split,
      split_pos,
      displacement_normal_method,
      |_, _, _, _| {},
    )
  }

  pub fn get_edge_key(&self, mut vertices: [VertexKey; 2]) -> Option<EdgeKey> {
    sort_edge_inplace::<false>(&mut vertices);
    self.get_edge_key_from_sorted::<false>(vertices)
  }

  pub fn get_edge_key_from_sorted<const DENSE: bool>(
    &self,
    sorted_vtx_keys: [VertexKey; 2],
  ) -> Option<EdgeKey> {
    // iter through the vtx with fewer edges
    let vtx0 = unsafe { self.vertices.get_unchecked(sorted_vtx_keys[0]) };
    let vtx1 = unsafe { self.vertices.get_unchecked(sorted_vtx_keys[1]) };
    let vtx = if vtx0.edges.len() < vtx1.edges.len() {
      vtx0
    } else {
      vtx1
    };

    vtx
      .edges
      .iter()
      .find(|&&edge_key| {
        let edge = unsafe { self.edges.get_unchecked(edge_key) };
        if DENSE {
          vkey_ix(&edge.vertices[0]) == vkey_ix(&sorted_vtx_keys[0])
            && vkey_ix(&edge.vertices[1]) == vkey_ix(&sorted_vtx_keys[1])
        } else {
          edge.vertices == sorted_vtx_keys
        }
      })
      .copied()
  }

  /// Returns the key of the edge between the two provided vertices if it exists and creates it if
  /// not.
  ///
  /// If a new edge was created, the vertices of the edge will be updated to reference it.
  ///
  /// If `DENSE` is true, assumes that the `edges` slotmap is dense, meaning that it has never had
  /// elements removed from it, and uses a fastpath method for inserting new edges.
  pub fn get_or_create_edge<const DENSE: bool>(&mut self, mut vertices: [VertexKey; 2]) -> EdgeKey {
    sort_edge_inplace::<DENSE>(&mut vertices);

    match self.get_edge_key_from_sorted::<DENSE>(vertices) {
      Some(edge_key) => edge_key,
      None => {
        let new_edge = Edge {
          vertices,
          faces: SmallVec::new(),
          sharp: false,
        };
        let edge_key = if DENSE {
          slotmap_insert_dense(&mut self.edges, new_edge)
        } else {
          self.edges.insert(new_edge)
        };

        for vert_key in vertices {
          let vert = unsafe { self.vertices.get_unchecked_mut(vert_key) };
          vert.edges.push(edge_key);
        }

        edge_key
      }
    }
  }

  pub fn add_face<const DENSE: bool>(
    &mut self,
    vertices: [VertexKey; 3],
    data: FaceData,
  ) -> FaceKey {
    let edges = [(0, 1), (1, 2), (2, 0)];
    let edge_keys =
      edges.map(|(a, b)| self.get_or_create_edge::<DENSE>([vertices[a], vertices[b]]));

    let cb = |face_key| {
      for edge_key in edge_keys {
        let edge = unsafe { self.edges.get_unchecked_mut(edge_key) };
        edge.faces.push(face_key);
      }

      Face {
        vertices,
        edges: edge_keys,
        data,
      }
    };

    if DENSE {
      slotmap_insert_dense_with_key(&mut self.faces, cb)
    } else {
      self.faces.insert_with_key(cb)
    }
  }

  pub fn remove_face(&mut self, face_key: FaceKey) -> FaceData {
    let face = self
      .faces
      .remove(face_key)
      .unwrap_or_else(|| panic!("Tried to remove face that doesn't exist; key={face_key:?}"));

    for edge_key in face.edges {
      let edge = self.edges.get_mut(edge_key).unwrap_or_else(|| {
        panic!(
          "Tried to get edge that doesn't exist; key={edge_key:?}. Was referenced by removed face \
           with key={face_key:?}",
        )
      });
      swap_retain_sv(&mut edge.faces, |&mut f| f != face_key);
      if edge.faces.is_empty() {
        let edge = self.remove_edge(edge_key).unwrap();
        for vert_key in edge.vertices {
          let vert = &mut self.vertices[vert_key];
          swap_retain_sv(&mut vert.edges, |&mut e| e != edge_key);
        }
      }
    }

    face.data
  }

  pub fn connected_components(&self) -> Vec<Vec<FaceKey>> {
    let mut visited =
      FxHashSet::with_capacity_and_hasher(self.faces.len(), FxBuildHasher::default());
    let mut components = Vec::new();

    for face_key in self.faces.keys() {
      if visited.contains(&face_key) {
        continue;
      }

      let mut component = Vec::new();
      let mut stack = vec![face_key];

      while let Some(cur_face_key) = stack.pop() {
        if visited.insert(cur_face_key) {
          component.push(cur_face_key);
          let face = &self.faces[cur_face_key];
          for &edge_key in &face.edges {
            for &adj_face_key in &self.edges[edge_key].faces {
              if !visited.contains(&adj_face_key) {
                stack.push(adj_face_key);
              }
            }
          }
        }
      }

      components.push(component);
    }

    components
  }

  /// For each connected component in the mesh, checks if the component is closed.  If it is, it
  /// determines if the component is oriented outward (if the face normals are facing the right
  /// way).  If they are not, the winding order and vertex normals for the component are flipped in
  /// place.
  pub fn orient_connected_components_outward(&mut self) {
    let components = self.connected_components();
    for component in components {
      let is_closed = component
        .iter()
        .flat_map(|&face_key| {
          let face = &self.faces[face_key];
          face.edges.iter().map(|&edge_key| &self.edges[edge_key])
        })
        .all(|edge| edge.faces.len() == 2);

      if !is_closed {
        continue;
      }

      let uniq_vtx_keys: FxHashSet<_> = component
        .iter()
        .flat_map(|&face_key| {
          let face = &self.faces[face_key];
          face.vertices.iter().copied()
        })
        .collect();

      let signed_volume_of_component = {
        let mut v = 0.0f32;
        for &face_key in &component {
          let face = &self.faces[face_key];
          let [i0, i1, i2] = face.vertices;
          let a = self.vertices[i0].position;
          let b = self.vertices[i1].position;
          let c = self.vertices[i2].position;
          v += a.cross(&b).dot(&c);
        }
        v / 6.0
      };
      let is_outward_facing = signed_volume_of_component > 0.;

      if is_outward_facing {
        continue;
      }

      // Flip the winding order and vertex normals of all faces in this component
      for &face_key in &component {
        let face = &mut self.faces[face_key];
        face.vertices.swap(0, 2);
      }

      for vtx_key in uniq_vtx_keys {
        if let Some(n) = self.shading_normals.get_mut(vtx_key) {
          *n = -*n;
        }
        if let Some(n) = self.displacement_normals.get_mut(vtx_key) {
          *n = -*n;
        }
        for ch in self.vertex_channels.values_mut() {
          ch.store.flip_slot(vtx_key, ch.flip);
        }
      }
    }
  }

  /// Removes all faces from the mesh that contain an edge is not shared by exactly two faces.
  pub fn remove_nonmanifold_faces(&mut self) -> usize {
    let mut total_removed_count = 0;

    loop {
      let nonmanifold_face_keys: Vec<FaceKey> = self
        .faces
        .iter()
        .filter_map(|(face_key, face)| {
          let has_nonmanifold_edge = face.edges.iter().any(|&edge_key| {
            let edge = &self.edges[edge_key];
            edge.faces.len() != 2
          });
          if has_nonmanifold_edge {
            Some(face_key)
          } else {
            None
          }
        })
        .collect();

      let removed_count = nonmanifold_face_keys.len();
      if removed_count == 0 {
        break;
      }
      total_removed_count += removed_count;

      for face_key in nonmanifold_face_keys {
        self.remove_face(face_key);
      }
    }

    total_removed_count
  }

  /// Removes all degenerate faces from the mesh.
  pub fn cleanup_degenerate_triangles_cb(&mut self, mut cb: impl FnMut(FaceKey, FaceData)) {
    let all_face_keys = self.faces.keys().collect::<Vec<_>>();
    for face_key in all_face_keys {
      let face = &self.faces[face_key];
      if face.is_degenerate(&self.vertices) {
        let face_data = self.remove_face(face_key);
        cb(face_key, face_data);
      }
    }
  }

  pub fn cleanup_degenerate_triangles(&mut self) {
    self.cleanup_degenerate_triangles_cb(|_, _| {});
  }

  pub fn is_empty(&self) -> bool {
    self.vertices.is_empty()
  }

  pub fn subdivide_by_plane(&mut self, plane: &Plane) {
    let face_keys: Vec<FaceKey> = self.faces.keys().collect();
    let mut created_vertices = FxHashMap::<[VertexKey; 2], VertexKey>::default();

    for face_key in face_keys {
      if !self.faces.contains_key(face_key) {
        continue;
      }

      let face = &self.faces[face_key];
      let mut polygon_type = PolygonClass::Coplanar;
      let mut types = [PolygonClass::Coplanar; 3];
      let mut dists = [0.; 3];

      for (i, &v_key) in face.vertices.iter().enumerate() {
        let dist = plane.normal.dot(&self.vertices[v_key].position) - plane.w;
        dists[i] = dist;
        types[i] = if dist < -EPSILON {
          PolygonClass::Back
        } else if dist > EPSILON {
          PolygonClass::Front
        } else {
          PolygonClass::Coplanar
        };
        polygon_type = polygon_type | types[i];
      }

      if polygon_type == PolygonClass::Spanning {
        let mut f: ArrayVec<VertexKey, 4> = ArrayVec::new();
        let mut b: ArrayVec<VertexKey, 4> = ArrayVec::new();
        let face_verts = self.faces[face_key].vertices;

        for i in 0..3 {
          let j = (i + 1) % 3;
          let ti = types[i];
          let tj = types[j];
          let vi_key = face_verts[i];
          let vj_key = face_verts[j];

          if ti != PolygonClass::Back {
            f.push(vi_key);
          }
          if ti != PolygonClass::Front {
            b.push(vi_key);
          }

          if (ti | tj) == PolygonClass::Spanning {
            let edge_key = sort_edge(vi_key, vj_key);
            let new_v_key = if let Some(&key) = created_vertices.get(&edge_key) {
              key
            } else {
              let t = dists[i] / (dists[i] - dists[j]);
              let new_v_pos = self.vertices[vi_key]
                .position
                .lerp(&self.vertices[vj_key].position, t);

              let lerp_normalize =
                |ch: &SecondaryMap<VertexKey, Vec3>| match (ch.get(vi_key), ch.get(vj_key)) {
                  (Some(n0), Some(n1)) => Some(n0.lerp(n1, t).normalize()),
                  _ => None,
                };
              let shading_normal = lerp_normalize(&self.shading_normals);
              let displacement_normal = lerp_normalize(&self.displacement_normals);

              let key = self.add_vertex_interpolated(&[(vi_key, 1. - t), (vj_key, t)], new_v_pos);
              self.set_shading_normal(key, shading_normal);
              self.set_displacement_normal(key, displacement_normal);
              created_vertices.insert(edge_key, key);
              key
            };

            f.push(new_v_key);
            b.push(new_v_key);
          }
        }

        self.remove_face(face_key);

        if f.len() >= 3 {
          self.add_face::<false>([f[0], f[1], f[2]], Default::default());
          if f.len() == 4 {
            self.add_face::<false>([f[0], f[2], f[3]], Default::default());
          }
        }

        if b.len() >= 3 {
          self.add_face::<false>([b[0], b[1], b[2]], Default::default());
          if b.len() == 4 {
            self.add_face::<false>([b[0], b[2], b[3]], Default::default());
          }
        }
      }
    }
  }

  /// Skewers the mesh with an infinite line.  `eps` is a world-space snap tolerance for the
  /// hit point against nearby vertices/edges.  No-op for lines coplanar with a triangle's plane.
  pub fn subdivide_by_line(&mut self, origin: Vec3, dir: Vec3, eps: f32) {
    if dir.magnitude_squared() < 1e-12 || !eps.is_finite() || eps <= 0. {
      return;
    }
    let dir = dir.normalize();

    let face_keys: Vec<FaceKey> = self.faces.keys().collect();
    let mut edge_events: FxHashMap<EdgeKey, (f32, u32)> = FxHashMap::default();
    let mut interior_events: Vec<(FaceKey, Vec3, [f32; 3])> = Vec::new();

    for face_key in face_keys {
      let face = &self.faces[face_key];
      let [va_key, vb_key, vc_key] = face.vertices;
      let a = self.vertices[va_key].position;
      let b = self.vertices[vb_key].position;
      let c = self.vertices[vc_key].position;

      use crate::triangle_intersection::{classify_line_triangle, LineTriHit};
      match classify_line_triangle(origin, dir, a, b, c, eps) {
        LineTriHit::None | LineTriHit::Vertex => {}
        LineTriHit::Edge {
          face_edge_ix,
          t_along,
        } => {
          let edge_key = face.edges[face_edge_ix];
          let face_v0 = face.vertices[face_edge_ix];
          let edge_v0 = self.edges[edge_key].vertices[0];
          let t_canonical = if face_v0 == edge_v0 {
            t_along
          } else {
            1. - t_along
          };
          let entry = edge_events.entry(edge_key).or_insert((0., 0));
          entry.0 += t_canonical;
          entry.1 += 1;
        }
        LineTriHit::Interior { point, bary } => {
          interior_events.push((face_key, point, bary));
        }
      }
    }

    for (edge_key, (sum_t, count)) in edge_events {
      let t_avg = (sum_t / count as f32).clamp(0., 1.);
      let start_vtx_key = self.edges[edge_key].vertices[0];
      self.split_edge(
        edge_key,
        EdgeSplitPos {
          pos: t_avg,
          start_vtx_key,
        },
        DisplacementNormalMethod::Interpolate,
      );
    }

    for (face_key, point, bary) in interior_events {
      // Skip if a prior edge_event already removed this face (adjacent-face classifier
      // disagreement near a shared edge under f32 jitter).
      if self.faces.contains_key(face_key) {
        self.split_face_interior(face_key, point, bary);
      }
    }
  }

  fn split_face_interior(&mut self, face_key: FaceKey, point: Vec3, bary: [f32; 3]) {
    let face = &self.faces[face_key];
    let face_normal = face.normal(&self.vertices);
    let [va, vb, vc] = face.vertices;
    let [alpha, beta, gamma] = bary;
    let outer_edge_data: [(Option<Vec3>, bool); 3] = std::array::from_fn(|i| {
      let edge_key = face.edges[i];
      (
        self.edge_displacement_normals.get(edge_key).copied(),
        self.edges[edge_key].sharp,
      )
    });

    let lerp3 = |na: Vec3, nb: Vec3, nc: Vec3| {
      (na * alpha + nb * beta + nc * gamma)
        .try_normalize(0.)
        .unwrap_or(face_normal)
    };
    let blend = |ch: &SecondaryMap<VertexKey, Vec3>| {
      ch.get(va)
        .copied()
        .zip(ch.get(vb).copied())
        .zip(ch.get(vc).copied())
        .map(|((na, nb), nc)| lerp3(na, nb, nc))
    };
    let shading_normal = blend(&self.shading_normals);
    let displacement_normal = blend(&self.displacement_normals);

    let m_key = self.add_vertex_interpolated(&[(va, alpha), (vb, beta), (vc, gamma)], point);
    self.set_shading_normal(m_key, shading_normal);
    self.set_displacement_normal(m_key, displacement_normal);

    self.remove_face(face_key);

    for (i, (v0, v1)) in [(va, vb), (vb, vc), (vc, va)].into_iter().enumerate() {
      let new_face_key = self.add_face::<false>([v0, v1, m_key], Default::default());
      let new_face_edges = self.faces[new_face_key].edges;
      let outer_key = new_face_edges[0];
      let (outer_n, outer_sharp) = outer_edge_data[i];
      if !self.edge_displacement_normals.contains_key(outer_key) {
        set_channel(&mut self.edge_displacement_normals, outer_key, outer_n);
      }
      self.edges[outer_key].sharp |= outer_sharp;
      for &edge_key in &[new_face_edges[1], new_face_edges[2]] {
        if !self.edge_displacement_normals.contains_key(edge_key) {
          self.edge_displacement_normals.insert(edge_key, face_normal);
        }
      }
    }
  }

  pub fn compute_aabb(&self, transform: &Matrix4<f32>) -> Aabb {
    Aabb::from_points(
      self
        .vertices
        .values()
        .map(|v| (transform * v.position.push(1.)).xyz()),
    )
  }

  /// World-space BVH over the faces in `faces.values()` order (degenerate ones included).
  pub fn build_bvh(&self, transform: &Matrix4<f32>) -> TriBvh {
    let xf = |v: Vec3| (transform * v.push(1.)).xyz();
    TriBvh::build(
      self
        .faces
        .values()
        .map(|f| f.vertices.map(|k| xf(self.vertices[k].position)))
        .collect(),
    )
  }

  pub fn new_box(width: f32, height: f32, depth: f32) -> Self {
    let half_width = width / 2.;
    let half_height = height / 2.;
    let half_depth = depth / 2.;

    LinkedMesh::from_indexed_vertices(
      &[
        Vec3::new(-half_width, -half_height, -half_depth),
        Vec3::new(half_width, -half_height, -half_depth),
        Vec3::new(half_width, half_height, -half_depth),
        Vec3::new(-half_width, half_height, -half_depth),
        Vec3::new(-half_width, -half_height, half_depth),
        Vec3::new(half_width, -half_height, half_depth),
        Vec3::new(half_width, half_height, half_depth),
        Vec3::new(-half_width, half_height, half_depth),
      ],
      &[
        0, 2, 1, 0, 3, 2, //
        4, 5, 6, 4, 6, 7, //
        0, 1, 5, 0, 5, 4, //
        2, 3, 7, 2, 7, 6, //
        0, 7, 3, 0, 4, 7, //
        1, 2, 6, 1, 6, 5, //
      ],
      None,
      None,
    )
  }

  /// Generates an icosphere mesh with the given `radius` and number of recursive `subdivisions`.
  pub fn new_icosphere(radius: f32, subdivisions: u32) -> Self {
    // let t = (1. + 5.0_f32.sqrt()) / 2.;
    // let mut vertices = vec![
    //   Vec3::new(-1., t, 0.).normalize(),
    //   Vec3::new(1., t, 0.).normalize(),
    //   Vec3::new(-1., -t, 0.).normalize(),
    //   Vec3::new(1., -t, 0.).normalize(),
    //   Vec3::new(0., -1., t).normalize(),
    //   Vec3::new(0., 1., t).normalize(),
    //   Vec3::new(0., -1., -t).normalize(),
    //   Vec3::new(0., 1., -t).normalize(),
    //   Vec3::new(t, 0., -1.).normalize(),
    //   Vec3::new(t, 0., 1.).normalize(),
    //   Vec3::new(-t, 0., -1.).normalize(),
    //   Vec3::new(-t, 0., 1.).normalize(),
    // ];

    let total_vtx_count = 10 * (4usize).pow(subdivisions) + 2;

    let mut faces: Vec<[usize; 3]> = Vec::new();
    let mut vertices: Vec<Vec3> = Vec::with_capacity(total_vtx_count);
    vertices.extend_from_slice(ICOSAHEDRON_VERTICES);

    if subdivisions == 0 {
      faces = ICOSAHEDRON_FACES.to_vec();
    }

    for _ in 0..subdivisions {
      let old_faces: &[_] = if faces.is_empty() {
        ICOSAHEDRON_FACES
      } else {
        &faces
      };
      let new_face_count = old_faces.len() * 4;
      let mut new_faces = Vec::with_capacity(new_face_count);
      let mut midpoint_cache = FxHashMap::<(usize, usize), usize>::default();

      let mut get_midpoint = |i1: usize, i2: usize| {
        let key = if i1 < i2 { (i1, i2) } else { (i2, i1) };
        if let Some(&idx) = midpoint_cache.get(&key) {
          return idx;
        }
        let v1 = *unsafe { vertices.get_unchecked(i1) };
        let v2 = *unsafe { vertices.get_unchecked(i2) };
        let midpoint = ((v1 + v2) * 0.5).normalize();
        let vtx_ix = vertices.len();
        vertices.push(midpoint);
        midpoint_cache.insert(key, vtx_ix);
        vtx_ix
      };

      for &[i0, i1, i2] in old_faces {
        let a = get_midpoint(i0, i1);
        let b = get_midpoint(i1, i2);
        let c = get_midpoint(i2, i0);

        new_faces.push([i0, a, c]);
        new_faces.push([i1, b, a]);
        new_faces.push([i2, c, b]);
        new_faces.push([a, b, c]);
      }

      assert_eq!(new_faces.len(), new_face_count);
      faces = new_faces;
    }

    assert_eq!(vertices.len(), total_vtx_count);
    for v in &mut vertices {
      *v *= radius;
    }

    if std::mem::size_of::<usize>() == std::mem::size_of::<u32>() {
      let indices: &[u32] =
        unsafe { std::slice::from_raw_parts(faces.as_ptr() as *const u32, faces.len() * 3) };
      return LinkedMesh::from_indexed_vertices(&vertices, indices, None, None);
    }

    let indices: Vec<u32> = faces
      .into_iter()
      .flat_map(|[i0, i1, i2]| [i0 as u32, i1 as u32, i2 as u32])
      .collect();

    LinkedMesh::from_indexed_vertices(&vertices, &indices, None, None)
  }

  /// Regular octahedron with circumradius `radius`; vertices sit on the ±X/±Y/±Z axes, so a vertex
  /// points straight up.
  pub fn new_octahedron(radius: f32) -> Self {
    let verts = [
      Vec3::new(radius, 0., 0.),
      Vec3::new(-radius, 0., 0.),
      Vec3::new(0., radius, 0.),
      Vec3::new(0., -radius, 0.),
      Vec3::new(0., 0., radius),
      Vec3::new(0., 0., -radius),
    ];
    let mut indices = Vec::with_capacity(24);
    for &x in &[0usize, 1] {
      for &y in &[2usize, 3] {
        for &z in &[4usize, 5] {
          indices.extend_from_slice(&outward_tri(&verts, x, y, z));
        }
      }
    }
    LinkedMesh::from_indexed_vertices(&verts, &indices, None, None)
  }

  /// Regular tetrahedron with circumradius `radius`, inscribed in the cube (vertices at four
  /// alternating cube corners; 2-fold symmetric about each axis).
  pub fn new_tetrahedron(radius: f32) -> Self {
    let k = radius / 3f32.sqrt();
    let verts = [
      Vec3::new(k, k, k),
      Vec3::new(k, -k, -k),
      Vec3::new(-k, k, -k),
      Vec3::new(-k, -k, k),
    ];
    let mut indices = Vec::with_capacity(12);
    for [a, b, c] in [[0, 1, 2], [0, 1, 3], [0, 2, 3], [1, 2, 3]] {
      indices.extend_from_slice(&outward_tri(&verts, a, b, c));
    }
    LinkedMesh::from_indexed_vertices(&verts, &indices, None, None)
  }

  /// Regular icosahedron with circumradius `radius` (the icosphere's base mesh).
  pub fn new_icosahedron(radius: f32) -> Self {
    Self::new_icosphere(radius, 0)
  }

  /// Regular dodecahedron with circumradius `radius`, built as the dual of the icosahedron: dual
  /// vertices are the icosahedron's face centroids; dual faces are the pentagons ringing each
  /// icosahedron vertex.
  pub fn new_dodecahedron(radius: f32) -> Self {
    let verts: Vec<Vec3> = ICOSAHEDRON_FACES
      .iter()
      .map(|&[a, b, c]| {
        (ICOSAHEDRON_VERTICES[a] + ICOSAHEDRON_VERTICES[b] + ICOSAHEDRON_VERTICES[c]).normalize()
          * radius
      })
      .collect();
    let mut indices: Vec<u32> = Vec::with_capacity(108);
    for vi in 0..ICOSAHEDRON_VERTICES.len() {
      let mut ring: Vec<usize> = (0..ICOSAHEDRON_FACES.len())
        .filter(|&f| ICOSAHEDRON_FACES[f].contains(&vi))
        .collect();
      // order the ring cyclically by angle in the plane perpendicular to the vertex direction
      let axis = ICOSAHEDRON_VERTICES[vi];
      let seed = if axis.x.abs() < 0.9 {
        Vec3::x()
      } else {
        Vec3::y()
      };
      let u = (seed - axis * seed.dot(&axis)).normalize();
      let w = axis.cross(&u);
      ring.sort_by(|&f1, &f2| {
        let a1 = verts[f1].dot(&w).atan2(verts[f1].dot(&u));
        let a2 = verts[f2].dot(&w).atan2(verts[f2].dot(&u));
        a1.partial_cmp(&a2).unwrap()
      });
      for k in 1..ring.len() - 1 {
        indices.extend_from_slice(&outward_tri(&verts, ring[0], ring[k], ring[k + 1]));
      }
    }
    LinkedMesh::from_indexed_vertices(&verts, &indices, None, None)
  }

  /// N-gonal bipyramid: a regular `n`-gon equator of circumradius `radius` in the XZ plane with
  /// apexes at ±`height` on the Y axis.  `new_bipyramid(4, r, r)` is the axis-aligned octahedron.
  pub fn new_bipyramid(n: usize, radius: f32, height: f32) -> Self {
    let mut verts: Vec<Vec3> = (0..n)
      .map(|i| {
        let a = std::f32::consts::TAU * (i as f32) / (n as f32);
        Vec3::new(radius * a.cos(), 0., radius * a.sin())
      })
      .collect();
    let top = verts.len();
    verts.push(Vec3::new(0., height, 0.));
    let bot = verts.len();
    verts.push(Vec3::new(0., -height, 0.));
    let mut indices = Vec::with_capacity(n * 6);
    for i in 0..n {
      let j = (i + 1) % n;
      indices.extend_from_slice(&outward_tri(&verts, top, i, j));
      indices.extend_from_slice(&outward_tri(&verts, bot, i, j));
    }
    LinkedMesh::from_indexed_vertices(&verts, &indices, None, None)
  }

  pub fn new_cylinder(
    radius: f32,
    height: f32,
    radial_segments: usize,
    height_segments: usize,
  ) -> Self {
    let mut vertices = Vec::new();
    let mut indices = Vec::new();

    generate_radial_vertices(
      &mut vertices,
      radial_segments,
      height_segments,
      |_h| radius,
      |h| -height / 2. + h * height,
    );

    generate_radial_side_indices(&mut indices, radial_segments, height_segments);

    let bottom_center_ix = vertices.len() as u32;
    vertices.push(Vec3::new(0., -height / 2., 0.));
    let top_center_ix = vertices.len() as u32;
    vertices.push(Vec3::new(0., height / 2., 0.));

    for (center_ix, reverse_winding, ring_offset) in [
      (bottom_center_ix, false, 0usize),
      (top_center_ix, true, height_segments * radial_segments),
    ] {
      generate_cap_indices(
        &mut indices,
        center_ix,
        radial_segments,
        ring_offset,
        reverse_winding,
      );
    }

    LinkedMesh::from_indexed_vertices(&vertices, &indices, None, None)
  }

  /// Generates a capsule mesh (cylinder with hemispherical caps).  The capsule is centered at the
  /// origin and extends along the Y axis.  `height` is the height of the cylindrical middle
  /// section; the total height is `height + 2 * radius`.
  ///
  /// Matches the semantics of Three.js `CapsuleGeometry`.
  pub fn new_capsule(
    radius: f32,
    height: f32,
    cap_segments: usize,
    radial_segments: usize,
    height_segments: usize,
  ) -> Self {
    let mut vertices = Vec::new();
    let mut indices = Vec::new();

    let half_height = height / 2.;
    // Number of vertex rings (excluding the two pole vertices)
    let total_rings = 2 * (cap_segments - 1) + height_segments + 1;

    for ring in 0..total_rings {
      let (y, r) = if ring < cap_segments - 1 {
        // Bottom hemisphere rings (excluding pole and equator)
        let angle = ((ring + 1) as f32 / cap_segments as f32) * (PI / 2.);
        (-half_height - radius * angle.cos(), radius * angle.sin())
      } else if ring < cap_segments - 1 + height_segments + 1 {
        // Cylinder section (including both equator rings)
        let t = (ring - (cap_segments - 1)) as f32 / height_segments as f32;
        (-half_height + t * height, radius)
      } else {
        // Top hemisphere rings (excluding pole and equator)
        let j = ring - (cap_segments - 1) - height_segments;
        let angle = (j as f32 / cap_segments as f32) * (PI / 2.);
        (half_height + radius * angle.sin(), radius * angle.cos())
      };

      for seg in 0..radial_segments {
        let theta = (seg as f32 / radial_segments as f32) * 2. * PI;
        vertices.push(Vec3::new(r * theta.cos(), y, r * theta.sin()));
      }
    }

    // Connect adjacent rings with quad strips
    generate_radial_side_indices(&mut indices, radial_segments, total_rings - 1);

    // Add pole vertices and connect with triangle fans
    let bottom_pole_ix = vertices.len() as u32;
    vertices.push(Vec3::new(0., -half_height - radius, 0.));
    let top_pole_ix = vertices.len() as u32;
    vertices.push(Vec3::new(0., half_height + radius, 0.));

    generate_cap_indices(&mut indices, bottom_pole_ix, radial_segments, 0, false);
    let last_ring_offset = (total_rings - 1) * radial_segments;
    generate_cap_indices(
      &mut indices,
      top_pole_ix,
      radial_segments,
      last_ring_offset,
      true,
    );

    LinkedMesh::from_indexed_vertices(&vertices, &indices, None, None)
  }

  /// Generates a cone mesh.  The base of the cone is centered at the origin, and the cone extends
  /// upwards along the positive Y axis.  The apex (point) of the cone is at the top.
  pub fn new_cone(
    radius: f32,
    height: f32,
    radial_segments: usize,
    height_segments: usize,
  ) -> Self {
    let mut vertices = Vec::new();
    let mut indices = Vec::new();

    for h in 0..height_segments {
      let h_norm = h as f32 / height_segments as f32;
      let y = h_norm * height;
      let r = radius * (1.0 - h_norm);
      for seg in 0..radial_segments {
        let theta = (seg as f32 / radial_segments as f32) * 2. * PI;
        let x = r * theta.cos();
        let z = r * theta.sin();
        vertices.push(Vec3::new(x, y, z));
      }
    }

    for h in 0..(height_segments - 1) {
      for r in 0..radial_segments {
        let v0 = h * radial_segments + r;
        let v1 = v0 + radial_segments;
        let next_r = (r + 1) % radial_segments;
        let v0_next = h * radial_segments + next_r;
        let v1_next = v0_next + radial_segments;

        indices.push(v0 as u32);
        indices.push(v1 as u32);
        indices.push(v0_next as u32);
        indices.push(v1 as u32);
        indices.push(v1_next as u32);
        indices.push(v0_next as u32);
      }
    }

    let apex_ix = vertices.len() as u32;
    vertices.push(Vec3::new(0., height, 0.));

    let top_ring_offset = (height_segments - 1) * radial_segments;
    generate_cap_indices(
      &mut indices,
      apex_ix,
      radial_segments,
      top_ring_offset,
      true,
    );

    let bottom_center_ix = vertices.len() as u32;
    vertices.push(Vec3::new(0., 0., 0.));
    generate_cap_indices(&mut indices, bottom_center_ix, radial_segments, 0, false);

    LinkedMesh::from_indexed_vertices(&vertices, &indices, None, None)
  }

  /// Generates a flat grid mesh in the XZ plane centered at the origin with the specified size and
  /// number of subdivisions.
  pub fn new_grid(
    size_x: f32,
    size_z: f32,
    subdivisions_x: usize,
    subdivisions_z: usize,
    flipped: bool,
  ) -> Self {
    let half_size_x = size_x / 2.;
    let half_size_z = size_z / 2.;

    let step_x = size_x / subdivisions_x as f32;
    let step_z = size_z / subdivisions_z as f32;

    let mut vertices = Vec::with_capacity((subdivisions_x + 1) * (subdivisions_z + 1));
    for z in 0..=subdivisions_z {
      for x in 0..=subdivisions_x {
        vertices.push(Vec3::new(
          -half_size_x + x as f32 * step_x,
          0.,
          -half_size_z + z as f32 * step_z,
        ));
      }
    }

    let mut indices = Vec::with_capacity(subdivisions_x * subdivisions_z * 6);
    for z in 0..subdivisions_z {
      for x in 0..subdivisions_x {
        let v0 = z * (subdivisions_x + 1) + x;
        let v1 = v0 + 1;
        let v2 = v0 + (subdivisions_x + 1);
        let v3 = v2 + 1;

        if flipped {
          indices.push(v0 as u32);
          indices.push(v1 as u32);
          indices.push(v2 as u32);

          indices.push(v1 as u32);
          indices.push(v3 as u32);
          indices.push(v2 as u32);
        } else {
          indices.push(v0 as u32);
          indices.push(v2 as u32);
          indices.push(v1 as u32);

          indices.push(v1 as u32);
          indices.push(v2 as u32);
          indices.push(v3 as u32);
        }
      }
    }

    LinkedMesh::from_indexed_vertices(&vertices, &indices, None, None)
  }

  #[cfg(not(target_arch = "wasm32"))]
  pub fn new_utah_teapot() -> Self {
    let verts = unsafe {
      std::slice::from_raw_parts(
        UTAH_TEAPOT_VERTICES.as_ptr() as *const Vec3,
        UTAH_TEAPOT_VERTICES.len(),
      )
    };
    let indices = UTAH_TEAPOT_INDICES;
    let indices = indices.iter().map(|&i| i as u32).collect::<Vec<_>>();

    Self::from_indexed_vertices(verts, &indices, None, None)
  }

  #[cfg(not(target_arch = "wasm32"))]
  pub fn new_stanford_bunny() -> Self {
    let verts = unsafe {
      std::slice::from_raw_parts(
        STANFORD_BUNNY_VERTICES.as_ptr() as *const Vec3,
        STANFORD_BUNNY_VERTICES.len(),
      )
    };
    let indices = STANFORD_BUNNY_INDICES;
    let indices = indices.iter().map(|&i| i as u32).collect::<Vec<_>>();

    Self::from_indexed_vertices(verts, &indices, None, None)
  }

  /// Actually flips the winding order of all faces in the mesh and updates all shading +
  /// displacement normals accordingly.
  pub fn flip_normals(&mut self) {
    for face in self.faces.values_mut() {
      face.vertices.swap(0, 2);
    }

    for n in self.shading_normals.values_mut() {
      *n = -*n;
    }
    for n in self.displacement_normals.values_mut() {
      *n = -*n;
    }
    for n in self.edge_displacement_normals.values_mut() {
      *n = -*n;
    }
    for ch in self.vertex_channels.values_mut() {
      ch.store.flip_all(ch.flip);
    }
  }

  /// Bakes affine transform `m` into the mesh in place: positions by `m`, shading/displacement
  /// normals by the inverse-transpose of its linear part, and `Direction` channels (tangents) by
  /// its linear part.  Orientation-reversing transforms (negative determinant) also reverse
  /// triangle winding and apply each passive channel's flip rule, keeping the surface
  /// consistently oriented.
  pub fn bake_transform(&mut self, m: Matrix4<f32>) {
    let linear = m.fixed_view::<3, 3>(0, 0).into_owned();
    let normal_mat = linear
      .try_inverse()
      .map(|inv| inv.transpose())
      .unwrap_or(linear);
    let renorm = |v: Vec3| {
      let n = v.norm();
      if n > 1e-12 {
        v / n
      } else {
        v
      }
    };

    for vtx in self.vertices.values_mut() {
      vtx.position = (m * vtx.position.push(1.)).xyz();
    }
    for n in self.shading_normals.values_mut() {
      *n = renorm(normal_mat * *n);
    }
    for n in self.displacement_normals.values_mut() {
      *n = renorm(normal_mat * *n);
    }
    for n in self.edge_displacement_normals.values_mut() {
      *n = renorm(normal_mat * *n);
    }
    for ch in self.vertex_channels.values_mut() {
      ch.transform_directions(&linear);
    }

    if linear.determinant() < 0. {
      for face in self.faces.values_mut() {
        face.vertices.swap(0, 2);
      }
      for ch in self.vertex_channels.values_mut() {
        ch.store.flip_all(ch.flip);
      }
    }
  }

  /// Mirrors the mesh in place across the plane `dot(normal, p) = offset` (`normal` must be unit).
  /// A Householder reflection baked via `bake_transform`, so positions, normals, and tangents are
  /// reflected and winding is reversed to keep the surface consistently oriented.
  pub fn reflect(&mut self, normal: Vec3, offset: f32) {
    let linear = Matrix3::identity() - 2. * normal * normal.transpose();
    let mut m = Matrix4::identity();
    m.fixed_view_mut::<3, 3>(0, 0).copy_from(&linear);
    m.fixed_view_mut::<3, 1>(0, 3)
      .copy_from(&(2. * offset * normal));
    self.bake_transform(m);
  }

  /// Appends `other`'s faces and the vertices they reference, with `xform` (other-local →
  /// self-local) baked per `bake_transform`. The channel set becomes the union, zero-filled where
  /// either side lacked a channel; flags are OR'd so authored seams survive.
  pub fn append(&mut self, other: &LinkedMesh<FaceData>, xform: &Mat4)
  where
    FaceData: Clone,
  {
    let other: Cow<LinkedMesh<FaceData>> = if *xform == Mat4::identity() {
      Cow::Borrowed(other)
    } else {
      let mut baked = other.clone();
      baked.bake_transform(*xform);
      Cow::Owned(baked)
    };

    let zero = [0f32; 4];
    for (name, ch) in &other.vertex_channels {
      if !self.vertex_channels.contains_key(name) {
        let mut mine = ch.empty_like();
        for k in self.vertices.keys() {
          mine.set(k, zero);
        }
        self.vertex_channels.insert(name.clone(), mine);
      }
    }

    let mut remap: FxHashMap<VertexKey, VertexKey> = FxHashMap::default();
    for face in other.faces.values() {
      let verts = face.vertices.map(|vk| {
        *remap.entry(vk).or_insert_with(|| {
          let nk = self
            .vertices
            .insert(Vertex::new(other.vertices[vk].position));
          if let Some(n) = other.shading_normals.get(vk) {
            self.shading_normals.insert(nk, *n);
          }
          if let Some(n) = other.displacement_normals.get(vk) {
            self.displacement_normals.insert(nk, *n);
          }
          nk
        })
      });
      self.add_face::<false>(verts, face.data.clone());
    }
    for (name, ch) in self.vertex_channels.iter_mut() {
      let src = other.vertex_channels.get(name);
      for (&ok, &nk) in &remap {
        ch.set(nk, src.and_then(|c| c.get(ok)).unwrap_or(zero));
      }
    }
    self.flags |= other.flags;
  }

  /// Builds `num_partitions` independent sub-meshes; face `fk` is placed into mesh
  /// `partition_of(fk)`.  Shared-vertex connectivity, per-vertex normals, and all passive vertex
  /// channels are preserved within each partition.
  pub fn extract_face_partitions(
    &self,
    num_partitions: usize,
    partition_of: impl Fn(FaceKey) -> usize,
  ) -> Vec<LinkedMesh<FaceData>>
  where
    FaceData: Clone,
  {
    let mut outs: Vec<LinkedMesh<FaceData>> = (0..num_partitions)
      .map(|_| {
        let mut out = LinkedMesh::new(0, 0, None);
        for (name, ch) in &self.vertex_channels {
          out.vertex_channels.insert(name.clone(), ch.empty_like());
        }
        out
      })
      .collect();
    let mut remaps: Vec<FxHashMap<VertexKey, VertexKey>> =
      vec![FxHashMap::default(); num_partitions];

    for (face_key, face) in self.faces.iter() {
      let idx = partition_of(face_key);
      let verts = face.vertices;
      let data = face.data.clone();
      let out = &mut outs[idx];
      let remap = &mut remaps[idx];
      let new_verts = verts.map(|vk| {
        *remap.entry(vk).or_insert_with(|| {
          let nk = out.vertices.insert(Vertex::new(self.vertices[vk].position));
          if let Some(n) = self.shading_normals.get(vk) {
            out.shading_normals.insert(nk, *n);
          }
          if let Some(n) = self.displacement_normals.get(vk) {
            out.displacement_normals.insert(nk, *n);
          }
          for (name, ch) in &self.vertex_channels {
            if let Some(v) = ch.get(vk) {
              out.vertex_channels.get_mut(name).unwrap().set(nk, v);
            }
          }
          nk
        })
      });
      out.add_face::<true>(new_verts, data);
    }

    outs
  }

  /// Splits the faces into independent sub-meshes keyed by `assign`'s returned index (`0..`),
  /// evaluated once per face with its three vertex positions in CCW winding order.  At least
  /// `min_partitions` meshes are always returned (trailing partitions may be empty).  Connectivity
  /// and attributes are preserved per `extract_face_partitions`.
  pub fn partition_faces<E>(
    &self,
    min_partitions: usize,
    mut assign: impl FnMut([Vec3; 3]) -> Result<usize, E>,
  ) -> Result<Vec<LinkedMesh<FaceData>>, E>
  where
    FaceData: Clone,
  {
    let mut idx_by_face: SecondaryMap<FaceKey, usize> = SecondaryMap::new();
    let mut num_partitions = min_partitions;
    for (face_key, face) in self.faces.iter() {
      let [a, b, c] = face.vertices;
      let idx = assign([
        self.vertices[a].position,
        self.vertices[b].position,
        self.vertices[c].position,
      ])?;
      num_partitions = num_partitions.max(idx + 1);
      idx_by_face.insert(face_key, idx);
    }
    Ok(self.extract_face_partitions(num_partitions, |fk| idx_by_face[fk]))
  }

  /// Computes the centroid of the mesh surface, weighted by triangle area.  Used as an
  /// approximation of the center of mass for non-2-manifold meshes.
  pub fn compute_surface_area_centroid(&self) -> Vec3 {
    let mut weighted_sum = Vec3::zeros();
    let mut total_area = 0.0f32;

    for face in self.faces.values() {
      let area = face.area(&self.vertices);
      let center = face.center(&self.vertices);
      weighted_sum += center * area;
      total_area += area;
    }

    if total_area > 0.0 {
      weighted_sum / total_area
    } else {
      Vec3::zeros()
    }
  }

  /// Computes the centroid of the volume enclosed by a 2-manifold mesh using the divergence
  /// theorem.
  ///
  /// Returns `None` if the mesh is not a closed 2-manifold.
  pub fn compute_volume_centroid(&self) -> Option<Vec3> {
    if self.check_is_manifold::<true>().is_err() {
      return None;
    }

    let mut total_volume = 0.0f32;
    let mut centroid_accum = Vec3::zeros();

    for face in self.faces.values() {
      let [p0, p1, p2] = face.vertex_positions(&self.vertices);

      // signed volume of tetrahedron formed by triangle and origin
      let cross = p1.cross(&p2);
      let signed_volume_6x = p0.dot(&cross);

      total_volume += signed_volume_6x;

      // Centroid of tetrahedron is at (p0 + p1 + p2 + origin) / 4 = (p0 + p1 + p2) / 4
      // Weighted by volume: contribution = V_tet * centroid_tet
      // = (signed_volume_6x / 6) * (p0 + p1 + p2) / 4
      // = signed_volume_6x * (p0 + p1 + p2) / 24
      centroid_accum += (p0 + p1 + p2) * signed_volume_6x;
    }

    // total_volume is 6x the actual volume, centroid_accum needs to be divided by 24
    // Final centroid = centroid_accum / 24 / (total_volume / 6) = centroid_accum / (4 *
    // total_volume)

    if total_volume.abs() < 1e-10 {
      // degenerate case - zero volume
      return None;
    }

    Some(centroid_accum / (4. * total_volume))
  }

  /// Computes the centroid of the mesh using the best available method based on mesh topology.
  ///
  /// For closed 2-manifold meshes, uses volume centroid (center of mass of enclosed solid).
  /// For other meshes, uses surface area weighted centroid.
  pub fn compute_centroid(&self) -> Vec3 {
    self
      .compute_volume_centroid()
      .unwrap_or_else(|| self.compute_surface_area_centroid())
  }

  /// Translates all vertices so that the mesh centroid is at the origin.
  ///
  /// For closed 2-manifold meshes, uses volume centroid (center of mass of enclosed solid).
  /// For other meshes, uses surface area weighted centroid.
  pub fn origin_to_geometry(&mut self) {
    let centroid = self.compute_centroid();
    for vtx in self.vertices.values_mut() {
      vtx.position -= centroid;
    }
  }

  /// Searches for the first pair of self-intersecting (non-adjacent) triangles in this mesh.
  ///
  /// Uses a BVH for broad-phase acceleration and the Guigue-Devillers algorithm for
  /// narrow-phase triangle-triangle intersection testing.  Face pairs sharing an edge are
  /// skipped.  Face pairs sharing a single vertex are tested, but point-only intersections
  /// at the shared vertex are filtered out.
  ///
  /// Returns `None` if no self-intersection is found.
  pub fn find_first_self_intersection(&self) -> Option<SelfIntersectionResult> {
    use crate::triangle_intersection::tri_tri_intersection;

    if self.faces.len() < 2 {
      return None;
    }
    let face_keys: Vec<FaceKey> = self.faces.keys().collect();
    let bvh = self.build_bvh(&Matrix4::identity());
    let mut result = None;
    bvh.for_each_pair(&bvh, 1e-6, |a, pos_a, b, pos_b| {
      let (fk_a, fk_b) = (face_keys[a as usize], face_keys[b as usize]);
      let (face_a, face_b) = (&self.faces[fk_a], &self.faces[fk_b]);
      let shared = face_a
        .vertices
        .iter()
        .filter(|v| face_b.vertices.contains(v))
        .count();
      // edge-adjacent faces always meet along their shared edge
      if shared >= 2 {
        return false;
      }
      let Some(isect) =
        tri_tri_intersection(pos_a[0], pos_a[1], pos_a[2], pos_b[0], pos_b[1], pos_b[2])
      else {
        return false;
      };
      result = Some(SelfIntersectionResult {
        face_a: fk_a,
        face_b: fk_b,
        positions_a: *pos_a,
        positions_b: *pos_b,
        point: isect.point,
        intersection_type: isect.intersection_type,
      });
      true
    });
    result
  }
}

/// Result of a self-intersection query on a mesh.
#[derive(Debug, Clone)]
pub struct SelfIntersectionResult {
  pub face_a: FaceKey,
  pub face_b: FaceKey,
  /// The three vertex positions of the first intersecting triangle
  pub positions_a: [Vec3; 3],
  /// The three vertex positions of the second intersecting triangle
  pub positions_b: [Vec3; 3],
  /// A representative point on or near the intersection
  pub point: Vec3,
  /// The type of intersection (segment or coplanar)
  pub intersection_type: crate::triangle_intersection::TriTriIntersectionType,
}

#[cfg(test)]
mod tests {
  use crate::slotmap_utils::fkey;

  use super::*;

  #[test]
  fn basic_edge_split() {
    let verts = [
      Vec3::new(-1., 0., 0.),
      Vec3::new(0., 1., 0.),
      Vec3::new(1., 0., 0.),
      Vec3::new(0., -1., 0.),
    ];

    // 1 -> 0 -> 3
    // 2 -> 1 -> 3
    let indices = [1, 0, 3, 2, 1, 3];

    let mut mesh: LinkedMesh<()> = LinkedMesh::from_indexed_vertices(&verts, &indices, None, None);
    let middle_edge_key = mesh
      .iter_edges()
      .find(|(_, e)| e.faces.len() == 2)
      .unwrap()
      .0;
    mesh.split_edge(
      middle_edge_key,
      EdgeSplitPos::middle(),
      DisplacementNormalMethod::EdgeNormal,
    );

    assert_eq!(mesh.iter_vertices().count(), 5);
    let center_vtx_key = mesh
      .iter_vertices()
      .find(|(_, v)| v.position == Vec3::zeros())
      .unwrap()
      .0;
    let center_vtx = &mesh.vertices[center_vtx_key];
    assert_eq!(center_vtx.edges.len(), 4);

    for &edge_key in &center_vtx.edges {
      if edge_key == middle_edge_key {
        continue;
      }

      let edge = &mesh.edges[edge_key];
      assert!(edge.faces.len() <= 2);

      for vtx_key in edge.vertices {
        if vtx_key == center_vtx_key {
          continue;
        }

        let vtx = &mesh.vertices[vtx_key];
        assert!(vtx.edges.len() <= 3);
      }
    }

    /// Returns true if the two triangles have the same vertices, regardless of
    /// order (as long as the winding order is consistent)
    fn tris_equal(t0: [[i32; 3]; 3], t1: [[i32; 3]; 3]) -> bool {
      return t0 == t1 || t0 == [t1[1], t1[2], t1[0]] || t0 == [t1[2], t1[0], t1[1]];
    }

    let expected_unindexed_verts = [
      [[0, 0, 0], [0, 1, 0], [-1, 0, 0]],
      [[1, 0, 0], [0, 1, 0], [0, 0, 0]],
      [[1, 0, 0], [0, 0, 0], [0, -1, 0]],
      [[0, -1, 0], [0, 0, 0], [-1, 0, 0]],
    ];
    let actual_indexed_mesh = mesh.to_raw_indexed(true, true, false);
    assert_eq!(actual_indexed_mesh.indices.len(), 12);
    let actual_verts = actual_indexed_mesh
      .indices
      .chunks(3)
      .map(|v| {
        let v0 = &actual_indexed_mesh.vertices[(v[0] * 3)..(v[0] * 3 + 3)];
        let v1 = &actual_indexed_mesh.vertices[(v[1] * 3)..(v[1] * 3 + 3)];
        let v2 = &actual_indexed_mesh.vertices[(v[2] * 3)..(v[2] * 3 + 3)];
        [
          [v0[0] as i32, v0[1] as i32, v0[2] as i32],
          [v1[0] as i32, v1[1] as i32, v1[2] as i32],
          [v2[0] as i32, v2[1] as i32, v2[2] as i32],
        ]
      })
      .collect::<Vec<_>>();

    assert_eq!(actual_verts.len(), expected_unindexed_verts.len());
    for expected_tri in &expected_unindexed_verts {
      assert!(
        actual_verts
          .iter()
          .any(|actual_tri| tris_equal(*actual_tri, *expected_tri)),
        "Expected triangle {:?} not found in actual triangles: {:?}",
        expected_tri,
        actual_verts
      );
    }
  }

  #[test]
  fn basic_smooth_fan() {
    // Smooth fan that consists of 3 edges and 2 faces, smooth all the way around
    let verts = [
      // vertex that we'll fan around
      Vec3::new(0., 0., 0.),
      Vec3::new(-1., 1., 0.),
      Vec3::new(0., 1., 0.),
      Vec3::new(1., 1., 0.),
    ];
    let indices = [0, 2, 1, 0, 3, 2];

    let mut mesh: LinkedMesh<()> = LinkedMesh::from_indexed_vertices(&verts, &indices, None, None);
    let vtx_key = vkey(1, 1);
    let mut visited_edges = bitarr![0; 1024];
    let visited_edges = &mut visited_edges[..mesh.vertices[vtx_key].edges.len()];
    assert_eq!(visited_edges.len(), 3);
    let mut smooth_fan_faces = SmallVec::<[_; 16]>::new();
    let mut visited_faces = SmallVec::<[_; 16]>::new();
    let (face_data, _) = mesh.precompute_face_shading::<false>();

    mesh.walk_one_smooth_fan(
      vtx_key,
      &face_data,
      visited_edges,
      &mut visited_faces,
      &mut smooth_fan_faces,
      &mut NormalAcc::new(),
    );

    assert!(visited_edges.all());
    let smooth_fan_faces: FxHashSet<_> = smooth_fan_faces.iter().copied().collect();
    assert_eq!(smooth_fan_faces.len(), 2);
  }

  fn test_walk_smooth_fan(
    mesh: &mut LinkedMesh,
    center_vtx_key: VertexKey,
  ) -> Vec<FxHashSet<FaceKey>> {
    let center_vtx_edge_count = mesh.vertices[center_vtx_key].edges.len();
    let mut visited_edges = bitarr![0; 1024];
    let visited_edges = &mut visited_edges[..center_vtx_edge_count];
    let mut smooth_fan_faces = SmallVec::<[_; 16]>::new();
    let mut visited_faces = SmallVec::<[_; 16]>::new();
    let (face_data, _) = mesh.precompute_face_shading::<false>();

    let mut smooth_fans = Vec::new();
    loop {
      mesh.walk_one_smooth_fan(
        center_vtx_key,
        &face_data,
        visited_edges,
        &mut visited_faces,
        &mut smooth_fan_faces,
        &mut NormalAcc::new(),
      );

      let uniq_smooth_fan_faces: FxHashSet<_> = smooth_fan_faces.iter().copied().collect();
      smooth_fans.push(uniq_smooth_fan_faces);
      smooth_fan_faces.clear();

      if visited_edges.all() {
        break;
      }
    }

    smooth_fans
  }

  #[test]
  fn two_triangles_joined_at_one_point() {
    let verts = [
      Vec3::new(0., 0., 0.),
      Vec3::new(-1., 1., 0.),
      Vec3::new(-1., 0., 0.),
      Vec3::new(1., 1., 0.),
      Vec3::new(1., 0., 0.),
    ];
    let indices = [0, 1, 2, 0, 4, 3];

    let mut mesh = LinkedMesh::from_indexed_vertices(&verts, &indices, None, None);
    let vtx_key = vkey(1, 1);

    let smooth_fans = test_walk_smooth_fan(&mut mesh, vtx_key);

    assert_eq!(smooth_fans.len(), 2);
    let expected_smooth_fans = [[fkey(1, 1)], [fkey(2, 1)]]
      .into_iter()
      .map(|fan| fan.into_iter().collect::<FxHashSet<_>>())
      .collect::<Vec<_>>();
    for smooth_fan in &smooth_fans {
      assert!(expected_smooth_fans.contains(smooth_fan));
    }

    // walking from any other vertex should only have a single smooth fan with a
    // single face
    for vtx_key in [
      unsafe { std::mem::transmute((1u32, 2u32)) },
      unsafe { std::mem::transmute((1u32, 3u32)) },
      unsafe { std::mem::transmute((1u32, 4u32)) },
      unsafe { std::mem::transmute((1u32, 5u32)) },
    ] {
      let mut visited_edges = bitarr![0; 1024];
      let visited_edges = &mut visited_edges[..mesh.vertices[vtx_key].edges.len()];
      let mut smooth_fan_faces = SmallVec::<[_; 16]>::new();
      let mut visited_faces = SmallVec::<[_; 16]>::new();
      let (face_data, _) = mesh.precompute_face_shading::<false>();

      let mut smooth_fans = Vec::new();
      loop {
        mesh.walk_one_smooth_fan(
          vtx_key,
          &face_data,
          visited_edges,
          &mut visited_faces,
          &mut smooth_fan_faces,
          &mut NormalAcc::new(),
        );

        let uniq_smooth_fan_faces: Vec<_> = smooth_fan_faces.iter().copied().collect();
        smooth_fans.push(uniq_smooth_fan_faces);
        smooth_fan_faces.clear();

        if visited_edges.all() {
          break;
        }
      }

      assert_eq!(smooth_fans.len(), 1);
      assert_eq!(smooth_fans[0].len(), 1);
    }
  }

  #[test]
  fn one_sharp_edge() {
    // Same as `basic_smooth_fan` but with a sharp edge between the two
    // triangles.  The result should be two separate smooth fans, each
    // containing one face.
    let verts = [
      // vertex that we'll fan around
      Vec3::new(0., 0., 0.),
      Vec3::new(-1., 1., 0.),
      Vec3::new(0., 1., 0.),
      Vec3::new(1., 1., 0.),
    ];
    let indices = [0, 2, 1, 0, 3, 2];

    let mut mesh = LinkedMesh::from_indexed_vertices(&verts, &indices, None, None);
    // mark the edge between the two triangles as sharp
    let sharp_edge_key = mesh
      .edges
      .iter_mut()
      .find_map(|(edge_key, edge)| {
        if edge.faces.len() == 2 {
          Some(edge_key)
        } else {
          None
        }
      })
      .unwrap();
    mesh.edges[sharp_edge_key].sharp = true;

    let vtx_key = vkey(1, 1);
    let smooth_fans = test_walk_smooth_fan(&mut mesh, vtx_key);

    assert_eq!(smooth_fans.len(), 2);
  }

  fn build_full_fan_mesh_with_sharp_edges(sharp_edges_vtx_coords: &[[[i32; 2]; 2]]) -> LinkedMesh {
    let verts = [[0, 0], [-1, 0], [0, 1], [1, 0], [0, -1]]
      .into_iter()
      .map(|[x, y]| Vec3::new(x as f32, y as f32, 0.))
      .collect::<Vec<_>>();
    let indices = [0, 2, 1, 0, 3, 2, 0, 1, 4, 0, 4, 3];

    let mut mesh = LinkedMesh::from_indexed_vertices(&verts, &indices, None, None);
    let sharp_edge_keys = mesh
      .edges
      .iter()
      .filter_map(|(edge_key, edge)| {
        let [v0_key, v1_key] = edge.vertices;
        let [v0_pos, v1_pos] = [
          [
            mesh.vertices[v0_key].position.x as i32,
            mesh.vertices[v0_key].position.y as i32,
          ],
          [
            mesh.vertices[v1_key].position.x as i32,
            mesh.vertices[v1_key].position.y as i32,
          ],
        ];

        for &[sharp_edge_v0_pos, sharp_edge_v1_pos] in sharp_edges_vtx_coords {
          if (sharp_edge_v0_pos == v0_pos && sharp_edge_v1_pos == v1_pos)
            || (sharp_edge_v0_pos == v1_pos && sharp_edge_v1_pos == v0_pos)
          {
            return Some(edge_key);
          }
        }

        None
      })
      .collect::<Vec<_>>();
    assert_eq!(sharp_edge_keys.len(), sharp_edges_vtx_coords.len());

    for sharp_edge_key in sharp_edge_keys {
      mesh.edges.get_mut(sharp_edge_key).unwrap().sharp = true;
    }

    mesh
  }

  #[test]
  fn full_fan_two_sharp_edges() {
    let mut mesh = build_full_fan_mesh_with_sharp_edges(&[[[0, 0], [0, 1]], [[0, 0], [0, -1]]]);

    let center_vtx_key = vkey(1, 1);
    let smooth_fans = test_walk_smooth_fan(&mut mesh, center_vtx_key);

    // faces 1 + 3 should be together, and 2 + 4 should be together
    assert_eq!(smooth_fans.len(), 2);
    let expected_smooth_fans = [[fkey(1, 1), fkey(3, 1)], [fkey(2, 1), fkey(4, 1)]]
      .into_iter()
      .map(|v| v.into_iter().collect::<FxHashSet<_>>())
      .collect::<Vec<_>>();

    for smooth_fan in &smooth_fans {
      assert!(expected_smooth_fans
        .iter()
        .any(|expected| expected == smooth_fan));
    }
  }

  /// A full fan with one sharp edges will actually result in only a single
  /// smooth fan because there's no way to split the vertex without breaking
  /// topology
  #[test]
  fn full_fan_one_sharp_edge() {
    let mut mesh = build_full_fan_mesh_with_sharp_edges(&[[[0, 0], [0, 1]]]);

    let center_vtx_key = vkey(1, 1);
    let smooth_fans = test_walk_smooth_fan(&mut mesh, center_vtx_key);

    assert_eq!(smooth_fans.len(), 1);
    let expected_smooth_fan = [fkey(1, 1), fkey(2, 1), fkey(3, 1), fkey(4, 1)]
      .into_iter()
      .collect::<FxHashSet<_>>();

    assert_eq!(smooth_fans[0], expected_smooth_fan);
  }

  #[test]
  fn cube_two_manifold() {
    let mesh: LinkedMesh<()> = LinkedMesh::new_box(1., 1., 1.);

    mesh
      .check_is_manifold::<true>()
      .expect("basic cube mesh should be two-manifold");
  }

  #[test]
  fn subdivide_by_plane_manifold() {
    let mut mesh: LinkedMesh<()> = LinkedMesh::new_box(1., 1., 1.);
    let plane = Plane {
      normal: Vec3::new(1., 1., 0.).normalize(),
      w: 0.,
    };

    mesh.subdivide_by_plane(&plane);

    mesh
      .check_is_manifold::<true>()
      .expect("Mesh should remain manifold after plane subdivision");

    assert!(mesh.faces.len() > 12);
  }

  #[test]
  fn skewer_cube_through_origin_hits_diagonal_edges() {
    let mut mesh: LinkedMesh<()> = LinkedMesh::new_box(1., 1., 1.);
    mesh.subdivide_by_line(Vec3::zeros(), Vec3::new(0., 1., 0.), 1e-4);
    mesh
      .check_is_manifold::<true>()
      .expect("Mesh should remain manifold after edge-hit skewer");
    // 12 base tris, top+bottom diagonal edge splits each turn 2 tris into 4 (+2 each, +4 total)
    assert_eq!(mesh.faces.len(), 16);
  }

  #[test]
  fn skewer_cube_through_corners_no_change() {
    let mut mesh: LinkedMesh<()> = LinkedMesh::new_box(1., 1., 1.);
    let face_count_before = mesh.faces.len();
    mesh.subdivide_by_line(Vec3::zeros(), Vec3::new(1., 1., 1.).normalize(), 1e-4);
    mesh
      .check_is_manifold::<true>()
      .expect("Mesh should remain manifold after vertex-hit skewer");
    assert_eq!(mesh.faces.len(), face_count_before);
  }

  #[test]
  fn skewer_cube_interior_hits() {
    let mut mesh: LinkedMesh<()> = LinkedMesh::new_box(1., 1., 1.);
    // (0.2, 0, 0.2) is well inside one of the top/bottom face tris, not on the diagonal (z = -x)
    mesh.subdivide_by_line(Vec3::new(0.2, 0., 0.2), Vec3::new(0., 1., 0.), 1e-4);
    mesh
      .check_is_manifold::<true>()
      .expect("Mesh should remain manifold after interior-hit skewer");
    // Two interior fan splits: each removes 1 tri, adds 3 (+2 each, +4 total)
    assert_eq!(mesh.faces.len(), 16);
  }

  #[test]
  fn volume_centroid_unit_cube() {
    let mesh: LinkedMesh<()> = LinkedMesh::new_box(2., 2., 2.);

    let centroid = mesh.compute_volume_centroid();
    assert!(centroid.is_some(), "Unit cube should be 2-manifold");

    let centroid = centroid.unwrap();
    assert!(
      centroid.magnitude() < 1e-5,
      "Centroid of centered cube should be at origin, got {centroid:?}",
    );
  }

  #[test]
  fn volume_centroid_offset_cube() {
    let mut mesh: LinkedMesh<()> = LinkedMesh::new_box(2., 2., 2.);

    let offset = Vec3::new(1., 2., 3.);
    for vtx in mesh.vertices.values_mut() {
      vtx.position += offset;
    }

    let centroid = mesh.compute_volume_centroid().unwrap();
    let diff = (centroid - offset).magnitude();
    assert!(
      diff < 1e-5,
      "Centroid should be at offset {offset:?}, got {centroid:?}",
    );
  }

  #[test]
  fn surface_area_centroid_unit_cube() {
    let mesh: LinkedMesh<()> = LinkedMesh::new_box(2., 2., 2.);

    let centroid = mesh.compute_surface_area_centroid();
    assert!(
      centroid.magnitude() < 1e-5,
      "Surface centroid of centered cube should be at origin, got {centroid:?}",
    );
  }

  #[test]
  fn origin_to_geometry_centers_mesh() {
    let mut mesh: LinkedMesh<()> = LinkedMesh::new_box(2., 2., 2.);

    let offset = Vec3::new(5., -3., 7.);
    for vtx in mesh.vertices.values_mut() {
      vtx.position += offset;
    }

    let centroid_before = mesh.compute_centroid();
    assert!(
      (centroid_before - offset).magnitude() < 1e-5,
      "Mesh should be offset before origin_to_geometry"
    );

    mesh.origin_to_geometry();

    let centroid_after = mesh.compute_centroid();
    assert!(
      centroid_after.magnitude() < 1e-5,
      "Mesh should be centered after origin_to_geometry, got {centroid_after:?}",
    );
  }

  #[test]
  fn capsule_is_manifold_and_centered() {
    let mesh: LinkedMesh<()> = LinkedMesh::new_capsule(1., 2., 4, 8, 1);

    mesh
      .check_is_manifold::<true>()
      .expect("Capsule should be 2-manifold");

    let centroid = mesh.compute_volume_centroid();
    assert!(
      centroid.is_some(),
      "Capsule should be 2-manifold (volume centroid)"
    );

    let centroid = centroid.unwrap();
    assert!(
      centroid.magnitude() < 1e-4,
      "Centroid of centered capsule should be at origin, got {centroid:?}",
    );
  }

  #[test]
  fn capsule_minimal_cap_segments() {
    let mesh: LinkedMesh<()> = LinkedMesh::new_capsule(0.5, 1., 1, 6, 1);

    mesh
      .check_is_manifold::<true>()
      .expect("Capsule with cap_segments=1 should be 2-manifold");
  }
}