nightshade-renderer 0.57.0

use crate::geometry::Mesh;
use crate::meshlet::asset::{
    BvhNode, Meshlet, MeshletAabb, MeshletAabbErrorOffset, MeshletBoundingSphere, MeshletCullData,
    MeshletMesh,
};
use bitvec::order::Lsb0;
use bitvec::vec::BitVec;
use bitvec::view::BitView;
use itertools::Itertools;
use meshopt::ffi::meshopt_Meshlet;
use meshopt::{
    Meshlets, SimplifyOptions, VertexDataAdapter, build_meshlets, compute_meshlet_bounds,
    generate_position_remap, simplify_with_attributes_and_locks,
};
use metis::Graph;
use metis::option::Opt;
use nalgebra_glm::{Vec2, Vec3};
use std::collections::HashMap;
use std::ops::Range;

const TARGET_MESHLETS_PER_GROUP: usize = 8;
const SIMPLIFICATION_FAILURE_PERCENTAGE: f32 = 0.60;
const CENTIMETERS_PER_METER: f32 = 100.0;
const PACKED_VERTEX_STRIDE: usize = 32;

/// Default vertex position quantization factor for [`from_mesh`].
///
/// Snaps vertices to the nearest 1/16th of a centimeter (1/2^4).
pub const MESHLET_DEFAULT_VERTEX_POSITION_QUANTIZATION_FACTOR: u8 = 4;

/// Processes a [`Mesh`] into a [`MeshletMesh`] hierarchy of triangle clusters and
/// a BVH8 level-of-detail DAG. This is slow and meant to run ahead of time.
///
/// `vertex_position_quantization_factor` controls how much precision to keep when
/// quantizing vertex positions: vertices snap to the nearest (1/2^x)th of a
/// centimeter, where x is the factor. Use
/// [`MESHLET_DEFAULT_VERTEX_POSITION_QUANTIZATION_FACTOR`] as a starting point.
pub fn from_mesh(
    mesh: &Mesh,
    vertex_position_quantization_factor: u8,
) -> Result<MeshletMesh, MeshToMeshletMeshConversionError> {
    validate_input_mesh(mesh)?;

    let mut vertex_buffer = Vec::with_capacity(mesh.vertices.len() * PACKED_VERTEX_STRIDE);
    for vertex in &mesh.vertices {
        vertex_buffer.extend_from_slice(bytemuck::bytes_of(&vertex.position));
        vertex_buffer.extend_from_slice(bytemuck::bytes_of(&vertex.normal));
        vertex_buffer.extend_from_slice(bytemuck::bytes_of(&vertex.tex_coords));
    }

    let vertices = VertexDataAdapter::new(&vertex_buffer, PACKED_VERTEX_STRIDE, 0)?;
    let vertex_normal_attributes: Vec<f32> = mesh
        .vertices
        .iter()
        .flat_map(|vertex| vertex.normal)
        .collect();

    let position_only_vertex_remap = generate_position_remap(&vertices);
    let indices = drop_degenerate_triangles(&mesh.indices, &position_only_vertex_remap);
    if indices.is_empty() {
        return Err(MeshToMeshletMeshConversionError::MeshMissingIndices);
    }

    let (mut meshlets, mut cull_data) =
        compute_meshlets(&indices, &vertices, &position_only_vertex_remap, None)?;

    let mut vertex_locks = vec![false; vertices.vertex_count];

    let mut bvh_builder = BvhBuilder::default();
    let mut all_groups = Vec::new();
    let mut simplification_queue: Vec<u32> = (0..meshlets.len() as u32).collect();
    let mut stuck = Vec::new();
    while !simplification_queue.is_empty() {
        let connected_meshlets_per_meshlet = find_connected_meshlets(
            &simplification_queue,
            &meshlets,
            &position_only_vertex_remap,
        );

        let groups = group_meshlets(
            &simplification_queue,
            &cull_data,
            &connected_meshlets_per_meshlet,
        )?;
        simplification_queue.clear();

        lock_group_borders(
            &mut vertex_locks,
            &groups,
            &meshlets,
            &position_only_vertex_remap,
        );

        let mut simplified = Vec::with_capacity(groups.len());
        for mut group in groups {
            if group.meshlets.len() == 1 {
                simplified.push(Err(group));
                continue;
            }

            let Some((simplified_group_indices, mut group_error)) = simplify_meshlet_group(
                &group,
                &meshlets,
                &vertices,
                &vertex_normal_attributes,
                &vertex_locks,
            ) else {
                simplified.push(Err(group));
                continue;
            };

            for &meshlet_id in group.meshlets.iter() {
                group_error = group_error.max(cull_data[meshlet_id as usize].error);
            }
            group.parent_error = group_error;

            let new_meshlets = compute_meshlets(
                &simplified_group_indices,
                &vertices,
                &position_only_vertex_remap,
                Some((group.lod_bounds, group.parent_error)),
            )?;

            simplified.push(Ok((group, new_meshlets)));
        }

        let first_group = all_groups.len() as u32;
        let mut passed_triangles = 0;
        let mut stuck_triangles = 0;
        for entry in simplified {
            match entry {
                Ok((group, (new_meshlets, new_cull_data))) => {
                    let start = meshlets.len();
                    merge_meshlets(&mut meshlets, new_meshlets);
                    cull_data.extend(new_cull_data);
                    let end = meshlets.len();
                    let new_meshlet_ids = start as u32..end as u32;

                    passed_triangles += triangles_in_meshlets(&meshlets, new_meshlet_ids.clone());
                    simplification_queue.extend(new_meshlet_ids);
                    all_groups.push(group);
                }
                Err(group) => {
                    stuck_triangles +=
                        triangles_in_meshlets(&meshlets, group.meshlets.iter().copied());
                    stuck.push(group);
                }
            }
        }

        if passed_triangles > stuck_triangles / 3 {
            simplification_queue.extend(stuck.drain(..).flat_map(|group| group.meshlets));
        }

        bvh_builder.add_lod(first_group, &all_groups);
    }

    if !stuck.is_empty() {
        let first_group = all_groups.len() as u32;
        all_groups.extend(stuck);
        bvh_builder.add_lod(first_group, &all_groups);
    }

    let (bvh, aabb, bvh_depth) = bvh_builder.build(&mut meshlets, all_groups, &mut cull_data);

    let mut accumulator = VertexDataAccumulator {
        vertex_positions: BitVec::<u32, Lsb0>::new(),
        vertex_normals: Vec::new(),
        vertex_uvs: Vec::new(),
        meshlets: Vec::with_capacity(meshlets.len()),
    };
    for index in 0..meshlets.meshlets.len() {
        let meshlet = meshlets.meshlets[index];
        accumulator.append_meshlet(
            &meshlet,
            meshlets.get(index).vertices,
            &vertex_buffer,
            PACKED_VERTEX_STRIDE,
            vertex_position_quantization_factor,
        );
    }
    accumulator.vertex_positions.set_uninitialized(false);

    Ok(MeshletMesh {
        vertex_positions: accumulator.vertex_positions.into_vec().into(),
        vertex_normals: accumulator.vertex_normals.into(),
        vertex_uvs: accumulator.vertex_uvs.into(),
        indices: meshlets.triangles.into(),
        bvh: bvh.into(),
        meshlets: accumulator.meshlets.into(),
        meshlet_cull_data: cull_data
            .into_iter()
            .map(|cull_data| MeshletCullData {
                aabb: aabb_to_meshlet(cull_data.aabb, cull_data.error, 0),
                lod_group_sphere: sphere_to_meshlet(cull_data.lod_group_sphere),
            })
            .collect(),
        aabb,
        bvh_depth,
    })
}

