cg-math 0.1.2

use crate::geometry::{WeightedEdge, WeightedTriangle, WeightedVertex};

fn setup() -> WeightedTriangle {
    pub const LARGE_NUM: f64 = 100.0; // NOTE: this parameter highly influences the convex hull property of the triangulation

    let p_2_lw: WeightedVertex = WeightedVertex::new(0.0, LARGE_NUM, 0.0);

    let p_1_lw: WeightedVertex = WeightedVertex::new(LARGE_NUM, -LARGE_NUM, 0.0);

    let p_0_lw: WeightedVertex = WeightedVertex::new(-LARGE_NUM, -LARGE_NUM, 0.0);

    let _e_21_lw: WeightedEdge = WeightedEdge::new(p_2_lw, p_1_lw);

    let _e_10_lw: WeightedEdge = WeightedEdge::new(p_1_lw, p_0_lw);

    let _e_02_lw: WeightedEdge = WeightedEdge::new(p_0_lw, p_2_lw);

    WeightedTriangle::new(p_2_lw, p_1_lw, p_0_lw)
}

/// **R**andomized **I**ncremental **C**onstruction algorithm for (weighted) delaunay triangulations.
pub mod ric {
    use crate::{
        algos::setup,
        geometry::{TriangleLogic, WeightedEdge, WeightedTriangle, WeightedVertex, Triangulation},
        utils::{point_in_tri, c_hull, is_flippable_strict, induced_subcomplex},
    };
    use rand::{seq::SliceRandom, thread_rng};

    /// Given a triangulation and a vertex, find the idx and triangle that contains this vertex. <br>
    /// Also removes the conflicting triangle from the triangulation.
    fn find_containing_triangle(
        triangulation: &mut [WeightedTriangle],
        vertex: &WeightedVertex,
    ) -> Option<WeightedTriangle> {
        triangulation
            .iter()
            .find(|triangle| point_in_tri(vertex, triangle))
            .copied()
    }

    /// Given an *edge* and a *vertex*, i.e. a triangle, find the optional neighboring triangle. <br>
    /// Return the neighboring triangle and the point (_pk_) of the triangle, that is **not** on the edge.
    fn find_neighboring_triangle(
        edge: &WeightedEdge,
        vertex: &WeightedVertex,
        triangulation: &[WeightedTriangle],
    ) -> Option<(WeightedTriangle, WeightedVertex)> {
        triangulation.iter().find_map(|triangle| {
            if triangle.has_edge(edge) && !triangle.has_vertex(vertex) {
                Some((*triangle, triangle.get_third_point(edge)))
            } else {
                None
            }
        })
    }

    /// Creates three new triangles from a given triangle and a point that lies within.
    fn one_to_three_split(triangle: &WeightedTriangle, vertex: &WeightedVertex) -> Vec<WeightedTriangle> {
        triangle
            .edges()
            .iter()
            .map(|edge| WeightedTriangle::new(*vertex, edge.a, edge.b))
            .collect()
    }

    /// Creates a new triangle from the four given points, which form a concave quadrilateral.
    fn three_to_one_split(points: &[WeightedVertex; 4]) -> WeightedTriangle {
        let c_hull = c_hull(&points.to_vec());
        WeightedTriangle::new(c_hull[0], c_hull[1], c_hull[2])
    }

    /// Given and edge and a vertex, i.e. a triangle, and a neighboring triangle, i.e. they share the edge <br>
    /// Return the two triangles two remove & add from / to the triangulation by flipping the edge, as (to remove, to add) <br>
    /// *pk* is the neighbors point that is not on the common edge
    fn two_flip(
        edge: &WeightedEdge,
        vertex: &WeightedVertex,
        neighbor_triangle: &WeightedTriangle,
    ) -> ([WeightedTriangle; 2], Vec<WeightedTriangle>) {
        let pk = neighbor_triangle.get_third_point(edge);
        (
            [
                *neighbor_triangle,
                WeightedTriangle::new(edge.a, edge.b, *vertex),
            ],
            vec![
                WeightedTriangle::new(edge.a, *vertex, pk),
                WeightedTriangle::new(edge.b, *vertex, pk),
            ],
        )
    }

