cg-math 0.1.2

extern crate nalgebra as na;
extern crate parry2d_f64 as parry;

use crate::geometry::{WeightedEdge, WeightedTriangle, WeightedVertex, Triangulation, TriangleLogic};
use na::Point2;
use rand::{distributions::Uniform, rngs::ThreadRng, thread_rng, Rng};
use std::collections::HashMap;

/// Sample n vertices (i.e. w_i = 0.0), with coordinates drawn form a Uniform distribution
#[allow(dead_code)]
pub fn generate_vertices(n: u32) -> Vec<WeightedVertex> {
    let mut rng: ThreadRng = thread_rng();
    let distribution: Uniform<f64> = Uniform::new(0.0, 10.0);

    (0..n)
        .map(|_| WeightedVertex::new(rng.sample(distribution), rng.sample(distribution), 0.0))
        .collect()
}

/// Sample n weighted vertices, with coordinates and weights drawn form a Uniform distribution
#[allow(dead_code)]
pub fn generate_weighted_vertices(n: u32) -> Vec<WeightedVertex> {
    let mut rng: ThreadRng = thread_rng();
    let distribution: Uniform<f64> = Uniform::new(0.0, 10.0);

    (0..n)
        .map(|_| {
            WeightedVertex::new(
                rng.sample(distribution),
                rng.sample(distribution),
                rng.sample(distribution),
            )
        })
        .collect()
}

/// Given an edge of the triangle and, retrieve the point not on the edge
#[allow(dead_code)]
pub fn get_fourth_point(triangle: &WeightedTriangle, edge: &WeightedEdge) -> WeightedVertex {
    if triangle.a != edge.a && triangle.a != edge.b {
        WeightedVertex::new(triangle.a.x(), triangle.a.y(), triangle.a.w())
    } else if triangle.b != edge.a && triangle.b != edge.b {
        WeightedVertex::new(triangle.b.x(), triangle.b.y(), triangle.b.w())
    } else {
        WeightedVertex::new(triangle.c.x(), triangle.c.y(), triangle.c.w())
    }
}

/// Compute a points barycentrics coordinates w.r.t. the given triangle
pub fn barycentric_coords(vertex: &WeightedVertex, triangle: &WeightedTriangle) -> (f64, f64, f64) {
    // adapted from: https://gamedev.stackexchange.com/questions/23743/whats-the-most-efficient-way-to-find-barycentric-coordinates
    let v0 = triangle.b.as_vec2() - triangle.a.as_vec2();
    let v1 = triangle.c.as_vec2() - triangle.a.as_vec2();
    let v2 = vertex.as_vec2() - triangle.a.as_vec2();

    let d00 = v0.dot(&v0);
    let d01 = v0.dot(&v1);
    let d11 = v1.dot(&v1);
    let d20 = v2.dot(&v0);
    let d21 = v2.dot(&v1);

    let inv_denom = 1. / (d00 * d11 - d01 * d01); // less divisions for optimization

    let v = (d11 * d20 - d01 * d21) * inv_denom;
    let w = (d00 * d21 - d01 * d20) * inv_denom;
    let u = 1. - v - w;

    assert!(v + w + u - 1. < 0.0001); // sanity check

    (v, w, u)
}

/// Check whether the vertex lies inside the given triangle
pub fn point_in_tri(vertex: &WeightedVertex, triangle: &WeightedTriangle) -> bool {
    let (v, w, u) = barycentric_coords(vertex, triangle);
    v >= 0.0 && w >= 0.0 && u >= 0.0
}

/// Compute the convex hull of the given points in 2D
pub fn c_hull(points: &Vec<WeightedVertex>) -> Vec<WeightedVertex> {
    let points2d: Vec<Point2<f64>> = points.iter().map(|v| Point2::new(v.x(), v.y())).collect();
    let c_hull = parry::transformation::convex_hull(&points2d);
    
    c_hull
        .iter()
        .map(|p| {
            // since order is not invariant, we need to map the weights back to the correct vertices.
            for v in points {
                if v.x() == p.x && v.y() == p.y {
                    return WeightedVertex::new(p.x, p.y, v.w());
                }
            }
            panic!("One of the vertices should match the convex hull vertices.")
        })
        .collect()
}

