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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

//! methods for 2D poly loop

use num_traits::AsPrimitive;

/// area
pub fn area<T>(
    vtx2xy: &[T]) -> T
    where T: num_traits::Float + Copy + 'static + std::ops::AddAssign,
          f64: AsPrimitive<T>
{
    let num_vtx = vtx2xy.len() / 2;
    assert_eq!(vtx2xy.len(), num_vtx * 2);
    let zero = [T::zero(), T::zero()];
    let mut area = T::zero();
    for i_edge in 0..num_vtx {
        let i0 = i_edge;
        let i1 = (i_edge + 1) % num_vtx;
        let p0 = &vtx2xy[i0 * 2..i0 * 2 + 2];
        let p1 = &vtx2xy[i1 * 2..i1 * 2 + 2];
        area += del_geo::tri2::area_(&zero, p0, p1);
    }
    area
}

pub fn from_circle(
    rad: f32,
    n: usize) -> nalgebra::Matrix2xX::<f32> {
    let mut vtx2xy = nalgebra::Matrix2xX::<f32>::zeros(n);
    for i in 0..n {
        let theta = std::f32::consts::PI * 2_f32 * i as f32 / n as f32;
        vtx2xy.column_mut(i).x = rad * f32::cos(theta);
        vtx2xy.column_mut(i).y = rad * f32::sin(theta);
    }
    vtx2xy
}

pub fn from_pentagram<Real>(
    center: &[Real],
    scale: Real) -> Vec<Real>
    where Real: num_traits::Float + 'static + Copy,
          f64: AsPrimitive<Real>,
          usize: AsPrimitive<Real>
{
    let dt: Real = (std::f64::consts::PI / 5_f64).as_();
    let hp: Real = (std::f64::consts::FRAC_PI_2).as_();
    let ratio = (2f64 / (3f64 + 5f64.sqrt())).as_();
    let mut xys = Vec::<Real>::new();
    for i in 0..10usize {
        let rad = if i % 2 == 0 { scale } else { ratio * scale };
        xys.push((dt * i.as_()+hp).cos() * rad + center[0]);
        xys.push((dt * i.as_()+hp).sin() * rad + center[1]);
    }
    xys
}

pub fn is_inside<Real>(
    vtx2xy: &[Real],
    p: &[Real]) -> bool
    where Real: num_traits::Float + Copy + 'static + std::ops::AddAssign,
          f64: AsPrimitive<Real>
{
    let num_vtx = vtx2xy.len() / 2;
    let mut wn: Real = Real::zero();
    for i in 0..num_vtx {
        let j = (i + 1) % num_vtx;
        wn += del_geo::edge2::winding_number_(
            &vtx2xy[i * 2..(i + 1) * 2],
            &vtx2xy[j * 2..(j + 1) * 2], p);
    }
    if (wn - Real::one()).abs() < 0.1.as_() { return true; }
    false
}

pub fn to_uniform_density_random_points<Real>(
    vtx2xy: &[Real],
    cell_len: Real,
    rng: &mut rand::rngs::StdRng) -> Vec<Real>
    where Real: num_traits::Float + std::ops::AddAssign + Copy + 'static + AsPrimitive<usize>,
          rand::distributions::Standard: rand::prelude::Distribution<Real>,
          f64: AsPrimitive<Real>,
          usize: AsPrimitive<Real>
{
    let aabb = del_geo::aabb2::from_vtx2xy(vtx2xy);
    use rand::Rng;
    let base_pos = [aabb[0] - cell_len * rng.gen::<Real>(), aabb[1] - cell_len * rng.gen::<Real>()];
    let nx = ((aabb[2] - base_pos[0]) / cell_len).as_() + 1;
    let ny = ((aabb[3] - base_pos[1]) / cell_len).as_() + 1;
    let mut res = vec!();
    for ix in 0..nx {
        for iy in 0..ny {
            let x = base_pos[0] + (ix.as_() + rng.gen::<Real>()) * cell_len;
            let y = base_pos[1] + (iy.as_() + rng.gen::<Real>()) * cell_len;
            let is_inside = is_inside(vtx2xy, &[x, y]);
            if !is_inside { continue; }
            res.push(x);
            res.push(y);
        }
    }
    res
}

// -----------------------------------------------------

#[test]
fn test_circle() {
    let vtx2xy0 = from_circle(1.0, 300);
    let arclen0 = crate::polyloop::arclength::<f32, 2>(vtx2xy0.as_slice());
    assert!((arclen0 - 2. * std::f32::consts::PI).abs() < 1.0e-3);
    //
    {
        let ndiv1 = 330;
        let vtx2xy1 = crate::polyloop::resample::<f32, 2>(vtx2xy0.as_slice(), ndiv1);
        assert_eq!(vtx2xy1.len(), ndiv1 * 2);
        let arclen1 = crate::polyloop::arclength::<f32, 2>(vtx2xy1.as_slice());
        assert!((arclen0 - arclen1).abs() < 1.0e-3);
        let edge2length1 = crate::polyloop::edge2length::<f32, 2>(vtx2xy1.as_slice());
        let min_edge_len1 = edge2length1.iter()
            .min_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
        assert!((min_edge_len1 - arclen1 / ndiv1 as f32).abs() < 1.0e-3);
    }
    {
        let ndiv2 = 156;
        let vtx2xy2 = crate::polyloop::resample::<f32, 2>(vtx2xy0.as_slice(), ndiv2);
        assert_eq!(vtx2xy2.len(), ndiv2 * 2);
        let arclen2 = crate::polyloop::arclength::<f32, 2>(vtx2xy2.as_slice());
        assert!((arclen0 - arclen2).abs() < 1.0e-3);
        let edge2length2 = crate::polyloop::edge2length::<f32, 2>(vtx2xy2.as_slice());
        let min_edge_len2 = edge2length2.iter()
            .min_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
        assert!((min_edge_len2 - arclen2 / ndiv2 as f32).abs() < 1.0e-3);
    }
}