fn validate_input_mesh(mesh: &Mesh) -> Result<(), MeshToMeshletMeshConversionError> {
    if mesh.indices.is_empty() {
        return Err(MeshToMeshletMeshConversionError::MeshMissingIndices);
    }
    if !mesh.indices.len().is_multiple_of(3) {
        return Err(MeshToMeshletMeshConversionError::WrongMeshPrimitiveTopology);
    }
    // Everything downstream indexes the vertex list by these without checking,
    // so an index naming a vertex that is not there has to be caught before it
    // becomes a panic somewhere less obvious.
    if let Some(vertex) = mesh
        .indices
        .iter()
        .copied()
        .find(|index| *index as usize >= mesh.vertices.len())
    {
        return Err(MeshToMeshletMeshConversionError::IndexOutOfBounds {
            vertex,
            vertex_count: mesh.vertices.len(),
        });
    }
    Ok(())
}

/// Drops triangles that enclose no area, judged by position rather than by
/// index.
///
/// Two vertices can be distinct entries and still sit in the same place: a
/// sphere's pole is one vertex per meridian, all coincident, and the triangles
/// spanning them are slivers with no area. Meshopt answers a cluster of those
/// with zeroed bounds rather than an error, which reads as a bounding sphere at
/// the origin with no radius, and the level of detail hierarchy built on top of
/// it violates its own monotonicity. Art exported from a modelling tool carries
/// these routinely, so they are dropped rather than diagnosed: they cover no
/// pixels, and nothing downstream misses them.
fn drop_degenerate_triangles(indices: &[u32], position_only_vertex_remap: &[u32]) -> Vec<u32> {
    let mut kept = Vec::with_capacity(indices.len());
    for triangle in indices.chunks_exact(3) {
        let a = position_only_vertex_remap[triangle[0] as usize];
        let b = position_only_vertex_remap[triangle[1] as usize];
        let c = position_only_vertex_remap[triangle[2] as usize];
        if a == b || b == c || a == c {
            continue;
        }
        kept.extend_from_slice(triangle);
    }
    kept
}

fn triangles_in_meshlets(meshlets: &Meshlets, ids: impl IntoIterator<Item = u32>) -> u32 {
    ids.into_iter()
        .map(|id| meshlets.get(id as usize).triangles.len() as u32 / 3)
        .sum()
}

