vox_geometry_rust/

triangle_mesh3.rs

/*
 * // Copyright (c) 2021 Feng Yang
 * //
 * // I am making my contributions/submissions to this project solely in my
 * // personal capacity and am not conveying any rights to any intellectual
 * // property of any third parties.
 */

use crate::vector3::Vector3D;
use crate::vector2::Vector2D;
use crate::usize3::USize3;
use crate::transform3::Transform3;
use crate::surface3::*;
use crate::bvh3::Bvh3;
use crate::triangle3::Triangle3;
use crate::quaternion::QuaternionD;
use crate::ray3::Ray3D;
use crate::bounding_box3::BoundingBox3D;
use crate::nearest_neighbor_query_engine3::NearestNeighborQueryEngine3;
use crate::intersection_query_engine3::IntersectionQueryEngine3;
use crate::matrix3x3::Matrix3x3D;
use rayon::iter::ParallelIterator;
use rayon::prelude::IntoParallelRefMutIterator;
use std::mem::swap;
use std::sync::{RwLock, Arc};
use std::ops::Add;
use std::io::{Cursor, BufWriter, Write, Result};
use std::borrow::BorrowMut;
use std::fs::OpenOptions;

const K_DEFAULT_FAST_WINDING_NUMBER_ACCURACY: f64 = 2.0;

#[derive(Clone)]
struct WindingNumberGatherData {
    area_sums: f64,
    area_weighted_normal_sums: Vector3D,
    area_weighted_position_sums: Vector3D,
}

impl WindingNumberGatherData {
    pub fn new() -> WindingNumberGatherData {
        return WindingNumberGatherData {
            area_sums: 0.0,
            area_weighted_normal_sums: Vector3D::new_default(),
            area_weighted_position_sums: Vector3D::new_default(),
        };
    }
}

impl Add for WindingNumberGatherData {
    type Output = WindingNumberGatherData;

    fn add(self, rhs: Self) -> Self::Output {
        let mut sum = WindingNumberGatherData::new();
        sum.area_sums = self.area_sums + rhs.area_sums;
        sum.area_weighted_normal_sums =
            self.area_weighted_normal_sums + rhs.area_weighted_normal_sums;
        sum.area_weighted_position_sums =
            self.area_weighted_position_sums + rhs.area_weighted_position_sums;

        return sum;
    }
}

fn post_order_traversal<GatherFunc, LeafGatherFunc>(bvh: &Bvh3<usize>, node_index: usize, visitor_func: &mut GatherFunc,
                                                    leaf_func: &mut LeafGatherFunc,
                                                    init_gather_data: &WindingNumberGatherData) -> WindingNumberGatherData
    where GatherFunc: FnMut(usize, &WindingNumberGatherData),
          LeafGatherFunc: FnMut(usize) -> WindingNumberGatherData {
    let mut data = init_gather_data.clone();

    if bvh.is_leaf(node_index) {
        data = leaf_func(node_index);
    } else {
        let children = bvh.children(node_index);
        data = data + post_order_traversal(bvh, children.0, visitor_func,
                                           leaf_func, init_gather_data);
        data = data + post_order_traversal(bvh, children.1, visitor_func,
                                           leaf_func, init_gather_data);
    }
    visitor_func(node_index, &data);

    return data;
}

//--------------------------------------------------------------------------------------------------
///
/// # 3-D triangle mesh geometry.
///
/// This class represents 3-D triangle mesh geometry which extends Surface3 by
/// overriding surface-related queries. The mesh structure stores point,
/// normals, and UV coordinates.
///
pub struct TriangleMesh3 {
    _points: Vec<Vector3D>,
    _normals: Vec<Vector3D>,
    _uvs: Vec<Vector2D>,
    _point_indices: Vec<USize3>,
    _normal_indices: Vec<USize3>,
    _uv_indices: Vec<USize3>,

    _bvh: RwLock<Bvh3<usize>>,
    _bvh_invalidated: RwLock<bool>,

    _wn_area_weighted_normal_sums: RwLock<Vec<Vector3D>>,
    _wn_area_weighted_avg_positions: RwLock<Vec<Vector3D>>,
    _wn_invalidated: RwLock<bool>,

    /// data from surface3
    pub surface_data: Surface3Data,
}

impl TriangleMesh3 {
    /// Default constructor.
    /// ```
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// let mesh1 = TriangleMesh3::new_default(None, None);
    /// assert_eq!(0, mesh1.number_of_points());
    /// assert_eq!(0, mesh1.number_of_normals());
    /// assert_eq!(0, mesh1.number_of_uvs());
    /// assert_eq!(0, mesh1.number_of_triangles());
    /// ```
    pub fn new_default(transform: Option<Transform3>,
                       is_normal_flipped: Option<bool>) -> TriangleMesh3 {
        return TriangleMesh3 {
            _points: vec![],
            _normals: vec![],
            _uvs: vec![],
            _point_indices: vec![],
            _normal_indices: vec![],
            _uv_indices: vec![],
            _bvh: RwLock::new(Bvh3::new()),
            _bvh_invalidated: RwLock::new(true),
            _wn_area_weighted_normal_sums: RwLock::new(vec![]),
            _wn_area_weighted_avg_positions: RwLock::new(vec![]),
            _wn_invalidated: RwLock::new(true),
            surface_data: Surface3Data::new(transform, is_normal_flipped),
        };
    }

