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

use std::cmp::max;

use fj_math::{Circle, Point, Scalar};

use crate::objects::{Curve, CurveKind, GlobalCurve};

use super::{Approx, Local, Tolerance};

impl Approx for Curve {
    type Approximation = Vec<Local<Point<1>>>;

    fn approx(&self, tolerance: Tolerance) -> Self::Approximation {
        self.global().approx(tolerance)
    }
}

impl Approx for GlobalCurve {
    type Approximation = Vec<Local<Point<1>>>;

    /// Approximate the global curve
    ///
    /// # Implementation Note
    ///
    /// This only works as-is, because only circles need to be approximated
    /// right now and because only edges that are full circles are supported, as
    /// opposed to edges that only inhabit part of the circle.
    ///
    /// To support that, we will need additional information here, to define
    /// between which points the curve needs to be approximated.
    fn approx(&self, tolerance: Tolerance) -> Self::Approximation {
        match self.kind() {
            CurveKind::Circle(curve) => {
                approx_circle(curve, [[Scalar::ZERO], [Scalar::TAU]], tolerance)
            }
            CurveKind::Line(_) => Vec::new(),
        }
    }
}

/// Approximate a circle
///
/// `tolerance` specifies how much the approximation is allowed to deviate
/// from the circle.
pub fn approx_circle(
    circle: &Circle<3>,
    between: [impl Into<Point<1>>; 2],
    tolerance: Tolerance,
) -> Vec<Local<Point<1>>> {
    let mut points = Vec::new();

    let radius = circle.a().magnitude();

    let [start, end] = between.map(Into::into);
    let range = (end - start).t;

    // To approximate the circle, we use a regular polygon for which
    // the circle is the circumscribed circle. The `tolerance`
    // parameter is the maximum allowed distance between the polygon
    // and the circle. This is the same as the difference between
    // the circumscribed circle and the incircle.

    let n = number_of_vertices_for_circle(tolerance, radius, range.abs());

    for i in 0..n {
        let angle =
            start.t + (Scalar::TAU / n as f64 * i as f64) * range.sign();
        let point = circle.point_from_circle_coords([angle]);
        points.push(Local::new([angle], point));
    }

    points
}

fn number_of_vertices_for_circle(
    tolerance: Tolerance,
    radius: Scalar,
    range: Scalar,
) -> u64 {
    let n = (range / (Scalar::ONE - (tolerance.inner() / radius)).acos() / 2.)
        .ceil()
        .into_u64();

    max(n, 3)
}

#[cfg(test)]
mod tests {
    use fj_math::Scalar;

    use crate::algorithms::approx::Tolerance;

    #[test]
    fn number_of_vertices_for_circle() {
        verify_result(50., 100., Scalar::TAU, 3);
        verify_result(50., 100., Scalar::PI, 3);
        verify_result(10., 100., Scalar::TAU, 7);
        verify_result(10., 100., Scalar::PI, 4);
        verify_result(1., 100., Scalar::TAU, 23);
        verify_result(1., 100., Scalar::PI, 12);

        fn verify_result(
            tolerance: impl Into<Tolerance>,
            radius: impl Into<Scalar>,
            range: impl Into<Scalar>,
            n: u64,
        ) {
            let tolerance = tolerance.into();
            let radius = radius.into();
            let range = range.into();

            assert_eq!(
                n,
                super::number_of_vertices_for_circle(tolerance, radius, range)
            );

            assert!(calculate_error(radius, range, n) <= tolerance.inner());
            if n > 3 {
                assert!(
                    calculate_error(radius, range, n - 1) >= tolerance.inner()
                );
            }
        }

        fn calculate_error(radius: Scalar, range: Scalar, n: u64) -> Scalar {
            radius - radius * (range / Scalar::from_u64(n) / 2.).cos()
        }
    }
}