/// A 1-simplex and its two incident 2-simplices yield four points, of which all 2-simplices of the given triangulation that
/// consist of any permutation of 3 of these vertices form the induced subcomplex.
/// 
/// More formally: <br>
/// Let e be and edge of an arbitrary triangulation T and t1, t2 be it's incident triangles (given they exist). <br>
/// The induced subcomplex of the four vertices of e, t1, t2, namely V, consists of all simplices in T <br>
/// spanned by vertices in V.
/// 
/// ```
/// use cg_math::{geometry::{WeightedVertex, WeightedTriangle, Triangulation}, utils::induced_subcomplex};
/// 
/// // these 4 vertices span two of the three triangles in the triangulation
/// let v0 = WeightedVertex::new(6.24, 8.38, 2.24);
/// let v1 = WeightedVertex::new(9.34, 9.19, 7.46);
/// let v2 = WeightedVertex::new(8.04, 5.53, 7.53);
/// let v3 = WeightedVertex::new(8.53, 8.43, 1.34);
/// 
/// let mut triangulation = Triangulation::empty();
/// triangulation.triangles_mut().append(&mut vec![
///     WeightedTriangle::new(v0, v1, v3),
///     WeightedTriangle::new(v2 ,v1, v3),
///     WeightedTriangle::new(v2 ,v1, v0),
/// ]);
/// 
/// let mut ics = induced_subcomplex(&vec![v0, v1, v2, v3], &triangulation);
/// ics.sort();
/// triangulation.triangles_mut().sort();
/// 
/// // Yet their induced subcomplex consists of all three triangles.
/// assert!(ics.len() == 3);
/// assert_eq!(ics, *triangulation.triangles());
/// ```
pub fn induced_subcomplex(vertices: &[WeightedVertex], triangulation: &Triangulation) -> Vec<WeightedTriangle> {
    let subcomplex: Vec<WeightedTriangle> = triangulation.triangles().iter().filter(|tri| {
        tri.has_vertex(&vertices[0]) && tri.has_vertex(&vertices[1]) && tri.has_vertex(&vertices[2]) 
        || tri.has_vertex(&vertices[0]) && tri.has_vertex(&vertices[1]) && tri.has_vertex(&vertices[3])
        || tri.has_vertex(&vertices[1]) && tri.has_vertex(&vertices[2]) && tri.has_vertex(&vertices[3])
        || tri.has_vertex(&vertices[0]) && tri.has_vertex(&vertices[2]) && tri.has_vertex(&vertices[3])
    }).copied().collect();
    subcomplex
}


/// An edge is flippable if the pointwise union of its induced subcomplex equals the convex hull of all the vertices in the induced subcomplex.
/// ```
/// use cg_math::{geometry::{WeightedVertex, WeightedTriangle, Triangulation}, utils::{induced_subcomplex, is_flippable2}};
/// 
/// // these 4 vertices span two that would need to be 3-to-1 flipped
/// let v0 = WeightedVertex::new(6.24, 8.38, 2.24); 
/// let v1 = WeightedVertex::new(9.34, 9.19, 7.46);
/// let v2 = WeightedVertex::new(8.04, 5.53, 7.53);
/// let v3 = WeightedVertex::new(8.53, 8.43, 1.34);
/// 
/// let mut triangulation = Triangulation::empty();
/// triangulation.triangles_mut().append(&mut vec![
///     WeightedTriangle::new(v0, v1, v3),
///     WeightedTriangle::new(v2 ,v1, v3),
///     WeightedTriangle::new(v2 ,v1, v0),
/// ]);
/// 
/// let mut ics = induced_subcomplex(&vec![v0, v1, v2, v3], &triangulation);
/// for tri in ics.iter() {
///     println!("{}", tri);
///  }
/// assert_eq!(is_flippable2(&ics), true);
/// ```
pub fn  is_flippable_strict(induced_subcomplex: &[WeightedTriangle]) -> bool {
    let points = induced_subcomplex
        .iter()
        .flat_map(|t| [Point2::new(t.a.x(), t.a.y()), Point2::new(t.b.x(), t.b.y()), Point2::new(t.c.x(), t.c.y())] )
        .collect::<Vec<Point2<f64>>>();

    let c_hull = parry::transformation::convex_hull(&points);
    // is even faster for low number of items, than hashing, s.: https://stackoverflow.com/questions/64226562/check-if-vec-contains-all-elements-from-another-vec
    points.iter().all(|item| c_hull.contains(item))
}