    /// Compute the delaunay triangulation for a point set in $\mathbb{R}^2$ by the method of randomized incremental construction.
    /// 
    /// # Example
    /// 
    /// ```
    /// use cg_math::{algos::ric, geometry::{WeightedVertex, WeightedTriangle, Triangulation}};
    /// 
    /// let v0 = WeightedVertex::new(2.5, 2.0, 0.0);
    /// let v1 = WeightedVertex::new(1.0, 2.5, 0.0);
    /// let v2 = WeightedVertex::new(2.0, 3.0, 0.0);
    /// 
    /// let triangulation = ric::compute2d(&mut vec![v0, v1, v2], 0.0);
    /// 
    /// assert_eq!(triangulation.triangles().len(), 1);
    /// assert_eq!(triangulation.triangles()[0], WeightedTriangle::new(v0, v1, v2));
    /// ```
    pub fn compute2d(vertices: &mut Vec<WeightedVertex>, epsilon: f64) -> Triangulation {
        let super_tri = setup();
        let mut rng = thread_rng();
        let mut triangulation = Triangulation::new(super_tri);
        triangulation.cache.push(vec![super_tri]);
        let mut edge_stack: Vec<WeightedEdge> = Vec::new();

        vertices.shuffle(&mut rng);

        for vertex in vertices.iter() {
            if let Some(invalid_triangle) = find_containing_triangle(triangulation.triangles_mut(), vertex) {
                if !invalid_triangle.is_eps_regular(vertex, epsilon) {
                    // if the point is not regular w.r.t to this triangle: remove invalid triangle & connect by three flip
                    triangulation.used_vertices_mut().push(*vertex);
                    triangulation.triangles_mut().retain(|tri| *tri != invalid_triangle);
                    triangulation.triangles_mut().append(&mut one_to_three_split(&invalid_triangle, vertex));
                    triangulation.cache.push(triangulation.triangles().clone());
                    edge_stack.append(&mut invalid_triangle.edges()); // append possibly illegal edges to the stack

                    while !edge_stack.is_empty() {
                        if let Some(edge) = edge_stack.pop() {
                            if let Some((neighbor_triangle, pk)) =
                                find_neighboring_triangle(&edge, vertex, triangulation.triangles())
                            {   
                                if !neighbor_triangle.is_regular(vertex) {   
                                    // get the induced subcomplex of the 4 points, i.e. all triangles that contain 3 of these points
                                    let subcomplex: Vec<WeightedTriangle> = induced_subcomplex(&[edge.a, edge.b, *vertex, pk], &triangulation);

                                    if subcomplex.len() == 2 && is_flippable_strict(&subcomplex) { // is_flippable_strict(&subcomplex) is_flippable(vertex, &edge, &pk)
                                        // Perform the 2 flip
                                        let (triangles_to_remove, mut triangles_to_add) =
                                            two_flip(&edge, vertex, &neighbor_triangle);
                                        triangulation.triangles_mut().retain(|tri| !triangles_to_remove.contains(tri));
                                        triangulation.triangles_mut().append(&mut triangles_to_add);
                                        triangulation.cache.push(triangulation.triangles().clone());
    
                                        // update edge stack
                                        edge_stack.append(&mut vec![
                                            WeightedEdge::new(edge.a, pk),
                                            WeightedEdge::new(edge.b, pk),
                                        ])
                                    } else if subcomplex.len() == 3 {
                                        let triangle = three_to_one_split(&[edge.a, edge.b, *vertex, pk]);
                                        // remove the 3 triangles from the triangulation
                                        triangulation.triangles_mut().retain(|tri| !subcomplex.contains(tri));
                                        triangulation.triangles_mut().push(triangle);
                                        edge_stack.append(&mut triangle.edges());
                                    }
                                }
                            }
                        }
                    }
                } else {
                    triangulation.ignored_vertices_mut().push(*vertex);
                }
            } else {
                panic!("The vertex should lie in one of the triangles!")
            }
            triangulation.cache.push(triangulation.triangles().clone());
        }
        // discard_super_triangle(&mut triangulation, &p_0, &p_1);
        triangulation.triangles_mut().retain(|triangle| {
            !triangle.has_vertex(&super_tri.a)
                && !triangle.has_vertex(&super_tri.b)
                && !triangle.has_vertex(&super_tri.c)
        });
        
        triangulation.cache.push(triangulation.triangles().clone());
    
        triangulation
    }
}

/// **B**owyer & **W**atson algorithm for (weighted) delaunay triangulations.
pub mod bw {

    extern crate nalgebra as na;

    use rand::seq::SliceRandom;
    use rand::thread_rng;

    use crate::geometry::{TriangleLogic, Triangulation};
    use crate::geometry::{WeightedEdge, WeightedTriangle, WeightedVertex};
    use crate::utils::{keep_unique_elements, point_in_tri};

    use super::setup;

    /// Compute the delaunay triangulation for a point set in $\mathbb{R}^2$ by the method of Bowyer & Watson.
    /// 
    /// # Example
    /// 
    /// ```
    /// use cg_math::{algos::bw, geometry::{WeightedVertex, WeightedTriangle, Triangulation}};
    /// 
    /// let v0 = WeightedVertex::new(2.5, 2.0, 0.0);
    /// let v1 = WeightedVertex::new(1.0, 2.5, 0.0);
    /// let v2 = WeightedVertex::new(2.0, 3.0, 0.0);
    /// 
    /// let triangulation = bw::compute2d(&mut vec![v0, v1, v2], 0.0);
    /// 
    /// assert_eq!(triangulation.triangles().len(), 1);
    /// assert_eq!(triangulation.triangles()[0], WeightedTriangle::new(v0, v1, v2));
    /// ```
    pub fn compute2d(vertices: &mut Vec<WeightedVertex>, epsilon: f64) -> Triangulation {
        // 1 initialize with super triangle
        let super_tri = setup();
        let mut triangulation = Triangulation::new(super_tri);
        triangulation.cache.push(vec![super_tri]);
        let mut rng = thread_rng();

        vertices.shuffle(&mut rng);


        // 2 triangulate each vertex
        for vertex in vertices {
            add_vertex(&mut triangulation, vertex, epsilon);
        }

        // 3 only keep the triangles that are not connected to the super triangle
        triangulation.triangles_mut().retain(|triangle| {
            !triangle.has_vertex(&super_tri.a)
                && !triangle.has_vertex(&super_tri.b)
                && !triangle.has_vertex(&super_tri.c)
        });
        triangulation.cache.push(triangulation.triangles().clone());
        triangulation
    }

