1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

#![allow(dead_code)] // TODO: remove this

use crate::mass_properties::MassProperties;
use crate::math::{Point, Real};
use crate::shape::Triangle;

impl MassProperties {
    /// Computes the mass properties of a convex polygon.
    pub fn from_convex_polygon(density: Real, vertices: &[Point<Real>]) -> MassProperties {
        let (area, com) = convex_polygon_area_and_center_of_mass(vertices);

        if area == 0.0 {
            return MassProperties::new(com, 0.0, 0.0);
        }

        let mut itot = 0.0;

        let mut iterpeek = vertices.iter().peekable();
        let first_element = *iterpeek.peek().unwrap(); // store first element to close the cycle in the end with unwrap_or
        while let Some(elem) = iterpeek.next() {
            let triangle = Triangle::new(com, *elem, **iterpeek.peek().unwrap_or(&first_element));
            let area = triangle.area();
            let ipart = triangle.unit_angular_inertia();
            itot += ipart * area;
        }

        Self::new(com, area * density, itot * density)
    }
}

/// Computes the area and center-of-mass of a convex polygon.
pub fn convex_polygon_area_and_center_of_mass(
    convex_polygon: &[Point<Real>],
) -> (Real, Point<Real>) {
    let geometric_center = convex_polygon
        .iter()
        .fold(Point::origin(), |e1, e2| e1 + e2.coords)
        / convex_polygon.len() as Real;
    let mut res = Point::origin();
    let mut areasum = 0.0;

    let mut iterpeek = convex_polygon.iter().peekable();
    let first_element = *iterpeek.peek().unwrap(); // Stores first element to close the cycle in the end with unwrap_or.
    while let Some(elem) = iterpeek.next() {
        let (a, b, c) = (
            elem,
            iterpeek.peek().unwrap_or(&first_element),
            &geometric_center,
        );
        let area = Triangle::new(*a, **b, *c).area();
        let center = (a.coords + b.coords + c.coords) / 3.0;

        res += center * area;
        areasum += area;
    }

    if areasum == 0.0 {
        (areasum, geometric_center)
    } else {
        (areasum, res / areasum)
    }
}