fn compute_meshlets(
    indices: &[u32],
    vertices: &VertexDataAdapter<'_>,
    position_only_vertex_remap: &[u32],
    previous_lod_data: Option<(BoundingSphere, f32)>,
) -> Result<(Meshlets, Vec<TempMeshletCullData>), MeshToMeshletMeshConversionError> {
    let mut vertices_to_triangles = vec![Vec::new(); position_only_vertex_remap.len()];
    for (index, vertex_index) in indices.iter().enumerate() {
        let vertex_id = position_only_vertex_remap[*vertex_index as usize];
        vertices_to_triangles[vertex_id as usize].push(index / 3);
    }

    let mut triangle_pair_to_shared_vertex_count: HashMap<(usize, usize), usize> = HashMap::new();
    for triangle_ids in vertices_to_triangles {
        for (first_triangle, second_triangle) in triangle_ids.into_iter().tuple_combinations() {
            let count = triangle_pair_to_shared_vertex_count
                .entry((
                    first_triangle.min(second_triangle),
                    first_triangle.max(second_triangle),
                ))
                .or_insert(0);
            *count += 1;
        }
    }

    let triangle_count = indices.len() / 3;
    let mut connected_triangles_per_triangle = vec![Vec::new(); triangle_count];
    for ((first_triangle, second_triangle), shared_vertex_count) in
        triangle_pair_to_shared_vertex_count
    {
        connected_triangles_per_triangle[first_triangle]
            .push((second_triangle, shared_vertex_count));
        connected_triangles_per_triangle[second_triangle]
            .push((first_triangle, shared_vertex_count));
    }
    for list in connected_triangles_per_triangle.iter_mut() {
        list.sort_unstable();
    }

    let mut xadj = Vec::with_capacity(triangle_count + 1);
    let mut adjncy = Vec::new();
    let mut adjwgt = Vec::new();
    for connected_triangles in connected_triangles_per_triangle {
        xadj.push(adjncy.len() as metis::Idx);
        for (connected_triangle, shared_vertex_count) in connected_triangles {
            adjncy.push(connected_triangle as metis::Idx);
            adjwgt.push(shared_vertex_count as metis::Idx);
        }
    }
    xadj.push(adjncy.len() as metis::Idx);

    let mut options = [-1; metis::NOPTIONS];
    options[metis::option::Seed::INDEX] = 17;
    options[metis::option::UFactor::INDEX] = 1;

    let mut meshlet_per_triangle = vec![0; triangle_count];
    let partition_count = triangle_count.div_ceil(126);
    Graph::new(1, partition_count as metis::Idx, &xadj, &adjncy)?
        .set_options(&options)
        .set_adjwgt(&adjwgt)
        .part_recursive(&mut meshlet_per_triangle)?;

    let mut indices_per_meshlet = vec![Vec::new(); partition_count];
    for (triangle_id, meshlet) in meshlet_per_triangle.into_iter().enumerate() {
        let base_index = triangle_id * 3;
        indices_per_meshlet[meshlet as usize]
            .extend_from_slice(&indices[base_index..base_index + 3]);
    }

    let mut meshlets = Meshlets {
        meshlets: Vec::new(),
        vertices: Vec::new(),
        triangles: Vec::new(),
    };
    let mut cull_data = Vec::new();
    for meshlet_indices in &indices_per_meshlet {
        let built = build_meshlets(meshlet_indices, vertices, 256, 128, 0.0);
        for meshlet in built.iter() {
            let positions: Vec<Vec3> = meshlet
                .vertices
                .iter()
                .map(|&vertex_id| read_vertex_position(vertices, vertex_id))
                .collect();
            let aabb = Aabb::from_points(positions.iter().copied());

            let (lod_group_sphere, error) = match previous_lod_data {
                Some(data) => data,
                None => {
                    let bounds = compute_meshlet_bounds(meshlet, vertices);
                    let sphere = if bounds.radius > 0.0 {
                        BoundingSphere::new(array_to_vec3(bounds.center), bounds.radius)
                    } else {
                        enclosing_sphere(aabb.center(), &positions)
                    };
                    (sphere, 0.0)
                }
            };

            cull_data.push(TempMeshletCullData {
                aabb,
                lod_group_sphere,
                error,
            });
        }
        merge_meshlets(&mut meshlets, built);
    }

    Ok((meshlets, cull_data))
}

fn find_connected_meshlets(
    simplification_queue: &[u32],
    meshlets: &Meshlets,
    position_only_vertex_remap: &[u32],
) -> Vec<Vec<(usize, usize)>> {
    let mut vertices_to_meshlets = vec![Vec::new(); position_only_vertex_remap.len()];
    for (local_index, &meshlet_id) in simplification_queue.iter().enumerate() {
        let meshlet = meshlets.get(meshlet_id as usize);
        for index in meshlet.triangles {
            let vertex_id = position_only_vertex_remap[meshlet.vertices[*index as usize] as usize];
            let vertex_to_meshlets = &mut vertices_to_meshlets[vertex_id as usize];
            if vertex_to_meshlets.last() != Some(&local_index) {
                vertex_to_meshlets.push(local_index);
            }
        }
    }

    let mut meshlet_pair_to_shared_vertex_count: HashMap<(usize, usize), usize> = HashMap::new();
    for meshlet_ids in vertices_to_meshlets {
        for (first_meshlet, second_meshlet) in meshlet_ids.into_iter().tuple_combinations() {
            let count = meshlet_pair_to_shared_vertex_count
                .entry((
                    first_meshlet.min(second_meshlet),
                    first_meshlet.max(second_meshlet),
                ))
                .or_insert(0);
            *count += 1;
        }
    }

    let mut connected_meshlets_per_meshlet = vec![Vec::new(); simplification_queue.len()];
    for ((first_meshlet, second_meshlet), shared_vertex_count) in
        meshlet_pair_to_shared_vertex_count
    {
        connected_meshlets_per_meshlet[first_meshlet].push((second_meshlet, shared_vertex_count));
        connected_meshlets_per_meshlet[second_meshlet].push((first_meshlet, shared_vertex_count));
    }
    for list in connected_meshlets_per_meshlet.iter_mut() {
        list.sort_unstable();
    }

    connected_meshlets_per_meshlet
}

