polyhedron 0.1.0

A half edge and radial edge implementation.
Documentation
use std::ops::Mul;

use crate::{MeshStorage, PolyHedron, VertData};

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
enum TagColor
{
    BLACK,
    BLUE,
}

/// Subdivide the faces of a mesh topologically. New vertices are placed
/// on the edges.
pub fn linear_subdivision<S, C: MeshStorage, const MANIFOLD: bool>(
    mesh: &mut PolyHedron<C, MANIFOLD>,
) where
    VertData<C>: linear_isomorphic::InnerSpace<S>,
    S: linear_isomorphic::RealField + Mul<VertData<C>, Output = VertData<C>>,
{
    let old_edge_count = mesh.edge_count();
    let cutoff_vert_index = mesh.vert_count();
    // For each existing *edge* (not half edge) we will add 6 new half edges.
    let mut edge_flags: Vec<TagColor> =
        Vec::with_capacity(old_edge_count + (old_edge_count / 2) * 6);
    for hedge_handle in mesh.iter_edges()
    {
        let v1 = hedge_handle.source().id();
        let v2 = hedge_handle.dest().id();

        let mid = hedge_handle.split();

        for h in mesh.vert_handle(mid).iter_star_half_edges()
        {
            let edge_index = h.edge_index();
            if edge_index >= edge_flags.len()
            {
                edge_flags.resize(edge_index + 1, TagColor::BLACK);
            }

            edge_flags[edge_index] = if h.dest().id() == v2 || h.dest().id() == v1
            {
                TagColor::BLACK
            }
            else
            {
                TagColor::BLUE
            };
        }
    }

    for hedge_handle in mesh.iter_edges().skip(old_edge_count)
    {
        if hedge_handle.is_in_boundary()
        {
            continue;
        }

        let v1_is_old = hedge_handle.source().id.to_index() < cutoff_vert_index;
        let v2_is_old = hedge_handle.dest().id.to_index() < cutoff_vert_index;

        if v1_is_old != v2_is_old
            && edge_flags[hedge_handle.edge_index()] == TagColor::BLUE
        {
            hedge_handle.flip();
        };
    }
}

/// Loop subdivision.
pub fn loop_subdivision<S, C: MeshStorage, const MANIFOLD: bool>(
    mesh: &mut PolyHedron<C, MANIFOLD>,
) where
    VertData<C>: linear_isomorphic::InnerSpace<S>,
    S: linear_isomorphic::RealField + Mul<VertData<C>, Output = VertData<C>>,
{
    // TODO(medium): Optimize this function to use less memory, ideally $O(1)$
    let mut vertex_new = Vec::<_>::with_capacity(mesh.vert_count() + mesh.edge_count());

    for vert in mesh.iter_verts()
    {
        let v_data = vert.read_data().clone();
        let mut average = v_data.clone() * S::from(0.0).unwrap();
        let mut valence = 0;
        for neighbour in vert.iter_neighbours()
        {
            average = average + neighbour.read_data().clone();
            valence += 1;
        }

        let n = S::from(valence).unwrap();
        let u = if valence == 3
        {
            S::from(3.0 / 16.0).unwrap()
        }
        else
        {
            S::from(3.0).unwrap() / (S::from(8.0).unwrap() * n)
        };

        average = v_data * (S::from(1.0).unwrap() - u * n) + average * u;
        vertex_new.push(average);
    }

    for handle in mesh.iter_edges()
    {
        let v1 = handle.source().read_data().clone();
        let v2 = handle.orbit_next().face_prev().source().read_data().clone();
        let v3 = handle.dest().read_data().clone();
        let v4 = handle.face_prev().source().read_data().clone();

        let average = (v1 + v3) * S::from(3.0 / 8.0).unwrap()
            + (v2 + v4) * S::from(1.0 / 8.0).unwrap();

        vertex_new.push(average);
    }

    linear_subdivision(mesh);

    for (mut vert, data) in mesh.iter_verts().zip(vertex_new.into_iter())
    {
        vert.update_data(move |v| *v = data.clone());
    }
}

#[cfg(test)]
mod tests
{
    use wavefront_loader::ObjData;

    use super::*;
    use crate::verification::verify_correctness;
    use crate::BasicPolyHedron;

    #[test]
    fn test_subdivision()
    {
        let ObjData {
            vertices,
            vertex_face_indices,
            ..
        } = ObjData::from_disk_file("../../../Assets/cube.obj");

        let mut test = BasicPolyHedron::<_, false>::new(
            vertices.clone(),
            (),
            vertex_face_indices
                .iter()
                .map(|l| l.iter().map(|i| *i as usize)),
        );
        verify_correctness(&test).unwrap();

        linear_subdivision(&mut test);

        let (vs, fs) = test.face_list();
        ObjData::export(&(&vs, &fs), "tmp/linear_subdivision_test.obj");

        let mut test = BasicPolyHedron::<_, false>::new(
            vertices.clone(),
            (),
            vertex_face_indices
                .iter()
                .map(|l| l.iter().map(|i| *i as usize)),
        );
        verify_correctness(&test).unwrap();

        for _ in 0..3
        {
            loop_subdivision(&mut test);
        }

        let (vs, fs) = test.face_list();
        ObjData::export(&(&vs, &fs), "tmp/loop_subdivision_test.obj");
    }
}