    /// Constructs mesh with points, normals, uvs, and their indices.
    pub fn new(points: Vec<Vector3D>, normals: Vec<Vector3D>,
               uvs: Vec<Vector2D>, point_indices: Vec<USize3>,
               normal_indices: Vec<USize3>, uv_indices: Vec<USize3>,
               transform: Option<Transform3>,
               is_normal_flipped: Option<bool>) -> TriangleMesh3 {
        return TriangleMesh3 {
            _points: points,
            _normals: normals,
            _uvs: uvs,
            _point_indices: point_indices,
            _normal_indices: normal_indices,
            _uv_indices: uv_indices,
            _bvh: RwLock::new(Bvh3::new()),
            _bvh_invalidated: RwLock::new(true),
            _wn_area_weighted_normal_sums: RwLock::new(vec![]),
            _wn_area_weighted_avg_positions: RwLock::new(vec![]),
            _wn_invalidated: RwLock::new(true),
            surface_data: Surface3Data::new(transform, is_normal_flipped),
        };
    }

    /// Returns builder fox TriangleMesh3.
    pub fn builder() -> Builder {
        return Builder::new();
    }

    /// Clears all content.
    pub fn clear(&mut self) {
        self._points.clear();
        self._normals.clear();
        self._uvs.clear();
        self._point_indices.clear();
        self._normal_indices.clear();
        self._uv_indices.clear();

        self.invalidate_cache();
    }

    /// Copies the contents from \p other mesh.
    pub fn set(&mut self, other: &TriangleMesh3) {
        self._points = other._points.clone();
        self._normals = other._normals.clone();
        self._uvs = other._uvs.clone();
        self._point_indices = other._point_indices.clone();
        self._normal_indices = other._normal_indices.clone();
        self._uv_indices = other._uv_indices.clone();

        self.invalidate_cache();
    }

    /// Swaps the contents with \p other mesh.
    pub fn swap(&mut self, other: &mut TriangleMesh3) {
        swap(&mut self._points, &mut other._points);
        swap(&mut self._normals, &mut other._normals);
        swap(&mut self._uvs, &mut other._uvs);
        swap(&mut self._point_indices, &mut other._point_indices);
        swap(&mut self._normal_indices, &mut other._normal_indices);
        swap(&mut self._uv_indices, &mut other._uv_indices);
    }

    /// Returns area of this mesh.
    pub fn area(&self) -> f64 {
        let mut a: f64 = 0.0;
        for i in 0..self.number_of_triangles() {
            let tri = self.triangle(i);
            a += tri.area();
        }
        return a;
    }

    /// Returns volume of this mesh.
    pub fn volume(&self) -> f64 {
        let mut vol: f64 = 0.0;
        for i in 0..self.number_of_triangles() {
            let tri = self.triangle(i);
            vol += tri.points[0].dot(&tri.points[1].cross(&tri.points[2])) / 6.0;
        }
        return vol;
    }

    /// Returns constant reference to the i-th point.
    pub fn point(&self, i: usize) -> Vector3D {
        return self._points[i];
    }

    /// Returns reference to the i-th point.
    pub fn point_mut(&mut self, i: usize) -> Vector3D {
        self.invalidate_cache();
        return self._points[i];
    }

    /// Returns constant reference to the i-th normal.
    pub fn normal(&self, i: usize) -> Vector3D {
        return self._normals[i];
    }

    /// Returns reference to the i-th normal.
    pub fn normal_mut(&mut self, i: usize) -> Vector3D {
        return self._normals[i];
    }

    /// Returns constant reference to the i-th UV coordinates.
    pub fn uv(&self, i: usize) -> Vector2D {
        return self._uvs[i];
    }

    /// Returns reference to the i-th UV coordinates.
    pub fn uv_mut(&mut self, i: usize) -> Vector2D {
        return self._uvs[i];
    }

    /// Returns constant reference to the point indices of i-th triangle.
    pub fn point_index(&self, i: usize) -> USize3 {
        return self._point_indices[i];
    }

    /// Returns reference to the point indices of i-th triangle.
    pub fn point_index_mut(&mut self, i: usize) -> USize3 {
        return self._point_indices[i];
    }

    /// Returns constant reference to the normal indices of i-th triangle.
    pub fn normal_index(&self, i: usize) -> USize3 {
        return self._normal_indices[i];
    }

    /// Returns reference to the normal indices of i-th triangle.
    pub fn normal_index_mut(&mut self, i: usize) -> USize3 {
        return self._normal_indices[i];
    }

    /// Returns constant reference to the UV indices of i-th triangle.
    pub fn uv_index(&self, i: usize) -> USize3 {
        return self._uv_indices[i];
    }

    /// Returns reference to the UV indices of i-th triangle.
    pub fn uv_index_mut(&mut self, i: usize) -> USize3 {
        return self._uv_indices[i];
    }

    /// Returns i-th triangle.
    pub fn triangle(&self, i: usize) -> Triangle3 {
        let mut tri = Triangle3::new_default(None, None);
        for j in 0..3 {
            tri.points[j] = self._points[self._point_indices[i][j]];
            if self.has_uvs() {
                tri.uvs[j] = self._uvs[self._uv_indices[i][j]];
            }
        }

        let n = tri.face_normal();

        for j in 0..3 {
            if self.has_normals() {
                tri.normals[j] = self._normals[self._normal_indices[i][j]];
            } else {
                tri.normals[j] = n;
            }
        }

        return tri;
    }

    /// Returns number of points.
    pub fn number_of_points(&self) -> usize {
        return self._points.len();
    }

    /// Returns number of normals.
    pub fn number_of_normals(&self) -> usize {
        return self._normals.len();
    }

    /// Returns number of UV coordinates.
    pub fn number_of_uvs(&self) -> usize {
        return self._uvs.len();
    }

    /// Returns number of triangles.
    pub fn number_of_triangles(&self) -> usize {
        return self._point_indices.len();
    }

    /// Returns true if the mesh has normals.
    pub fn has_normals(&self) -> bool {
        return self._normals.len() > 0;
    }

    /// Returns true if the mesh has UV coordinates.
    pub fn has_uvs(&self) -> bool {
        return self._uvs.len() > 0;
    }

    /// Adds a point.
    pub fn add_point(&mut self, pt: Vector3D) {
        self._points.push(pt);
    }