/// Check whether the edge is flippable in the ploygon described all four points. <br>
/// "*We call T and e flippable if conv(T) is the underlying space of the induced subcomplex of T.*"
pub fn is_flippable(point: &WeightedVertex, edge: &WeightedEdge, pk: &WeightedVertex) -> bool {
    // we could also check if the polygon is convex
    // hopefully more correct parry version
    let points: Vec<Point2<f64>> = vec![    
        Point2::new(point.x(), point.y()),
        Point2::new(edge.a.x(), edge.a.y()),
        Point2::new(edge.b.x(), edge.b.y()),
        Point2::new(pk.x(), pk.y()),
    ];

    //let points2d: Vec<Point2<f64>> = points.iter().map(|v| Point2::new(v.x, v.y)).collect();
    let c_hull = parry::transformation::convex_hull(&points);

    //println!("points: {:?}", points);
    //println!("c_hull: {:?}", c_hull);

    // is even faster for low number of items, than hashing, s.: https://stackoverflow.com/questions/64226562/check-if-vec-contains-all-elements-from-another-vec
    points.iter().all(|item| c_hull.contains(item))
}

/// Keep the elements that are unique in the given vector of elements, i.e. that appear exactly once.
pub fn keep_unique_elements<T: Eq + std::hash::Hash>(elements: Vec<T>) -> Vec<T> {
    elements
        .into_iter()
        .fold(HashMap::new(), |mut map, element| {
            *map.entry(element).or_insert(0) += 1;
            map
        })
        .into_iter()
        .filter_map(|(element, count)| (count == 1).then_some(element))
        .collect()
}

/// Older, less performant version than `keep_unique_elements`
/// as found out by Tbolt on discord in rust-questions-and-answers-2
#[allow(dead_code)]
pub fn keep_unique_elements2<T: Eq + std::hash::Hash + Clone>(elements: &mut Vec<T>) {
    let mut freq_map = HashMap::new();

    // count frequency of each element
    for element in elements.iter() {
        *freq_map.entry(element.clone()).or_insert(0) += 1;
    }

    elements.retain(|element| freq_map.get(element) == Some(&1));
}