fn group_meshlets(
    simplification_queue: &[u32],
    meshlet_cull_data: &[TempMeshletCullData],
    connected_meshlets_per_meshlet: &[Vec<(usize, usize)>],
) -> Result<Vec<TempMeshletGroup>, MeshToMeshletMeshConversionError> {
    let mut xadj = Vec::with_capacity(simplification_queue.len() + 1);
    let mut adjncy = Vec::new();
    let mut adjwgt = Vec::new();
    for connected_meshlets in connected_meshlets_per_meshlet {
        xadj.push(adjncy.len() as metis::Idx);
        for (connected_meshlet, shared_vertex_count) in connected_meshlets {
            adjncy.push(*connected_meshlet as metis::Idx);
            adjwgt.push(*shared_vertex_count as metis::Idx);
        }
    }
    xadj.push(adjncy.len() as metis::Idx);

    let mut options = [-1; metis::NOPTIONS];
    options[metis::option::Seed::INDEX] = 17;
    options[metis::option::UFactor::INDEX] = 200;

    let mut group_per_meshlet = vec![0; simplification_queue.len()];
    let partition_count = simplification_queue
        .len()
        .div_ceil(TARGET_MESHLETS_PER_GROUP);
    Graph::new(1, partition_count as metis::Idx, &xadj, &adjncy)?
        .set_options(&options)
        .set_adjwgt(&adjwgt)
        .part_recursive(&mut group_per_meshlet)?;

    let mut groups = vec![TempMeshletGroup::default(); partition_count];
    for (local_index, meshlet_group) in group_per_meshlet.into_iter().enumerate() {
        let group = &mut groups[meshlet_group as usize];
        let meshlet_id = simplification_queue[local_index];

        group.meshlets.push(meshlet_id);
        let data = &meshlet_cull_data[meshlet_id as usize];
        group.aabb = group.aabb.merge(&data.aabb);
        group.lod_bounds = merge_spheres(group.lod_bounds, data.lod_group_sphere);
    }

    Ok(groups)
}

fn lock_group_borders(
    vertex_locks: &mut [bool],
    groups: &[TempMeshletGroup],
    meshlets: &Meshlets,
    position_only_vertex_remap: &[u32],
) {
    let mut position_only_locks = vec![-1_i32; position_only_vertex_remap.len()];

    for (group_id, group) in groups.iter().enumerate() {
        for &meshlet_id in group.meshlets.iter() {
            let meshlet = meshlets.get(meshlet_id as usize);
            for index in meshlet.triangles {
                let vertex_id =
                    position_only_vertex_remap[meshlet.vertices[*index as usize] as usize] as usize;

                if position_only_locks[vertex_id] == -1
                    || position_only_locks[vertex_id] == group_id as i32
                {
                    position_only_locks[vertex_id] = group_id as i32;
                } else {
                    position_only_locks[vertex_id] = -2;
                }
            }
        }
    }

    for (lock, &remap) in vertex_locks
        .iter_mut()
        .zip(position_only_vertex_remap.iter())
    {
        *lock = position_only_locks[remap as usize] == -2;
    }
}

fn simplify_meshlet_group(
    group: &TempMeshletGroup,
    meshlets: &Meshlets,
    vertices: &VertexDataAdapter<'_>,
    vertex_normal_attributes: &[f32],
    vertex_locks: &[bool],
) -> Option<(Vec<u32>, f32)> {
    let group_indices: Vec<u32> = group
        .meshlets
        .iter()
        .flat_map(|&meshlet_id| {
            let meshlet = meshlets.get(meshlet_id as usize);
            meshlet
                .triangles
                .iter()
                .map(move |&meshlet_index| meshlet.vertices[meshlet_index as usize])
        })
        .collect();

    let mut error = 0.0;
    let simplified_group_indices = simplify_with_attributes_and_locks(
        &group_indices,
        vertices,
        vertex_normal_attributes,
        &[0.5; 3],
        std::mem::size_of::<[f32; 3]>(),
        vertex_locks,
        group_indices.len() / 2,
        f32::MAX,
        SimplifyOptions::Sparse | SimplifyOptions::ErrorAbsolute,
        Some(&mut error),
    );

    if simplified_group_indices.len() as f32 / group_indices.len() as f32
        > SIMPLIFICATION_FAILURE_PERCENTAGE
    {
        return None;
    }

    Some((simplified_group_indices, error))
}

fn merge_meshlets(meshlets: &mut Meshlets, merge: Meshlets) {
    let vertex_offset = meshlets.vertices.len() as u32;
    let triangle_offset = meshlets.triangles.len() as u32;
    meshlets.vertices.extend_from_slice(&merge.vertices);
    meshlets.triangles.extend_from_slice(&merge.triangles);
    meshlets
        .meshlets
        .extend(merge.meshlets.into_iter().map(|mut meshlet| {
            meshlet.vertex_offset += vertex_offset;
            meshlet.triangle_offset += triangle_offset;
            meshlet
        }));
}

fn read_vertex_position(vertices: &VertexDataAdapter<'_>, vertex_id: u32) -> Vec3 {
    let bytes = *vertices.reader.get_ref();
    let start = vertices.position_offset + vertex_id as usize * vertices.vertex_stride;
    read_vec3(&bytes[start..start + 12])
}

fn read_vec3(bytes: &[u8]) -> Vec3 {
    Vec3::new(
        f32::from_ne_bytes(bytes[0..4].try_into().unwrap()),
        f32::from_ne_bytes(bytes[4..8].try_into().unwrap()),
        f32::from_ne_bytes(bytes[8..12].try_into().unwrap()),
    )
}

fn read_vec2(bytes: &[u8]) -> Vec2 {
    Vec2::new(
        f32::from_ne_bytes(bytes[0..4].try_into().unwrap()),
        f32::from_ne_bytes(bytes[4..8].try_into().unwrap()),
    )
}