    /// Adds a normal.
    pub fn add_normal(&mut self, n: Vector3D) {
        self._normals.push(n);
    }

    /// Adds a UV.
    pub fn add_uv(&mut self, t: Vector2D) {
        self._uvs.push(t);
    }

    /// Adds a triangle with points.
    pub fn add_point_triangle(&mut self, new_point_indices: USize3) {
        self._point_indices.push(new_point_indices);
        self.invalidate_cache();
    }

    /// Adds a triangle with normal.
    pub fn add_normal_triangle(&mut self, new_normal_indices: USize3) {
        self._normal_indices.push(new_normal_indices);

        self.invalidate_cache();
    }

    /// Adds a triangle with UV.
    pub fn add_uv_triangle(&mut self, new_uv_indices: USize3) {
        self._uv_indices.push(new_uv_indices);

        self.invalidate_cache();
    }

    /// Adds a triangle with point and normal.
    pub fn add_point_normal_triangle(&mut self, new_point_indices: USize3,
                                     new_normal_indices: USize3) {
        self._point_indices.push(new_point_indices);
        self._normal_indices.push(new_normal_indices);

        self.invalidate_cache();
    }

    /// Adds a triangle with point, normal, and UV.
    pub fn add_point_uv_normal_triangle(&mut self, new_point_indices: USize3,
                                        new_uv_indices: USize3,
                                        new_normal_indices: USize3) {
        self._point_indices.push(new_point_indices);
        self._normal_indices.push(new_normal_indices);
        self._uv_indices.push(new_uv_indices);

        self.invalidate_cache();
    }

    /// Adds a triangle with point and UV.
    pub fn add_point_uv_triangle(&mut self, new_point_indices: USize3,
                                 new_uv_indices: USize3) {
        self._point_indices.push(new_point_indices);
        self._uv_indices.push(new_uv_indices);

        self.invalidate_cache();
    }

    /// Add a triangle.
    pub fn add_triangle(&mut self, tri: Triangle3) {
        let v_start = self._points.len();
        let n_start = self._normals.len();
        let t_start = self._uvs.len();
        let mut new_point_indices = USize3::new_default();
        let mut new_normal_indices = USize3::new_default();
        let mut new_uv_indices = USize3::new_default();
        for i in 0..3 {
            self._points.push(tri.points[i]);
            self._normals.push(tri.normals[i]);
            self._uvs.push(tri.uvs[i]);
            new_point_indices[i] = v_start + i;
            new_normal_indices[i] = n_start + i;
            new_uv_indices[i] = t_start + i;
        }
        self._point_indices.push(new_point_indices);
        self._normal_indices.push(new_normal_indices);
        self._uv_indices.push(new_uv_indices);

        self.invalidate_cache();
    }

    /// Sets entire normals to the face normals.
    pub fn set_face_normal(&mut self) {
        self._normals.resize(self._points.len(), Vector3D::new_default());
        self._normal_indices = self._point_indices.clone();

        for i in 0..self.number_of_triangles() {
            let tri = self.triangle(i);
            let n = tri.face_normal();
            let f = self._point_indices[i];
            self._normals[f.x] = n;
            self._normals[f.y] = n;
            self._normals[f.z] = n;
        }
    }

    /// Sets angle weighted vertex normal.
    pub fn set_angle_weighted_vertex_normal(&mut self) {
        self._normals.clear();
        self._normal_indices.clear();

        let mut angle_weights: Vec<f64> = Vec::new();
        angle_weights.resize(self._points.len(), 0.0);
        let mut pseudo_normals: Vec<Vector3D> = Vec::new();
        pseudo_normals.resize(self._points.len(), Vector3D::new_default());

        for i in 0..self._points.len() {
            angle_weights[i] = 0.0;
            pseudo_normals[i] = Vector3D::new_default();
        }

        for i in 0..self.number_of_triangles() {
            let mut pts: [Vector3D; 3] = [Vector3D::new_default(); 3];
            let mut normal: Vector3D;
            let mut e0: Vector3D;
            let mut e1: Vector3D;
            let mut cosangle: f64;
            let mut angle: f64;
            let mut idx: [usize; 3] = [0; 3];

            // Quick references
            for j in 0..3 {
                idx[j] = self._point_indices[i][j];
                pts[j] = self._points[idx[j]];
            }

            // Angle for point 0
            e0 = pts[1] - pts[0];
            e1 = pts[2] - pts[0];
            e0.normalize();
            e1.normalize();
            normal = e0.cross(&e1);
            normal.normalize();
            cosangle = crate::math_utils::clamp(e0.dot(&e1), -1.0, 1.0);
            angle = f64::acos(cosangle);
            angle_weights[idx[0]] += angle;
            pseudo_normals[idx[0]] += normal * angle;

            // Angle for point 1
            e0 = pts[2] - pts[1];
            e1 = pts[0] - pts[1];
            e0.normalize();
            e1.normalize();
            normal = e0.cross(&e1);
            normal.normalize();
            cosangle = crate::math_utils::clamp(e0.dot(&e1), -1.0, 1.0);
            angle = f64::acos(cosangle);
            angle_weights[idx[1]] += angle;
            pseudo_normals[idx[1]] += normal * angle;

            // Angle for point 2
            e0 = pts[0] - pts[2];
            e1 = pts[1] - pts[2];
            e0.normalize();
            e1.normalize();
            normal = e0.cross(&e1);
            normal.normalize();
            cosangle = crate::math_utils::clamp(e0.dot(&e1), -1.0, 1.0);
            angle = f64::acos(cosangle);
            angle_weights[idx[2]] += angle;
            pseudo_normals[idx[2]] += normal * angle;
        }

        for i in 0..self._points.len() {
            if angle_weights[i] > 0.0 {
                pseudo_normals[i] /= angle_weights[i];
            }
        }

        swap(&mut pseudo_normals, &mut self._normals);
        self._normal_indices = self._point_indices.clone();
    }