    /// Add a vertex to the triangulation.
    fn add_vertex(
        triangulation: &mut Triangulation,
        vertex: &WeightedVertex,    
        epsilon: f64,
    ) {
        let mut edges_to_retriangulate: Vec<WeightedEdge> = Vec::new();

        // 2.x find triangle that contains the vertex and if e-regular just return
        // with DAG we can do this in n log n
        for triangle in triangulation.triangles().iter() {
            if point_in_tri(vertex, triangle) {
                if triangle.is_eps_regular(vertex, epsilon) {
                    triangulation.ignored_vertices_mut().push(*vertex);
                    return;
                }
                break;
            }
        }
        triangulation.used_vertices_mut().push(*vertex);

        // 2.1 keep only the triangles that don't contain the point, for the other trianlges extract their edges
        triangulation.triangles_mut().retain(|triangle| {
            let regular = triangle.is_regular(vertex);
            if !regular {
                edges_to_retriangulate.append(&mut triangle.edges());
            }
            regular
        });

        // 2.2 keep only unique edges -> the ones appearing exactly one time, all other are discarded
        edges_to_retriangulate = keep_unique_elements(edges_to_retriangulate);

        // 2.3 to fill the polygonal hole produced in step 2.1, build new triangles with the current point and the unique edges
        for edge in edges_to_retriangulate {
            triangulation.triangles_mut().push(WeightedTriangle::new(edge.a, edge.b, *vertex));
        }
        triangulation.cache.push(triangulation.triangles().clone());
    }
}

#[cfg(test)]
mod tests {
    use crate::{
        algos::{bw, ric},
        geometry::{WeightedTriangle, WeightedVertex, WeightedTriangulation}, validation::all_globally_regular,
    };

    #[test]
    fn test_algos() {
        let vertices = vec![
            WeightedVertex::new(2.0, 3.0, 0.0),
            WeightedVertex::new(2.5, 2.0, 0.0),
            WeightedVertex::new(1.0, 2.5, 0.0),
            WeightedVertex::new(0.5, 0.5, 0.0),
        ];
        let mut true_triangulation = vec![
            WeightedTriangle::new(
                WeightedVertex::new(2.5, 2.0, 0.0),
                WeightedVertex::new(1.0, 2.5, 0.0),
                WeightedVertex::new(2.0, 3.0, 0.0),
            ),
            WeightedTriangle::new(
                WeightedVertex::new(2.5, 2.0, 0.0),
                WeightedVertex::new(1.0, 2.5, 0.0),
                WeightedVertex::new(0.5, 0.5, 0.0),
            ),
        ];
        let eps = 0.0;
        let mut triangulation_bw = bw::compute2d(&mut vertices.clone(), eps);
        let mut triangulation_ric = ric::compute2d(&mut vertices.clone(), eps);

        assert_eq!(true_triangulation.len(), triangulation_bw.triangles().len());
        assert_eq!(true_triangulation.len(), triangulation_ric.triangles().len());

        true_triangulation.sort();
        triangulation_bw.triangles_mut().sort();
        triangulation_ric.triangles_mut().sort();
        assert_eq!(true_triangulation, *triangulation_bw.triangles());
        assert_eq!(true_triangulation, *triangulation_ric.triangles());
    }

    #[test]
    fn test_debug_ric() {
        let vertices = vec![
            WeightedVertex::new(5.27, 1.20, 2.95),
            WeightedVertex::new(3.35, 7.15, 4.69),
            WeightedVertex::new(6.13, 6.62, 1.76),
            WeightedVertex::new(4.56, 8.19, 1.13),
            WeightedVertex::new(7.44, 5.76, 8.45),
            WeightedVertex::new(3.56, 7.33, 7.93),
        ];

        let triangulation = ric::compute2d(&mut vertices.clone(), 0.0);
        let regularity = all_globally_regular(&vertices, &triangulation.triangles());
        println!("\n\nFinal triangulation:\n{}", WeightedTriangulation(triangulation.triangles()));
        println!("Regularity: {}", regularity);
        println!();
        for tris in triangulation.cache {
            println!("\n\nTriangulation:\n{}", WeightedTriangulation(&tris));
        }
    }
}