fn array_to_vec3(array: [f32; 3]) -> Vec3 {
    Vec3::new(array[0], array[1], array[2])
}

fn vec3_to_array(vector: Vec3) -> [f32; 3] {
    [vector.x, vector.y, vector.z]
}

fn octahedral_encode(normal: Vec3) -> Vec2 {
    let normalized = normal / (normal.x.abs() + normal.y.abs() + normal.z.abs());
    let wrapped = Vec2::new(
        (1.0 - normalized.y.abs()) * if normalized.x >= 0.0 { 1.0 } else { -1.0 },
        (1.0 - normalized.x.abs()) * if normalized.y >= 0.0 { 1.0 } else { -1.0 },
    );
    if normalized.z >= 0.0 {
        Vec2::new(normalized.x, normalized.y)
    } else {
        wrapped
    }
}

fn pack2x16snorm(value: Vec2) -> u32 {
    let x = (value.x.clamp(-1.0, 1.0) * 32767.0 + 0.5).floor() as i16;
    let y = (value.y.clamp(-1.0, 1.0) * 32767.0 + 0.5).floor() as i16;
    (x as u16 as u32) | ((y as u16 as u32) << 16)
}

/// Bounds `points` with a sphere around `center`.
///
/// meshopt skips zero-area triangles when bounding a cluster and returns an
/// all-zero sphere for a cluster that has no non-degenerate triangles at all.
/// That sphere sits at the origin rather than on the geometry, so it would
/// break the nesting the lod hierarchy relies on. This bounds the meshlet's own
/// vertices instead, which stays on the geometry.
fn enclosing_sphere(center: Vec3, points: &[Vec3]) -> BoundingSphere {
    let radius = points
        .iter()
        .map(|point| (point - center).norm())
        .fold(0.0_f32, f32::max);
    BoundingSphere::new(center, radius)
}

fn merge_spheres(first: BoundingSphere, second: BoundingSphere) -> BoundingSphere {
    let smaller_radius = first.radius().min(second.radius());
    let larger_radius = first.radius().max(second.radius());
    let distance = first.center_distance(&second);
    if distance + smaller_radius <= larger_radius || smaller_radius <= 0.0 || distance <= 0.0 {
        if first.radius() > second.radius() {
            first
        } else {
            second
        }
    } else {
        let radius = (smaller_radius + larger_radius + distance) / 2.0;
        let center = (first.center()
            + second.center()
            + (first.radius() - second.radius()) * (first.center() - second.center()) / distance)
            / 2.0;
        BoundingSphere::new(center, radius)
    }
}

fn aabb_to_meshlet(aabb: Aabb, error: f32, child_offset: u32) -> MeshletAabbErrorOffset {
    MeshletAabbErrorOffset {
        center: vec3_to_array(aabb.center()),
        error,
        half_extent: vec3_to_array(aabb.half_extent()),
        child_offset,
    }
}

fn sphere_to_meshlet(sphere: BoundingSphere) -> MeshletBoundingSphere {
    MeshletBoundingSphere {
        center: vec3_to_array(sphere.center()),
        radius: sphere.radius(),
    }
}

#[derive(Copy, Clone)]
struct Aabb {
    min: Vec3,
    max: Vec3,
}

impl Aabb {
    fn empty() -> Self {
        Self {
            min: Vec3::new(f32::INFINITY, f32::INFINITY, f32::INFINITY),
            max: Vec3::new(f32::NEG_INFINITY, f32::NEG_INFINITY, f32::NEG_INFINITY),
        }
    }

    fn from_points(points: impl Iterator<Item = Vec3>) -> Self {
        let mut aabb = Self::empty();
        for point in points {
            aabb.min = aabb.min.inf(&point);
            aabb.max = aabb.max.sup(&point);
        }
        aabb
    }

    fn from_center_half_extent(center: Vec3, half_extent: Vec3) -> Self {
        Self {
            min: center - half_extent,
            max: center + half_extent,
        }
    }

    fn center(&self) -> Vec3 {
        (self.min + self.max) * 0.5
    }

    fn half_extent(&self) -> Vec3 {
        (self.max - self.min) * 0.5
    }

    fn merge(&self, other: &Aabb) -> Aabb {
        Aabb {
            min: self.min.inf(&other.min),
            max: self.max.sup(&other.max),
        }
    }

    fn visible_area(&self) -> f32 {
        let extent = self.max - self.min;
        extent.x * extent.y + extent.y * extent.z + extent.x * extent.z
    }
}

#[derive(Copy, Clone)]
struct BoundingSphere {
    center: Vec3,
    radius: f32,
}

impl BoundingSphere {
    fn new(center: Vec3, radius: f32) -> Self {
        Self { center, radius }
    }

    fn center(&self) -> Vec3 {
        self.center
    }

    fn radius(&self) -> f32 {
        self.radius
    }

    fn center_distance(&self, other: &BoundingSphere) -> f32 {
        (self.center - other.center).norm()
    }
}

#[derive(Copy, Clone)]
struct TempMeshletCullData {
    aabb: Aabb,
    lod_group_sphere: BoundingSphere,
    error: f32,
}

#[derive(Clone)]
struct TempMeshletGroup {
    aabb: Aabb,
    lod_bounds: BoundingSphere,
    parent_error: f32,
    meshlets: Vec<u32>,
}

impl Default for TempMeshletGroup {
    fn default() -> Self {
        Self {
            aabb: Aabb::empty(),
            lod_bounds: BoundingSphere::new(Vec3::new(0.0, 0.0, 0.0), 0.0),
            parent_error: f32::MAX,
            meshlets: Vec::new(),
        }
    }
}