    /// Scales the mesh by given factor.
    pub fn scale(&mut self, factor: f64) {
        self._points.par_iter_mut().for_each(|x| {
            *x *= factor;
        });

        self.invalidate_cache();
    }

    /// Translates the mesh.
    pub fn translate(&mut self, t: Vector3D) {
        self._points.par_iter_mut().for_each(|x| {
            *x += t;
        });

        self.invalidate_cache();
    }

    /// Rotates the mesh.
    pub fn rotate(&mut self, q: QuaternionD) {
        self._points.par_iter_mut().for_each(|x| {
            *x = q * (*x);
        });

        self._normals.par_iter_mut().for_each(|x| {
            *x = q * (*x);
        });

        self.invalidate_cache();
    }
}

impl TriangleMesh3 {
    fn invalidate_cache(&self) {
        *(self._bvh_invalidated.write().unwrap()) = true;
        *(self._wn_invalidated.write().unwrap()) = true;
    }

    fn build_bvh(&self) {
        if *self._bvh_invalidated.read().unwrap() {
            let n_tris = self.number_of_triangles();
            let mut ids: Vec<usize> = Vec::new();
            ids.resize(n_tris, 0);
            let mut bounds: Vec<BoundingBox3D> = Vec::new();
            bounds.resize(n_tris, BoundingBox3D::new_default());
            for i in 0..n_tris {
                ids[i] = i;
                bounds[i] = self.triangle(i).bounding_box();
            }
            self._bvh.write().unwrap().build(&ids, &bounds);
            *(self._bvh_invalidated.write().unwrap()) = false;
        }
    }

    fn build_winding_numbers(&self) {
        // Barill et al., Fast Winding Numbers for Soups and Clouds, ACM SIGGRAPH
        // 2018
        if *self._wn_invalidated.read().unwrap() {
            self.build_bvh();

            let n_nodes = self._bvh.read().unwrap().number_of_nodes();
            self._wn_area_weighted_normal_sums.write().unwrap().resize(n_nodes, Vector3D::new_default());
            self._wn_area_weighted_avg_positions.write().unwrap().resize(n_nodes, Vector3D::new_default());

            let mut visitor_func = |node_index: usize, data: &WindingNumberGatherData| {
                self._wn_area_weighted_normal_sums.write().unwrap()[node_index] = data.area_weighted_normal_sums;
                self._wn_area_weighted_avg_positions.write().unwrap()[node_index] = data.area_weighted_position_sums / data.area_sums;
            };
            let mut leaf_func = |node_index: usize| -> WindingNumberGatherData {
                let mut result = WindingNumberGatherData::new();

                let iter = self._bvh.read().unwrap().item_of_node(node_index);

                let tri = self.triangle(iter);
                let area = tri.area();
                result.area_sums = area;
                result.area_weighted_normal_sums = tri.face_normal() * area;
                result.area_weighted_position_sums =
                    (tri.points[0] + tri.points[1] + tri.points[2]) * area / 3.0;

                return result;
            };

            post_order_traversal(&*self._bvh.read().unwrap(), 0, &mut visitor_func, &mut leaf_func,
                                 &WindingNumberGatherData::new());

            *(self._wn_invalidated.write().unwrap()) = false;
        }
    }

    fn winding_number(&self, query_point: &Vector3D, tri_index: usize) -> f64 {
        // Jacobson et al., Robust Inside-Outside Segmentation using Generalized
        // Winding Numbers, ACM SIGGRAPH 2013.
        let vi = self._points[self._point_indices[tri_index][0]];
        let vj = self._points[self._point_indices[tri_index][1]];
        let vk = self._points[self._point_indices[tri_index][2]];
        let va = vi - *query_point;
        let vb = vj - *query_point;
        let vc = vk - *query_point;
        let a = va.length();
        let b = vb.length();
        let c = vc.length();

        let mat = Matrix3x3D::new(va.x, vb.x, vc.x,
                                  va.y, vb.y, vc.y,
                                  va.z, vb.z, vc.z);
        let det = mat.determinant();
        let denom =
            a * b * c + va.dot(&vb) * c + vb.dot(&vc) * a + vc.dot(&va) * b;

        let solid_angle = 2.0 * f64::atan2(det, denom);

        return solid_angle;
    }

    fn fast_winding_number(&self, query_point: &Vector3D, accuracy: f64) -> f64 {
        self.build_winding_numbers();

        return self.fast_winding_number_internal(query_point, 0, accuracy);
    }

    fn fast_winding_number_internal(&self, q: &Vector3D, root_node_index: usize,
                                    accuracy: f64) -> f64 {
        // Barill et al., Fast Winding Numbers for Soups and Clouds, ACM SIGGRAPH
        // 2018.
        let tree_p = self._wn_area_weighted_avg_positions.read().unwrap()[root_node_index];
        let q_to_p2 = q.distance_squared_to(tree_p);

        let tree_n = self._wn_area_weighted_normal_sums.read().unwrap()[root_node_index];
        let tree_bound = self._bvh.read().unwrap().node_bound(root_node_index);
        let tree_rvec =
            crate::vector3::max(&(tree_p - tree_bound.lower_corner), &(tree_bound.lower_corner - tree_p));
        let tree_r = tree_rvec.length();

        return if q_to_p2 > crate::math_utils::square(accuracy * tree_r) {
            // Case: q is sufficiently far from all elements in tree
            // TODO: This is zero-th order approximation. Higher-order approximation
            // from Section 3.2.1 could be implemented for better accuracy in the
            // future.
            (tree_p - *q).dot(&tree_n) / (crate::constants::K_FOUR_PI_D * crate::math_utils::cubic(f64::sqrt(q_to_p2)))
        } else {
            if self._bvh.read().unwrap().is_leaf(root_node_index) {
                // Case: q is nearby; use direct sum for tree’s elements
                let iter = self._bvh.read().unwrap().item_of_node(root_node_index);
                self.winding_number(q, iter) * crate::constants::K_FOUR_PI_D
            } else {
                // Case: Recursive call
                let children = self._bvh.read().unwrap().children(root_node_index);
                let mut wn = 0.0;
                wn += self.fast_winding_number_internal(q, children.0, accuracy);
                wn += self.fast_winding_number_internal(q, children.1, accuracy);
                wn
            }
        };
    }
}

