del-msh-cpu 0.1.49

//! method for 3D Bounding Volume Hierarchy

use num_traits::{AsPrimitive, PrimInt};

pub struct TriMesh3WithBvhRef<'a, Index> {
    pub tri2vtx: &'a [[Index; 3]],
    pub vtx2xyz: &'a [[f32; 3]],
    pub bvhnode2tri_btree: &'a [[Index; 3]],
    pub bvhnode2aabb: &'a [[f32; 6]],
}

pub fn intersections_ray<Index>(
    hits: &mut Vec<(f32, usize)>,
    ray_org: &[f32; 3],
    ray_dir: &[f32; 3],
    trimesh3: &TriMesh3WithBvhRef<Index>,
    i_bvhnode: usize,
) where
    Index: PrimInt + AsPrimitive<usize>,
{
    if del_geo_core::aabb3::from_aabbs(trimesh3.bvhnode2aabb.as_flattened(), i_bvhnode)
        .intersections_against_ray(ray_org, ray_dir)
        .is_none()
    {
        return;
    }
    assert_eq!(
        trimesh3.bvhnode2tri_btree.len(),
        trimesh3.bvhnode2aabb.len()
    );
    if trimesh3.bvhnode2tri_btree[i_bvhnode][2] == Index::max_value() {
        // leaf node
        let i_tri: usize = trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_();
        let Some((t, _bc)) = crate::trimesh3::to_tri3(trimesh3.tri2vtx, trimesh3.vtx2xyz, i_tri)
            .intersection_against_ray(ray_org, ray_dir)
        else {
            return;
        };
        hits.push((t, i_tri));
        return;
    }
    intersections_ray(
        hits,
        ray_org,
        ray_dir,
        trimesh3,
        trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_(),
    );
    intersections_ray(
        hits,
        ray_org,
        ray_dir,
        trimesh3,
        trimesh3.bvhnode2tri_btree[i_bvhnode][2].as_(),
    );
}

pub fn intersections_line<Index>(
    hits: &mut Vec<(f32, usize)>,
    line_org: &[f32; 3],
    line_dir: &[f32; 3],
    trimesh3: &TriMesh3WithBvhRef<Index>,
    i_bvhnode: usize,
) where
    Index: PrimInt + AsPrimitive<usize>,
{
    if del_geo_core::aabb3::from_aabbs(trimesh3.bvhnode2aabb.as_flattened(), i_bvhnode)
        .intersections_against_line(line_org, line_dir)
        .is_none()
    {
        return;
    }
    assert_eq!(
        trimesh3.bvhnode2tri_btree.len(),
        trimesh3.bvhnode2aabb.len()
    );
    if trimesh3.bvhnode2tri_btree[i_bvhnode][2] == Index::max_value() {
        // leaf node
        let i_tri: usize = trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_();
        let Some((t, _bc)) = crate::trimesh3::to_tri3(trimesh3.tri2vtx, trimesh3.vtx2xyz, i_tri)
            .intersection_against_line(line_org, line_dir)
        else {
            return;
        };
        hits.push((t, i_tri));
        return;
    }
    crate::search_bvh3::intersections_line(
        hits,
        line_org,
        line_dir,
        trimesh3,
        trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_(),
    );
    crate::search_bvh3::intersections_line(
        hits,
        line_org,
        line_dir,
        trimesh3,
        trimesh3.bvhnode2tri_btree[i_bvhnode][2].as_(),
    );
}