struct VertexDataAccumulator {
    vertex_positions: BitVec<u32, Lsb0>,
    vertex_normals: Vec<u32>,
    vertex_uvs: Vec<[f32; 2]>,
    meshlets: Vec<Meshlet>,
}

impl VertexDataAccumulator {
    fn append_meshlet(
        &mut self,
        meshlet: &meshopt_Meshlet,
        meshlet_vertex_ids: &[u32],
        vertex_buffer: &[u8],
        vertex_stride: usize,
        vertex_position_quantization_factor: u8,
    ) {
        let start_vertex_position_bit = self.vertex_positions.len() as u32;
        let start_vertex_attribute_id = self.vertex_normals.len() as u32;

        let quantization_factor =
            (1_i32 << vertex_position_quantization_factor) as f32 * CENTIMETERS_PER_METER;

        let mut min_channels = [i32::MAX; 3];
        let mut max_channels = [i32::MIN; 3];
        let mut quantized_positions = [[0_i32; 3]; 256];

        for (index, &vertex_id) in meshlet_vertex_ids.iter().enumerate() {
            let vertex_start = vertex_id as usize * vertex_stride;
            let vertex_data = &vertex_buffer[vertex_start..vertex_start + vertex_stride];
            let position = read_vec3(&vertex_data[0..12]);
            let normal = read_vec3(&vertex_data[12..24]);
            let texture_coordinates = read_vec2(&vertex_data[24..32]);

            self.vertex_uvs
                .push([texture_coordinates.x, texture_coordinates.y]);
            self.vertex_normals
                .push(pack2x16snorm(octahedral_encode(normal)));

            let quantized = [
                (position.x * quantization_factor + 0.5) as i32,
                (position.y * quantization_factor + 0.5) as i32,
                (position.z * quantization_factor + 0.5) as i32,
            ];
            quantized_positions[index] = quantized;
            for channel in 0..3 {
                min_channels[channel] = min_channels[channel].min(quantized[channel]);
                max_channels[channel] = max_channels[channel].max(quantized[channel]);
            }
        }

        let bits_per_channel = [
            ((max_channels[0] - min_channels[0] + 1) as f32)
                .log2()
                .ceil() as u8,
            ((max_channels[1] - min_channels[1] + 1) as f32)
                .log2()
                .ceil() as u8,
            ((max_channels[2] - min_channels[2] + 1) as f32)
                .log2()
                .ceil() as u8,
        ];

        for quantized in quantized_positions.iter().take(meshlet_vertex_ids.len()) {
            let remapped = [
                (quantized[0] - min_channels[0]) as u32,
                (quantized[1] - min_channels[1]) as u32,
                (quantized[2] - min_channels[2]) as u32,
            ];
            self.vertex_positions.extend_from_bitslice(
                &remapped[0].view_bits::<Lsb0>()[..bits_per_channel[0] as usize],
            );
            self.vertex_positions.extend_from_bitslice(
                &remapped[1].view_bits::<Lsb0>()[..bits_per_channel[1] as usize],
            );
            self.vertex_positions.extend_from_bitslice(
                &remapped[2].view_bits::<Lsb0>()[..bits_per_channel[2] as usize],
            );
        }

        self.meshlets.push(Meshlet {
            start_vertex_position_bit,
            start_vertex_attribute_id,
            start_index_id: meshlet.triangle_offset,
            vertex_count_minus_one: (meshlet.vertex_count - 1) as u8,
            triangle_count: meshlet.triangle_count as u8,
            padding: 0,
            bits_per_vertex_position_channel_x: bits_per_channel[0],
            bits_per_vertex_position_channel_y: bits_per_channel[1],
            bits_per_vertex_position_channel_z: bits_per_channel[2],
            vertex_position_quantization_factor,
            min_vertex_position_channel_x: min_channels[0] as f32,
            min_vertex_position_channel_y: min_channels[1] as f32,
            min_vertex_position_channel_z: min_channels[2] as f32,
        });
    }
}

struct TempBvhNode {
    group: u32,
    aabb: Aabb,
    children: Vec<u32>,
}

#[derive(Default)]
struct BvhBuilder {
    nodes: Vec<TempBvhNode>,
    lods: Vec<Range<u32>>,
}

impl BvhBuilder {
    fn add_lod(&mut self, offset: u32, all_groups: &[TempMeshletGroup]) {
        let first = self.nodes.len() as u32;
        self.nodes
            .extend(all_groups.iter().enumerate().skip(offset as usize).map(
                |(group_index, group)| TempBvhNode {
                    group: group_index as u32,
                    aabb: group.aabb,
                    children: Vec::new(),
                },
            ));
        let end = self.nodes.len() as u32;
        if first != end {
            self.lods.push(first..end);
        }
    }

    fn surface_area(&self, nodes: &[u32]) -> f32 {
        nodes
            .iter()
            .map(|&node| self.nodes[node as usize].aabb)
            .reduce(|accumulated, next| accumulated.merge(&next))
            .expect("cannot compute surface area of zero nodes")
            .visible_area()
    }

    fn node_center(&self, node: u32, axis: usize) -> f32 {
        self.nodes[node as usize].aabb.center()[axis]
    }