impl TriangleMesh3 {
    /// Writes the mesh in obj format to the file.
    pub fn write_obj(&self, filename: String) -> Result<()> {
        let f = OpenOptions::new().write(true).create(true)
            .open(filename).unwrap();
        let mut f = BufWriter::new(f);
        // vertex
        for pt in &self._points {
            write!(f, "v {}\n", pt)?;
        }

        // uv coords
        for uv in &self._uvs {
            write!(f, "vt {}\n", uv)?;
        }

        // normals
        for n in &self._normals {
            write!(f, "vn {}\n", n)?;
        }

        // faces
        let has_uvs = self.has_uvs();
        let has_normals = self.has_normals();
        for i in 0..self.number_of_triangles() {
            write!(f, "f ")?;
            for j in 0..3 {
                write!(f, "{}", self._point_indices[i][j] + 1)?;
                if has_normals || has_uvs {
                    write!(f, "/")?;
                }
                if has_uvs {
                    write!(f, "{}", self._uv_indices[i][j] + 1)?;
                }
                if has_normals {
                    write!(f, "/{}", self._normal_indices[i][j] + 1)?;
                }
                write!(f, " ")?;
            }
            write!(f, "\n")?;
        }

        f.flush()?;

        return Ok(());
    }


    /// Reads the mesh in obj format from the buffer.
    /// ```
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::unit_tests_utils::get_cube_tri_mesh3x3x3obj;
    /// use std::io::Cursor;
    /// let mut buffer = Cursor::new(get_cube_tri_mesh3x3x3obj());
    ///
    /// let mut mesh = TriangleMesh3::new_default(None, None);
    /// mesh.read_obj_buffer(&mut buffer);
    ///
    /// assert_eq!(56, mesh.number_of_points());
    /// assert_eq!(96, mesh.number_of_normals());
    /// assert_eq!(76, mesh.number_of_uvs());
    /// assert_eq!(108, mesh.number_of_triangles());
    /// ```
    pub fn read_obj_buffer(&mut self, strm: &mut Cursor<&[u8]>) -> bool {
        let (models, _) = tobj::load_obj_buf(
            &mut strm.borrow_mut(),
            &tobj::LoadOptions::default(),
            |_| { unreachable!(); })
            .expect("Failed to OBJ load buffer");

        return self.read_model(models);
    }

    /// Reads the mesh in obj format from the file.
    pub fn read_obj_file(&mut self, filename: String) -> bool {
        let (models, _) =
            tobj::load_obj(
                &filename,
                &tobj::LoadOptions::default(),
            ).expect("Failed to OBJ load file");

        return self.read_model(models);
    }

    fn read_model(&mut self, models: Vec<tobj::Model>) -> bool {
        for (_, m) in models.iter().enumerate() {
            let mesh = &m.mesh;
            // Read vertices
            for v in 0..mesh.positions.len() / 3 {
                let vx = mesh.positions[3 * v] as f64;
                let vy = mesh.positions[3 * v + 1] as f64;
                let vz = mesh.positions[3 * v + 2] as f64;
                self.add_point(Vector3D::new(vx, vy, vz));
            }

            // Read normals
            for v in 0..mesh.normals.len() / 3 {
                let vx = mesh.normals[3 * v] as f64;
                let vy = mesh.normals[3 * v + 1] as f64;
                let vz = mesh.normals[3 * v + 2] as f64;
                self.add_normal(Vector3D::new(vx, vy, vz));
            }

            // Read UVs
            for v in 0..mesh.texcoords.len() / 2 {
                let tu = mesh.texcoords[2 * v] as f64;
                let tv = mesh.texcoords[2 * v + 1] as f64;
                self.add_uv(Vector2D::new(tu, tv));
            }

            // the mesh consists only of triangles.
            let mut next_face = 0;
            if mesh.face_arities.is_empty() {
                let fv = 3;
                for _ in 0..mesh.indices.len() / 3 {
                    if !mesh.indices.is_empty() {
                        self.add_point_triangle(USize3::new(mesh.indices[next_face] as usize,
                                                            mesh.indices[next_face + 1] as usize,
                                                            mesh.indices[next_face + 2] as usize));
                    }

                    if !mesh.normal_indices.is_empty() {
                        self.add_normal_triangle(USize3::new(mesh.normal_indices[next_face] as usize,
                                                             mesh.normal_indices[next_face + 1] as usize,
                                                             mesh.normal_indices[next_face + 2] as usize));
                    }

                    if !mesh.texcoord_indices.is_empty() {
                        self.add_uv_triangle(USize3::new(mesh.texcoord_indices[next_face] as usize,
                                                         mesh.texcoord_indices[next_face + 1] as usize,
                                                         mesh.texcoord_indices[next_face + 2] as usize));
                    }
                    next_face += fv;
                }
            }
        }


        self.invalidate_cache();
        return true;
    }
}