// Todo: can we do this more simple / efficient ? -> yes this is O(n^2), O(n) version just above
/// Deprecated: more efficient and readable version is `keep_unique_elements`; this is kept for learning purposes <br>
/// Remove duplicates from a list of elements of type T.
/// So keeping only unique elemtents
#[allow(warnings)]
pub fn remove_duplicates<T: std::cmp::PartialEq<T>>(elements: &mut Vec<T>) {
    let mut indices_to_remove: Vec<usize> = vec![];

    for i in 0..elements.len() {
        let element: &T = &elements[i];
        for j in 0..elements.len() {
            let comp_element: &T = &elements[j];
            if element == comp_element && i != j {
                indices_to_remove.push(i);
                break;
            }
        }
    }
    // since idx has elements in ascending order, e.g. [2, 7, 13]
    // we can ierate in reverse and not conflic any array boundaries
    for idx in indices_to_remove.iter().rev() {
        elements.remove(*idx);
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use na::{Point2, Vector2};

    #[test]
    fn test_barycentric_coords() {
        let tri = WeightedTriangle::new(
            WeightedVertex::new(0.0, 0.0, 0.0),
            WeightedVertex::new(1.0, 0.0, 0.0),
            WeightedVertex::new(0.0, 1.0, 0.0),
        );

        let point = WeightedVertex::new(0.5, 0.5, 0.0);
        let (v, w, u) = barycentric_coords(&point, &tri);
        assert_eq!(v, 0.5);
        assert_eq!(w, 0.5);
        assert_eq!(u, 0.0);

        let point = WeightedVertex::new(0.0, 0.0, 0.0);
        let (v, w, u) = barycentric_coords(&point, &tri);
        assert_eq!(v, 0.0);
        assert_eq!(w, 0.0);
        assert_eq!(u, 1.0);

        let point = WeightedVertex::new(1.0, 0.0, 0.0);
        let (v, w, u) = barycentric_coords(&point, &tri);
        assert_eq!(v, 0.0);
        assert_eq!(w, 1.0);
        assert_eq!(u, 0.0);

        let point = WeightedVertex::new(0.0, 1.0, 0.0);
        let (v, w, u) = barycentric_coords(&point, &tri);
        assert_eq!(v, 1.0);
        assert_eq!(w, 0.0);
        assert_eq!(u, 0.0);
    }

    #[test]
    fn test_convex_hull() {
        let points = vec![
            Vector2::new(0.5, 0.5),
            Vector2::new(2.5, 2.0),
            Vector2::new(1.0, 2.5),
            Vector2::new(2.0, 3.0),
            Vector2::new(1.0, 1.0),
        ];

        let true_hull = vec![
            Point2::new(0.5, 0.5),
            Point2::new(2.5, 2.0),
            Point2::new(2.0, 3.0),
            Point2::new(1.0, 2.5),
        ];

        let points2d: Vec<Point2<f64>> = points.iter().map(|v| Point2::new(v.x, v.y)).collect();
        let c_hull = parry::transformation::convex_hull(&points2d);
        //println!("points: {:?}", points2d);
        //println!("c_hull: {:?}", c_hull);
        //println!("true_hull: {:?}", true_hull);
        assert!(c_hull != points2d);
        assert!(true_hull.iter().all(|item| true_hull.contains(item)));
    }

    #[test]
    fn test_edge_flippable() {
        pub const LARGE_NUM: f64 = 10000.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);

        // here the edge in question is the one between (2.5, 2.0) and (1.0, 2.5)
        // the two triangles are t1 = {(2.5, 2.0), (0.5, 0.5), (1.0, 2.5)} and t2 = {(2.5, 2.0),(1.0, 2.5), (2.0, 3.0)}
        assert!(is_flippable(
            &WeightedVertex::new(0.5, 0.5, 0.0),
            &WeightedEdge::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)
        ));

        // this is acutally flippable, try with low numbers for L and if necessary make a sketch
        assert!(is_flippable(
            &WeightedVertex::new(2.5, 2.0, 0.0),
            &WeightedEdge::new(WeightedVertex::new(2.0, 3.0, 0.0), p_2_lw),
            &p_1_lw
        ));
    }

    #[test]
    fn test_edge_not_flippable() {
        pub const LARGE_NUM: f64 = 10000.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);

        assert!(!is_flippable(
            &WeightedVertex::new(-5., 1., 0.0),
            &WeightedEdge::new(
                WeightedVertex::new(-4., 2., 0.0),
                WeightedVertex::new(-2.5, 2., 0.0),
            ),
            &WeightedVertex::new(-5., 3., 0.0)
        ));

        assert!(!is_flippable(
            &WeightedVertex::new(-1.5, 0., 0.0),
            &WeightedEdge::new(WeightedVertex::new(0., 1.5, 0.0), p_2_lw),
            &WeightedVertex::new(2., -0.5, 0.0)
        ));
    }
}