    fn sort_nodes_by_sah(&self, nodes: &mut [u32], splits: [usize; 8]) {
        for level in 0..3_usize {
            let parts = 1_usize << level;
            let nodes_per_split = 8_usize >> level;
            let half_count = nodes_per_split / 2;
            let mut offset = 0;
            for part in 0..parts {
                let first = part * nodes_per_split;
                let mut left_sum = 0;
                let mut right_sum = 0;
                for element in 0..half_count {
                    left_sum += splits[first + element];
                    right_sum += splits[first + half_count + element];
                }
                let total = left_sum + right_sum;
                let nodes = &mut nodes[offset..offset + total];
                offset += total;

                let mut cost = f32::MAX;
                let mut axis = 0;
                for candidate_axis in 0..3_usize {
                    nodes.sort_unstable_by(|&left, &right| {
                        self.node_center(left, candidate_axis)
                            .partial_cmp(&self.node_center(right, candidate_axis))
                            .unwrap()
                    });
                    let (left_nodes, right_nodes) = nodes.split_at(left_sum);
                    let candidate_cost =
                        self.surface_area(left_nodes) + self.surface_area(right_nodes);
                    if candidate_cost < cost {
                        axis = candidate_axis;
                        cost = candidate_cost;
                    }
                }
                if axis != 2 {
                    nodes.sort_unstable_by(|&left, &right| {
                        self.node_center(left, axis)
                            .partial_cmp(&self.node_center(right, axis))
                            .unwrap()
                    });
                }
            }
        }
    }

    fn build_temp_inner(&mut self, nodes: &mut [u32], optimize: bool) -> u32 {
        let count = nodes.len();
        if count == 1 {
            nodes[0]
        } else if count <= 8 {
            let node_index = self.nodes.len();
            self.nodes.push(TempBvhNode {
                group: u32::MAX,
                aabb: Aabb::empty(),
                children: nodes.to_vec(),
            });
            node_index as u32
        } else {
            let max_child_size = 1_usize << ((count.ilog2() / 3) * 3);
            let min_child_size = max_child_size >> 3;
            let max_extra_per_node = max_child_size - min_child_size;
            let mut extra = count - max_child_size;
            let splits: [usize; 8] = std::array::from_fn(|_| {
                let size = extra.min(max_extra_per_node);
                extra -= size;
                min_child_size + size
            });

            if optimize {
                self.sort_nodes_by_sah(nodes, splits);
            }

            let mut offset = 0;
            let children = splits
                .into_iter()
                .map(|size| {
                    let child = self.build_temp_inner(&mut nodes[offset..offset + size], optimize);
                    offset += size;
                    child
                })
                .collect();

            let node_index = self.nodes.len();
            self.nodes.push(TempBvhNode {
                group: u32::MAX,
                aabb: Aabb::empty(),
                children,
            });
            node_index as u32
        }
    }

    fn build_temp(&mut self) -> u32 {
        let mut lod_roots = Vec::with_capacity(self.lods.len());
        for lod in std::mem::take(&mut self.lods) {
            let mut lod: Vec<u32> = lod.collect();
            let root = self.build_temp_inner(&mut lod, true);
            let node = &self.nodes[root as usize];
            if node.group != u32::MAX || node.children.len() == 8 {
                lod_roots.push(root);
            } else {
                lod_roots.extend(node.children.iter().copied());
            }
        }
        self.build_temp_inner(&mut lod_roots, false)
    }

    fn build_inner(
        &self,
        groups: &[TempMeshletGroup],
        out: &mut Vec<BvhNode>,
        max_depth: &mut u32,
        node: u32,
        depth: u32,
    ) -> u32 {
        *max_depth = depth.max(*max_depth);
        let node_reference = &self.nodes[node as usize];
        let output_index = out.len();
        out.push(BvhNode::default());

        for (slot, &child_id) in node_reference.children.iter().enumerate() {
            let child = &self.nodes[child_id as usize];
            if child.group != u32::MAX {
                let group = &groups[child.group as usize];
                let output = &mut out[output_index];
                output.aabbs[slot] =
                    aabb_to_meshlet(group.aabb, group.parent_error, group.meshlets[0]);
                output.lod_bounds[slot] = sphere_to_meshlet(group.lod_bounds);
                output.child_counts[slot] = group.meshlets[1] as u8;
            } else {
                let child_output_index =
                    self.build_inner(groups, out, max_depth, child_id, depth + 1);
                let child_output = out[child_output_index as usize];
                let mut aabb = Aabb::empty();
                let mut parent_error = 0.0_f32;
                let mut lod_bounds = BoundingSphere::new(Vec3::new(0.0, 0.0, 0.0), 0.0);
                for child_slot in 0..8 {
                    if child_output.child_counts[child_slot] == 0 {
                        break;
                    }

                    aabb = aabb.merge(&Aabb::from_center_half_extent(
                        array_to_vec3(child_output.aabbs[child_slot].center),
                        array_to_vec3(child_output.aabbs[child_slot].half_extent),
                    ));
                    lod_bounds = merge_spheres(
                        lod_bounds,
                        BoundingSphere::new(
                            array_to_vec3(child_output.lod_bounds[child_slot].center),
                            child_output.lod_bounds[child_slot].radius,
                        ),
                    );
                    parent_error = parent_error.max(child_output.aabbs[child_slot].error);
                }

                let output = &mut out[output_index];
                output.aabbs[slot] = aabb_to_meshlet(aabb, parent_error, child_output_index);
                output.lod_bounds[slot] = sphere_to_meshlet(lod_bounds);
                output.child_counts[slot] = u8::MAX;
            }
        }

        output_index as u32
    }