impl Surface3 for TriangleMesh3 {
    /// ```
    /// use vox_geometry_rust::vector3::Vector3D;
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::unit_tests_utils::*;
    /// use std::io::Cursor;
    /// use vox_geometry_rust::surface3::Surface3;
    /// let mut buffer = Cursor::new(get_cube_tri_mesh3x3x3obj());
    ///
    /// let mut mesh = TriangleMesh3::new_default(None, None);
    /// mesh.read_obj_buffer(&mut buffer);
    ///
    /// let brute_force_search = |pt:&Vector3D| {
    ///     let mut min_dist2 = f64::MAX;
    ///     let mut result = Vector3D::new_default();
    ///     for i in 0..mesh.number_of_triangles() {
    ///         let tri = mesh.triangle(i);
    ///         let local_result = tri.closest_point(pt);
    ///         let local_dist2 = pt.distance_squared_to(local_result);
    ///         if local_dist2 < min_dist2 {
    ///             min_dist2 = local_dist2;
    ///             result = local_result;
    ///         }
    ///     }
    ///     return result;
    /// };
    ///
    /// let num_samples = get_number_of_sample_points3();
    /// for i in 0..num_samples {
    ///     let actual = mesh.closest_point(&Vector3D::new_lst(get_sample_points3()[i]));
    ///     let expected = brute_force_search(&Vector3D::new_lst(get_sample_points3()[i]));
    ///     assert_eq!(expected, actual);
    /// }
    /// ```
    fn closest_point_local(&self, other_point: &Vector3D) -> Vector3D {
        self.build_bvh();

        let mut distance_func = |tri_idx: &usize, pt: &Vector3D| {
            let tri = self.triangle(*tri_idx);
            return tri.closest_distance(pt);
        };

        let query_result = self._bvh.read().unwrap().nearest(other_point, &mut distance_func);
        match query_result.item {
            None => {
                panic!("No item found!");
            }
            Some(item) => {
                return self.triangle(item).closest_point(other_point);
            }
        }
    }

    /// ```
    /// use vox_geometry_rust::vector3::Vector3D;
    /// use vox_geometry_rust::bounding_box3::BoundingBox3D;
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::unit_tests_utils::*;
    /// use std::io::Cursor;
    /// use vox_geometry_rust::surface3::Surface3;
    /// let mut buffer = Cursor::new(get_cube_tri_mesh3x3x3obj());
    ///
    /// let mut mesh = TriangleMesh3::new_default(None, None);
    /// mesh.read_obj_buffer(&mut buffer);
    ///
    /// let aabb = BoundingBox3D::new(Vector3D::new(-0.5, -0.5, -0.5), Vector3D::new(0.5, 0.5, 0.5));
    /// assert_eq!(aabb.lower_corner, mesh.bounding_box().lower_corner);
    /// assert_eq!(aabb.upper_corner, mesh.bounding_box().upper_corner);
    /// ```
    fn bounding_box_local(&self) -> BoundingBox3D {
        self.build_bvh();

        return self._bvh.read().unwrap().bounding_box();
    }

    /// ```
    /// use vox_geometry_rust::vector3::Vector3D;
    /// use vox_geometry_rust::ray3::Ray3D;
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::unit_tests_utils::*;
    /// use std::io::Cursor;
    /// use vox_geometry_rust::surface3::*;
    /// use vox_geometry_rust::assert_delta;
    /// let mut buffer = Cursor::new(get_cube_tri_mesh3x3x3obj());
    ///
    /// let mut mesh = TriangleMesh3::new_default(None, None);
    /// mesh.read_obj_buffer(&mut buffer);
    ///
    /// let num_samples = get_number_of_sample_points3();
    ///
    /// let brute_force_test = |ray:&Ray3D| {
    ///     let mut result = SurfaceRayIntersection3::new();
    ///     for i in 0..mesh.number_of_triangles() {
    ///         let tri = mesh.triangle(i);
    ///         let local_result = tri.closest_intersection(ray);
    ///         if local_result.distance < result.distance {
    ///             result = local_result;
    ///         }
    ///     }
    ///     return result;
    /// };
    ///
    /// for i in 0..num_samples {
    ///     let ray = Ray3D::new(Vector3D::new_lst(get_sample_points3()[i]),
    ///                          Vector3D::new_lst(get_sample_dirs3()[i]));
    ///     let actual = mesh.closest_intersection(&ray);
    ///     let expected = brute_force_test(&ray);
    ///     assert_delta!(expected.distance, actual.distance, f64::EPSILON);
    ///     assert_eq!(expected.point.is_similar(&actual.point, None), true);
    ///     assert_eq!(expected.normal, actual.normal);
    ///     assert_eq!(expected.is_intersecting, actual.is_intersecting);
    /// }
    /// ```
    fn closest_intersection_local(&self, ray: &Ray3D) -> SurfaceRayIntersection3 {
        self.build_bvh();

        let mut test_func = |tri_idx: &usize, ray: &Ray3D| {
            let tri = self.triangle(*tri_idx);
            let result = tri.closest_intersection(ray);
            return result.distance;
        };

        let query_result = self._bvh.read().unwrap().closest_intersection(ray, &mut test_func);
        let mut result = SurfaceRayIntersection3::new();
        result.distance = query_result.distance;
        match query_result.item {
            None => {
                result.is_intersecting = false;
            }
            Some(item) => {
                result.is_intersecting = true;
                result.point = ray.point_at(query_result.distance);
                result.normal = self.triangle(item).closest_normal(&result.point);
            }
        }
        return result;
    }