///
/// # return
///  - the coefficient of the ray to triangle
///  - triangle index
///  - barycentric coordinate
pub fn first_intersection_ray<Index>(
    ray_org: &[f32; 3],
    ray_dir: &[f32; 3],
    trimesh3: &TriMesh3WithBvhRef<Index>,
    i_bvhnode: usize,
    dis: f32,
) -> Option<(f32, usize, [f32; 3])>
where
    Index: PrimInt + AsPrimitive<usize>,
{
    assert_eq!(
        trimesh3.bvhnode2tri_btree.len(),
        trimesh3.tri2vtx.len() * 2 - 1
    );
    assert_eq!(
        trimesh3.bvhnode2tri_btree.len(),
        trimesh3.bvhnode2aabb.len()
    );
    // culling the branch
    del_geo_core::aabb3::from_aabbs(trimesh3.bvhnode2aabb.as_flattened(), i_bvhnode)
        .intersections_against_ray(ray_org, ray_dir)?;
    if trimesh3.bvhnode2tri_btree[i_bvhnode][2] == Index::max_value() {
        // leaf node
        let i_tri: usize = trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_();
        return crate::trimesh3::to_tri3(trimesh3.tri2vtx, trimesh3.vtx2xyz, i_tri)
            .intersection_against_ray(ray_org, ray_dir)
            .filter(|(t, _bc)| *t < dis)
            .map(|(t, bc)| (t, i_tri, bc));
    }

    // near branch is checked first. Check which branch (left or right) is near
    let near_is_left = {
        let t_aabb_left = del_geo_core::aabb3::from_aabbs(
            trimesh3.bvhnode2aabb.as_flattened(),
            trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_(),
        )
        .intersections_against_ray(ray_org, ray_dir);
        let t_aabb_right = del_geo_core::aabb3::from_aabbs(
            trimesh3.bvhnode2aabb.as_flattened(),
            trimesh3.bvhnode2tri_btree[i_bvhnode][2].as_(),
        )
        .intersections_against_ray(ray_org, ray_dir);
        match (t_aabb_left, t_aabb_right) {
            (Some(tl), Some(tr)) => tl < tr,
            (None, None) => {
                return None;
            }
            (Some(_), None) => true,
            (None, Some(_)) => false,
        }
    };

    // check near branch
    let res_near = {
        let i_child_near: usize = if near_is_left {
            trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_()
        } else {
            trimesh3.bvhnode2tri_btree[i_bvhnode][2].as_()
        };
        first_intersection_ray(ray_org, ray_dir, trimesh3, i_child_near, dis)
    };

    if let Some((t_near, _i_tri_near, _bc_near)) = res_near {
        // if there is a hit branch in the left
        let i_bvhnode_far: usize = if near_is_left {
            trimesh3.bvhnode2tri_btree[i_bvhnode][2].as_()
        } else {
            trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_()
        };
        let aabb_far = &trimesh3.bvhnode2aabb[i_bvhnode_far];
        // the intersection point is closer than another bvhnode
        if let Some((t_aabb_far, _)) =
            del_geo_core::aabb::intersections_against_ray(aabb_far, ray_org, ray_dir)
        {
            if t_aabb_far > t_near {
                // the ray hit the far branch, but it is farther than the current nearest
                return res_near;
            }
        } else {
            // the ray does not hit the far-branch.
            return res_near;
        }
    }

    let res_far = {
        let i_child_far: usize = if near_is_left {
            trimesh3.bvhnode2tri_btree[i_bvhnode][2].as_()
        } else {
            trimesh3.bvhnode2tri_btree[i_bvhnode][1].as_()
        };
        first_intersection_ray(ray_org, ray_dir, trimesh3, i_child_far, dis)
    };

    if let Some((t_near, i_tri_near, bc_near)) = res_near {
        if let Some((t_far, i_tri_far, bc_far)) = res_far {
            if t_near < t_far {
                Some((t_near, i_tri_near, bc_near))
            } else {
                Some((t_far, i_tri_far, bc_far))
            }
        } else {
            Some((t_near, i_tri_near, bc_near))
        }
    } else if let Some((t_far, i_tri_far, bc_far)) = res_far {
        Some((t_far, i_tri_far, bc_far))
    } else {
        None
    }
}

/*
/// check if a point alone ray_dir is closer than an aabb
fn is_point_closer(aabb: &[f32; 6], ray_dir: &[f32; 3], t: f32) -> bool {
    // closest projection of aabb along ray
    let mut proj_aabb = 0.;
    for i_dim in 0..3 {
        if ray_dir[i_dim].abs() == 0. {
            continue;
        }

        if ray_dir[i_dim] > 0. {
            // min project
            proj_aabb += aabb[i_dim] * ray_dir[i_dim];
        } else {
            // max project
            proj_aabb += aabb[i_dim + 3] * ray_dir[i_dim];
        }
    }

    let sqlen = del_geo_core::vec3::dot(ray_dir, ray_dir);
    t * sqlen < proj_aabb
}
 */

#[test]
fn test_first_intersection_ray() {
    let trimesh3 = crate::trimesh3_primitive::sphere_yup::<usize, f32>(1.0, 128, 128);
    let bvhnodes = crate::bvhnodes_morton::from_triangle_mesh(&trimesh3.tri2vtx, &trimesh3.vtx2xyz);
    let bvhnode2aabb = crate::bvhnode2aabb3::from_uniform_mesh_with_bvh(
        0,
        &bvhnodes,
        &trimesh3.tri2vtx,
        &trimesh3.vtx2xyz,
        None,
    );
    use rand::RngExt;
    use rand::SeedableRng;
    let mut reng = rand_chacha::ChaChaRng::seed_from_u64(0u64);
    for _iter in 0..1000 {
        let x = reng.random::<f32>() * 2.0 - 1.;
        let y = reng.random::<f32>() * 2.0 - 1.;
        let ray_dir = [0., 0., -1.];
        let ray_org = [x, y, 1.];
        let res_first = first_intersection_ray(
            &ray_org,
            &ray_dir,
            &TriMesh3WithBvhRef {
                tri2vtx: &trimesh3.tri2vtx,
                vtx2xyz: &trimesh3.vtx2xyz,
                bvhnode2tri_btree: &bvhnodes,
                bvhnode2aabb: &bvhnode2aabb,
            },
            0,
            f32::INFINITY,
        );
        if let Some((t, _i_tri, _bc)) = res_first {
            let p = [
                ray_org[0] + t * ray_dir[0],
                ray_org[1] + t * ray_dir[1],
                ray_org[2] + t * ray_dir[2],
            ];
            assert!(del_geo_core::vec3::norm(&p) < 1.0);
        }
        {
            // intersection all
            let mut hits = Vec::<(f32, usize)>::new();
            intersections_ray::<usize>(
                &mut hits,
                &ray_org,
                &ray_dir,
                &TriMesh3WithBvhRef {
                    tri2vtx: &trimesh3.tri2vtx,
                    vtx2xyz: &trimesh3.vtx2xyz,
                    bvhnode2tri_btree: &bvhnodes,
                    bvhnode2aabb: &bvhnode2aabb,
                },
                0,
            );
            if let Some((t_first, _i_tri_first, _bc_first)) = res_first {
                assert_eq!(hits.len(), 2);
                assert!(hits[0].0 == t_first || hits[1].0 == t_first);
            } else {
                assert!(hits.is_empty());
            }
        }
    }
}