    fn build(
        mut self,
        meshlets: &mut Meshlets,
        mut groups: Vec<TempMeshletGroup>,
        cull_data: &mut Vec<TempMeshletCullData>,
    ) -> (Vec<BvhNode>, MeshletAabb, u32) {
        let mut remap = Vec::with_capacity(meshlets.meshlets.len());
        let mut remapped_cull_data = Vec::with_capacity(cull_data.len());
        for group in groups.iter_mut() {
            let first = remap.len() as u32;
            let count = group.meshlets.len() as u32;
            remap.extend(
                group
                    .meshlets
                    .iter()
                    .map(|&meshlet_id| meshlets.meshlets[meshlet_id as usize]),
            );
            remapped_cull_data.extend(
                group
                    .meshlets
                    .iter()
                    .map(|&meshlet_id| cull_data[meshlet_id as usize]),
            );
            assert!(
                count < u8::MAX as u32,
                "a meshlet group holds {count} meshlets, which a slot's u8 count cannot \
                 address: u8::MAX is the marker for a slot that points at another node, so \
                 such a group would be walked as one"
            );
            group.meshlets.resize(2, 0);
            group.meshlets[0] = first;
            group.meshlets[1] = count;
        }
        meshlets.meshlets = remap;
        *cull_data = remapped_cull_data;

        let mut out = Vec::new();
        let mut aabb = Aabb::empty();
        let mut max_depth = 0;

        if self.nodes.len() == 1 {
            let mut node = BvhNode::default();
            let group = &groups[0];
            node.aabbs[0] = aabb_to_meshlet(group.aabb, group.parent_error, group.meshlets[0]);
            node.lod_bounds[0] = sphere_to_meshlet(group.lod_bounds);
            node.child_counts[0] = group.meshlets[1] as u8;
            out.push(node);
            aabb = group.aabb;
            max_depth = 1;
        } else {
            let root = self.build_temp();
            let root = self.build_inner(&groups, &mut out, &mut max_depth, root, 1);
            assert_eq!(root, 0, "bvh root must be node zero");

            let root_node = out[0];
            for slot in 0..8 {
                if root_node.child_counts[slot] == 0 {
                    break;
                }

                aabb = aabb.merge(&Aabb::from_center_half_extent(
                    array_to_vec3(root_node.aabbs[slot].center),
                    array_to_vec3(root_node.aabbs[slot].half_extent),
                ));
            }
        }

        let mut reachable = vec![false; meshlets.meshlets.len()];
        verify_bvh(&out, cull_data, &mut reachable, 0);
        assert!(
            reachable.iter().all(|&value| value),
            "all meshlets must be reachable"
        );

        (
            out,
            MeshletAabb {
                center: vec3_to_array(aabb.center()),
                half_extent: vec3_to_array(aabb.half_extent()),
                ..Default::default()
            },
            max_depth,
        )
    }
}

fn verify_bvh(
    out: &[BvhNode],
    cull_data: &[TempMeshletCullData],
    reachable: &mut [bool],
    node: u32,
) {
    let node = out[node as usize];
    for slot in 0..8 {
        let sphere = node.lod_bounds[slot];
        let error = node.aabbs[slot].error;
        if node.child_counts[slot] == u8::MAX {
            let child_offset = node.aabbs[slot].child_offset;
            let child = out[child_offset as usize];
            for child_slot in 0..8 {
                if child.child_counts[child_slot] == 0 {
                    break;
                }
                assert!(
                    child.aabbs[child_slot].error <= error,
                    "bvh errors are not monotonic"
                );
                let sphere_error = (array_to_vec3(sphere.center)
                    - array_to_vec3(child.lod_bounds[child_slot].center))
                .norm()
                    - (sphere.radius - child.lod_bounds[child_slot].radius);
                assert!(sphere_error <= 0.0001, "bvh lod spheres are not monotonic");
            }
            verify_bvh(out, cull_data, reachable, child_offset);
        } else {
            for meshlet_offset in 0..node.child_counts[slot] as u32 {
                let meshlet_id = (meshlet_offset + node.aabbs[slot].child_offset) as usize;
                let meshlet = &cull_data[meshlet_id];
                assert!(meshlet.error <= error, "meshlet errors are not monotonic");
                let sphere_error =
                    (array_to_vec3(sphere.center) - meshlet.lod_group_sphere.center()).norm()
                        - (sphere.radius - meshlet.lod_group_sphere.radius());
                assert!(
                    sphere_error <= 0.0001,
                    "meshlet lod spheres are not monotonic"
                );
                reachable[meshlet_id] = true;
            }
        }
    }
}

/// An error produced by [`from_mesh`].
#[derive(Debug, thiserror::Error)]
pub enum MeshToMeshletMeshConversionError {
    /// The mesh index count is not a multiple of three, so it is not a triangle list.
    #[error("mesh index count is not divisible by three")]
    WrongMeshPrimitiveTopology,
    /// The mesh has no indices, or none that enclose any area.
    #[error("mesh has no indices")]
    MeshMissingIndices,
    /// An index names a vertex the mesh does not have.
    #[error("mesh index names vertex {vertex} but the mesh has {vertex_count} vertices")]
    IndexOutOfBounds { vertex: u32, vertex_count: usize },
    /// A meshopt operation failed.
    #[error("meshopt failed to process the mesh: {0}")]
    Meshopt(#[from] meshopt::Error),
    /// Building the metis partition graph failed.
    #[error("metis failed to build the partition graph: {0}")]
    MetisGraph(#[from] metis::NewGraphError),
    /// Partitioning the mesh with metis failed.
    #[error("metis failed to partition the mesh: {0}")]
    MetisPartition(#[from] metis::Error),
}