    /// ```
    /// use vox_geometry_rust::vector3::Vector3D;
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::unit_tests_utils::*;
    /// use std::io::Cursor;
    /// use vox_geometry_rust::surface3::Surface3;
    /// let mut buffer = Cursor::new(get_sphere_tri_mesh5x5obj());
    ///
    /// let mut mesh = TriangleMesh3::new_default(None, None);
    /// mesh.read_obj_buffer(&mut buffer);
    ///
    /// let brute_force_search = |pt:&Vector3D| {
    ///     let mut min_dist2 = f64::MAX;
    ///     let mut result = Vector3D::new_default();
    ///     for i in 0..mesh.number_of_triangles() {
    ///         let tri = mesh.triangle(i);
    ///         let local_result = tri.closest_normal(pt);
    ///         let closest_pt = tri.closest_point(pt);
    ///         let local_dist2 = pt.distance_squared_to(closest_pt);
    ///         if local_dist2 < min_dist2 {
    ///             min_dist2 = local_dist2;
    ///             result = local_result;
    ///         }
    ///     }
    ///     return result;
    /// };
    ///
    /// let num_samples = get_number_of_sample_points3();
    /// for i in 0..num_samples {
    ///     let actual = mesh.closest_normal(&Vector3D::new_lst(get_sample_points3()[i]));
    ///     let expected = brute_force_search(&Vector3D::new_lst(get_sample_points3()[i]));
    ///     assert_eq!(expected.is_similar(&actual, Some(1.0e-9)), true);
    /// }
    /// ```
    fn closest_normal_local(&self, other_point: &Vector3D) -> Vector3D {
        self.build_bvh();

        let mut distance_func = |tri_idx: &usize, pt: &Vector3D| {
            let tri = self.triangle(*tri_idx);
            return tri.closest_distance(pt);
        };

        let query_result = self._bvh.read().unwrap().nearest(other_point, &mut distance_func);
        match query_result.item {
            None => {
                panic!("No item found!");
            }
            Some(item) => {
                return self.triangle(item).closest_normal(other_point);
            }
        }
    }

    /// ```
    /// use vox_geometry_rust::vector3::Vector3D;
    /// use vox_geometry_rust::ray3::Ray3D;
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::unit_tests_utils::*;
    /// use std::io::Cursor;
    /// use vox_geometry_rust::surface3::Surface3;
    /// let mut buffer = Cursor::new(get_cube_tri_mesh3x3x3obj());
    ///
    /// let mut mesh = TriangleMesh3::new_default(None, None);
    /// mesh.read_obj_buffer(&mut buffer);
    ///
    /// let num_samples = get_number_of_sample_points3();
    ///
    /// let brute_force_test = |ray:&Ray3D| {
    ///     for i in 0..mesh.number_of_triangles() {
    ///         let tri = mesh.triangle(i);
    ///         if tri.intersects(ray) {
    ///             return true;
    ///         }
    ///     }
    ///     return false;
    /// };
    ///
    /// for i in 0..num_samples {
    ///     let ray = Ray3D::new(Vector3D::new_lst(get_sample_points3()[i]),
    ///                          Vector3D::new_lst(get_sample_dirs3()[i]));
    ///     let actual = mesh.intersects(&ray);
    ///     let expected = brute_force_test(&ray);
    ///     assert_eq!(expected, actual);
    /// }
    /// ```
    fn intersects_local(&self, ray_local: &Ray3D) -> bool {
        self.build_bvh();

        let mut test_func = |tri_idx: &usize, ray: &Ray3D| {
            let tri = self.triangle(*tri_idx);
            return tri.intersects(ray);
        };

        return self._bvh.read().unwrap().intersects_ray(ray_local, &mut test_func);
    }

    /// ```
    /// use vox_geometry_rust::vector3::Vector3D;
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::unit_tests_utils::*;
    /// use std::io::Cursor;
    /// use vox_geometry_rust::surface3::Surface3;
    /// let mut buffer = Cursor::new(get_cube_tri_mesh3x3x3obj());
    ///
    /// let mut mesh = TriangleMesh3::new_default(None, None);
    /// mesh.read_obj_buffer(&mut buffer);
    ///
    /// let brute_force_search = |pt:&Vector3D| {
    ///     let mut min_dist = f64::MAX;
    ///     for i in 0..mesh.number_of_triangles() {
    ///         let tri = mesh.triangle(i);
    ///         let local_result = tri.closest_distance(pt);
    ///         if local_result < min_dist {
    ///             min_dist = local_result;
    ///         }
    ///     }
    ///     return min_dist;
    /// };
    ///
    /// let num_samples = get_number_of_sample_points3();
    /// for i in 0..num_samples {
    ///     let actual = mesh.closest_distance(&Vector3D::new_lst(get_sample_points3()[i]));
    ///     let expected = brute_force_search(&Vector3D::new_lst(get_sample_points3()[i]));
    ///     assert_eq!(expected, actual);
    /// }
    /// ```
    fn closest_distance_local(&self, other_point_local: &Vector3D) -> f64 {
        self.build_bvh();

        let mut distance_func = |tri_idx: &usize, pt: &Vector3D| {
            let tri = self.triangle(*tri_idx);
            return tri.closest_distance(pt);
        };

        let query_result = self._bvh.read().unwrap().nearest(other_point_local, &mut distance_func);
        return query_result.distance;
    }

    /// ```
    /// use vox_geometry_rust::vector3::Vector3D;
    /// use vox_geometry_rust::ray3::Ray3D;
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::unit_tests_utils::*;
    /// use std::io::Cursor;
    /// use vox_geometry_rust::surface3::Surface3;
    /// let mut buffer = Cursor::new(get_cube_tri_mesh3x3x3obj());
    ///
    /// let mut mesh = TriangleMesh3::new_default(None, None);
    /// mesh.read_obj_buffer(&mut buffer);
    ///
    /// let num_samples = get_number_of_sample_points3();
    ///
    /// for i in 0..num_samples {
    ///     let p = Vector3D::new_lst(get_sample_points3()[i]);
    ///     let actual = mesh.is_inside(&p);
    ///     let expected = mesh.bounding_box().contains(&p);
    ///     assert_eq!(expected, actual);
    /// }
    /// ```
    fn is_inside_local(&self, other_point_local: &Vector3D) -> bool {
        return self.fast_winding_number(other_point_local, K_DEFAULT_FAST_WINDING_NUMBER_ACCURACY) >
            0.5;
    }

    fn update_query_engine(&self) {
        self.build_bvh();
        self.build_winding_numbers();
    }

    fn view(&self) -> &Surface3Data {
        return &self.surface_data;
    }
}

/// Shared pointer for the TriangleMesh3 type.
pub type TriangleMesh3Ptr = Arc<RwLock<TriangleMesh3>>;

//--------------------------------------------------------------------------------------------------
///
/// # Front-end to create TriangleMesh3 objects step by step.
///
pub struct Builder {
    _points: Vec<Vector3D>,
    _normals: Vec<Vector3D>,
    _uvs: Vec<Vector2D>,
    _point_indices: Vec<USize3>,
    _normal_indices: Vec<USize3>,
    _uv_indices: Vec<USize3>,

    _surface_data: Surface3Data,
}

impl Builder {
    /// Returns builder with points.
    pub fn with_points(&mut self, points: Vec<Vector3D>) -> &mut Self {
        self._points = points;
        return self;
    }

    /// Returns builder with normals.
    pub fn with_normals(&mut self, normals: Vec<Vector3D>) -> &mut Self {
        self._normals = normals;
        return self;
    }

    /// Returns builder with uvs.
    pub fn with_uvs(&mut self, uvs: Vec<Vector2D>) -> &mut Self {
        self._uvs = uvs;
        return self;
    }

    /// Returns builder with point indices.
    pub fn with_point_indices(&mut self, point_indices: Vec<USize3>) -> &mut Self {
        self._point_indices = point_indices;
        return self;
    }

    /// Returns builder with normal indices.
    pub fn with_normal_indices(&mut self, normal_indices: Vec<USize3>) -> &mut Self {
        self._normal_indices = normal_indices;
        return self;
    }

    /// Returns builder with uv indices.
    pub fn with_uv_indices(&mut self, uv_indices: Vec<USize3>) -> &mut Self {
        self._uv_indices = uv_indices;
        return self;
    }

    /// Builds TriangleMesh3.
    /// ```
    /// use vox_geometry_rust::vector3::Vector3D;
    /// use vox_geometry_rust::vector2::Vector2D;
    /// use vox_geometry_rust::triangle_mesh3::TriangleMesh3;
    /// use vox_geometry_rust::usize3::USize3;
    ///
    /// let points = vec![
    ///     Vector3D::new(1.0, 2.0, 3.0),
    ///     Vector3D::new(4.0, 5.0, 6.0),
    ///     Vector3D::new(7.0, 8.0, 9.0),
    ///     Vector3D::new(10.0, 11.0, 12.0)
    /// ];
    ///
    /// let normals = vec![
    ///     Vector3D::new(10.0, 11.0, 12.0),
    ///     Vector3D::new(7.0, 8.0, 9.0),
    ///     Vector3D::new(4.0, 5.0, 6.0),
    ///     Vector3D::new(1.0, 2.0, 3.0)
    /// ];
    ///
    /// let uvs = vec![
    ///     Vector2D::new(13.0, 14.0),
    ///     Vector2D::new(15.0, 16.0)
    /// ];
    ///
    /// let point_indices = vec![
    ///     USize3::new(0, 1, 2),
    ///     USize3::new(0, 1, 3)
    /// ];
    ///
    /// let normal_indices = vec![
    ///     USize3::new(1, 2, 3),
    ///     USize3::new(2, 1, 0)
    /// ];
    ///
    /// let uv_indices = vec![
    ///     USize3::new(1, 0, 2),
    ///     USize3::new(3, 1, 0)
    /// ];
    ///
    /// let mesh = TriangleMesh3::builder()
    ///     .with_points(points.clone())
    ///     .with_normals(normals.clone())
    ///     .with_uvs(uvs.clone())
    ///     .with_point_indices(point_indices.clone())
    ///     .with_normal_indices(normal_indices.clone())
    ///     .with_uv_indices(uv_indices.clone())
    ///     .build();
    ///
    /// assert_eq!(4, mesh.number_of_points());
    /// assert_eq!(4, mesh.number_of_normals());
    /// assert_eq!(2, mesh.number_of_uvs());
    /// assert_eq!(2, mesh.number_of_triangles());
    ///
    /// for i in 0..mesh.number_of_points() {
    ///     assert_eq!(points[i], mesh.point(i));
    /// }
    ///
    /// for i in 0..mesh.number_of_normals() {
    ///     assert_eq!(normals[i], mesh.normal(i));
    /// }
    ///
    /// for i in 0..mesh.number_of_uvs() {
    ///     assert_eq!(uvs[i], mesh.uv(i));
    /// }
    ///
    /// for i in 0..mesh.number_of_triangles() {
    ///     assert_eq!(point_indices[i], mesh.point_index(i));
    ///     assert_eq!(normal_indices[i], mesh.normal_index(i));
    ///     assert_eq!(uv_indices[i], mesh.uv_index(i));
    /// }
    /// ```
    pub fn build(&mut self) -> TriangleMesh3 {
        return TriangleMesh3::new(self._points.clone(),
                                  self._normals.clone(),
                                  self._uvs.clone(),
                                  self._point_indices.clone(),
                                  self._normal_indices.clone(),
                                  self._uv_indices.clone(),
                                  Some(self._surface_data.transform.clone()),
                                  Some(self._surface_data.is_normal_flipped));
    }

    /// Builds shared pointer of TriangleMesh3 instance.
    pub fn make_shared(&mut self) -> TriangleMesh3Ptr {
        return TriangleMesh3Ptr::new(RwLock::new(self.build()));
    }

    /// constructor
    pub fn new() -> Builder {
        return Builder {
            _points: vec![],
            _normals: vec![],
            _uvs: vec![],
            _point_indices: vec![],
            _normal_indices: vec![],
            _surface_data: Surface3Data::new(None, None),
            _uv_indices: vec![],
        };
    }
}

impl SurfaceBuilderBase3 for Builder {
    fn view(&mut self) -> &mut Surface3Data {
        return &mut self._surface_data;
    }
}