Skip to main content

kcl_lib/std/
utils.rs

1use std::f64::consts::PI;
2use std::f64::consts::TAU;
3
4use kcl_api::UnitLength;
5use kittycad_modeling_cmds::shared::Angle;
6
7use super::args::TyF64;
8use crate::execution::types::NumericType;
9use crate::execution::types::NumericTypeExt;
10use crate::util::MathExt;
11
12pub(crate) fn untype_point(p: [TyF64; 2]) -> ([f64; 2], NumericType) {
13    let (x, y, ty) = NumericType::combine_eq_coerce(p[0].clone(), p[1].clone(), None);
14    ([x, y], ty)
15}
16
17pub(crate) fn untype_array<const N: usize>(p: [TyF64; N]) -> ([f64; N], NumericType) {
18    let (vec, ty) = NumericType::combine_eq_array(&p);
19    (
20        vec.try_into()
21            .unwrap_or_else(|v: Vec<f64>| panic!("Expected a Vec of length {} but it was {}", N, v.len())),
22        ty,
23    )
24}
25
26pub(crate) fn point_to_mm(p: [TyF64; 2]) -> [f64; 2] {
27    [p[0].to_mm(), p[1].to_mm()]
28}
29
30pub(crate) fn untyped_point_to_mm(p: [f64; 2], units: UnitLength) -> [f64; 2] {
31    untyped_point_to_unit(p, units, UnitLength::Millimeters)
32}
33
34pub fn untyped_point_to_unit(point: [f64; 2], from_len_unit: UnitLength, to_len_unit: UnitLength) -> [f64; 2] {
35    [
36        crate::execution::types::adjust_length(from_len_unit, point[0], to_len_unit).0,
37        crate::execution::types::adjust_length(from_len_unit, point[1], to_len_unit).0,
38    ]
39}
40
41pub(crate) fn point_to_len_unit(p: [TyF64; 2], len: UnitLength) -> [f64; 2] {
42    [p[0].to_length_units(len), p[1].to_length_units(len)]
43}
44
45/// Precondition, `p` must be in `len` units (this function does no conversion).
46pub(crate) fn point_to_typed(p: [f64; 2], len: UnitLength) -> [TyF64; 2] {
47    [
48        TyF64::new(p[0], NumericType::length(len)),
49        TyF64::new(p[1], NumericType::length(len)),
50    ]
51}
52
53pub(crate) fn point_3d_to_mm(p: [TyF64; 3]) -> [f64; 3] {
54    [p[0].to_mm(), p[1].to_mm(), p[2].to_mm()]
55}
56
57/// Get the distance between two points.
58pub(crate) fn distance(a: Coords2d, b: Coords2d) -> f64 {
59    ((b[0] - a[0]).squared() + (b[1] - a[1]).squared()).sqrt()
60}
61
62/// Get the angle between these points
63pub(crate) fn between(a: Coords2d, b: Coords2d) -> Angle {
64    let x = b[0] - a[0];
65    let y = b[1] - a[1];
66    normalize(Angle::from_radians(libm::atan2(y, x)))
67}
68
69/// Normalize the angle
70pub(crate) fn normalize(angle: Angle) -> Angle {
71    let deg = angle.to_degrees();
72    let result = ((deg % 360.0) + 360.0) % 360.0;
73    Angle::from_degrees(if result > 180.0 { result - 360.0 } else { result })
74}
75
76/// Gives the ▲-angle between from and to angles (shortest path)
77///
78/// Sign of the returned angle denotes direction, positive means counterClockwise 🔄
79/// # Examples
80///
81/// ```
82/// use std::f64::consts::PI;
83///
84/// use kcl_lib::std::utils::Angle;
85///
86/// assert_eq!(
87///     Angle::delta(Angle::from_radians(PI / 8.0), Angle::from_radians(PI / 4.0)),
88///     Angle::from_radians(PI / 8.0)
89/// );
90/// ```
91pub(crate) fn delta(from_angle: Angle, to_angle: Angle) -> Angle {
92    let norm_from_angle = normalize_rad(from_angle.to_radians());
93    let norm_to_angle = normalize_rad(to_angle.to_radians());
94    let provisional = norm_to_angle - norm_from_angle;
95
96    if provisional > -PI && provisional <= PI {
97        return Angle::from_radians(provisional);
98    }
99    if provisional > PI {
100        return Angle::from_radians(provisional - TAU);
101    }
102    if provisional < -PI {
103        return Angle::from_radians(provisional + TAU);
104    }
105    Angle::default()
106}
107
108pub(crate) fn normalize_rad(angle: f64) -> f64 {
109    let draft = angle % (TAU);
110    if draft < 0.0 { draft + TAU } else { draft }
111}
112
113fn calculate_intersection_of_two_lines(line1: &[Coords2d; 2], line2_angle: f64, line2_point: Coords2d) -> Coords2d {
114    let line2_point_b = [
115        line2_point[0] + libm::cos(line2_angle.to_radians()) * 10.0,
116        line2_point[1] + libm::sin(line2_angle.to_radians()) * 10.0,
117    ];
118    intersect(line1[0], line1[1], line2_point, line2_point_b)
119}
120
121fn intersect(p1: Coords2d, p2: Coords2d, p3: Coords2d, p4: Coords2d) -> Coords2d {
122    let slope = |p1: Coords2d, p2: Coords2d| (p1[1] - p2[1]) / (p1[0] - p2[0]);
123    let constant = |p1: Coords2d, p2: Coords2d| p1[1] - slope(p1, p2) * p1[0];
124    let get_y = |for_x: f64, p1: Coords2d, p2: Coords2d| slope(p1, p2) * for_x + constant(p1, p2);
125
126    if p1[0] == p2[0] {
127        return [p1[0], get_y(p1[0], p3, p4)];
128    }
129    if p3[0] == p4[0] {
130        return [p3[0], get_y(p3[0], p1, p2)];
131    }
132
133    let x = (constant(p3, p4) - constant(p1, p2)) / (slope(p1, p2) - slope(p3, p4));
134    let y = get_y(x, p1, p2);
135    [x, y]
136}
137
138pub(crate) fn intersection_with_parallel_line(
139    line1: &[Coords2d; 2],
140    line1_offset: f64,
141    line2_angle: f64,
142    line2_point: Coords2d,
143) -> Coords2d {
144    calculate_intersection_of_two_lines(&offset_line(line1_offset, line1[0], line1[1]), line2_angle, line2_point)
145}
146
147fn offset_line(offset: f64, p1: Coords2d, p2: Coords2d) -> [Coords2d; 2] {
148    if p1[0] == p2[0] {
149        let direction = (p1[1] - p2[1]).signum();
150        return [[p1[0] + offset * direction, p1[1]], [p2[0] + offset * direction, p2[1]]];
151    }
152    if p1[1] == p2[1] {
153        let direction = (p2[0] - p1[0]).signum();
154        return [[p1[0], p1[1] + offset * direction], [p2[0], p2[1] + offset * direction]];
155    }
156    let x_offset = offset / libm::sin(libm::atan2(p1[1] - p2[1], p1[0] - p2[0]));
157    [[p1[0] + x_offset, p1[1]], [p2[0] + x_offset, p2[1]]]
158}
159
160pub(crate) fn get_y_component(angle: Angle, x: f64) -> Coords2d {
161    let normalised_angle = ((angle.to_degrees() % 360.0) + 360.0) % 360.0; // between 0 and 360
162    let y = x * libm::tan(normalised_angle.to_radians());
163    let sign = if normalised_angle > 90.0 && normalised_angle <= 270.0 {
164        -1.0
165    } else {
166        1.0
167    };
168    [x * sign, y * sign]
169}
170
171pub(crate) fn get_x_component(angle: Angle, y: f64) -> Coords2d {
172    let normalised_angle = ((angle.to_degrees() % 360.0) + 360.0) % 360.0; // between 0 and 360
173    let x = y / libm::tan(normalised_angle.to_radians());
174    let sign = if normalised_angle > 180.0 && normalised_angle <= 360.0 {
175        -1.0
176    } else {
177        1.0
178    };
179    [x * sign, y * sign]
180}
181
182pub(crate) fn arc_center_and_end(
183    from: Coords2d,
184    start_angle: Angle,
185    end_angle: Angle,
186    radius: f64,
187) -> (Coords2d, Coords2d) {
188    let start_angle = start_angle.to_radians();
189    let end_angle = end_angle.to_radians();
190
191    let center = [
192        -(radius * libm::cos(start_angle) - from[0]),
193        -(radius * libm::sin(start_angle) - from[1]),
194    ];
195
196    let end = [
197        center[0] + radius * libm::cos(end_angle),
198        center[1] + radius * libm::sin(end_angle),
199    ];
200
201    (center, end)
202}
203
204// Calculate the center of 3 points using an algebraic method
205// Handles if 3 points lie on the same line (collinear) by returning the average of the points (could return None instead..)
206pub(crate) fn calculate_circle_center(p1: [f64; 2], p2: [f64; 2], p3: [f64; 2]) -> [f64; 2] {
207    let (x1, y1) = (p1[0], p1[1]);
208    let (x2, y2) = (p2[0], p2[1]);
209    let (x3, y3) = (p3[0], p3[1]);
210
211    // Compute the determinant d = 2 * (x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))
212    // Visually d is twice the area of the triangle formed by the points,
213    // also the same as: cross(p2 - p1, p3 - p1)
214    let d = 2.0 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2));
215
216    // If d is nearly zero, the points are collinear, and a unique circle cannot be defined.
217    if d.abs() < f64::EPSILON {
218        return [(x1 + x2 + x3) / 3.0, (y1 + y2 + y3) / 3.0];
219    }
220
221    // squared lengths
222    let p1_sq = x1 * x1 + y1 * y1;
223    let p2_sq = x2 * x2 + y2 * y2;
224    let p3_sq = x3 * x3 + y3 * y3;
225
226    // This formula is derived from the circle equations:
227    //   (x - cx)^2 + (y - cy)^2 = r^2
228    // All 3 points will satisfy this equation, so we have 3 equations. Radius can be eliminated
229    // by subtracting one of the equations from the other two and the remaining 2 equations can
230    // be solved for cx and cy.
231    [
232        (p1_sq * (y2 - y3) + p2_sq * (y3 - y1) + p3_sq * (y1 - y2)) / d,
233        (p1_sq * (x3 - x2) + p2_sq * (x1 - x3) + p3_sq * (x2 - x1)) / d,
234    ]
235}
236
237pub struct CircleParams {
238    pub center: Coords2d,
239    pub radius: f64,
240}
241
242pub fn calculate_circle_from_3_points(points: [Coords2d; 3]) -> CircleParams {
243    let center = calculate_circle_center(points[0], points[1], points[2]);
244    CircleParams {
245        center,
246        radius: distance(center, points[1]),
247    }
248}
249
250#[cfg(test)]
251mod tests {
252    // Here you can bring your functions into scope
253    use std::f64::consts::TAU;
254
255    use approx::assert_relative_eq;
256    use pretty_assertions::assert_eq;
257
258    use super::Angle;
259    use super::calculate_circle_center;
260    use super::get_x_component;
261    use super::get_y_component;
262    use crate::util::MathExt;
263
264    static EACH_QUAD: [(i32, [i32; 2]); 12] = [
265        (-315, [1, 1]),
266        (-225, [-1, 1]),
267        (-135, [-1, -1]),
268        (-45, [1, -1]),
269        (45, [1, 1]),
270        (135, [-1, 1]),
271        (225, [-1, -1]),
272        (315, [1, -1]),
273        (405, [1, 1]),
274        (495, [-1, 1]),
275        (585, [-1, -1]),
276        (675, [1, -1]),
277    ];
278
279    #[test]
280    fn test_get_y_component() {
281        let mut expected = Vec::new();
282        let mut results = Vec::new();
283
284        for &(angle, expected_result) in EACH_QUAD.iter() {
285            let res = get_y_component(Angle::from_degrees(angle as f64), 1.0);
286            results.push([res[0].round() as i32, res[1].round() as i32]);
287            expected.push(expected_result);
288        }
289
290        assert_eq!(results, expected);
291
292        let result = get_y_component(Angle::zero(), 1.0);
293        assert_eq!(result[0] as i32, 1);
294        assert_eq!(result[1] as i32, 0);
295
296        let result = get_y_component(Angle::from_degrees(90.0), 1.0);
297        assert_eq!(result[0] as i32, 1);
298        assert!(result[1] > 100000.0);
299
300        let result = get_y_component(Angle::from_degrees(180.0), 1.0);
301        assert_eq!(result[0] as i32, -1);
302        assert!((result[1] - 0.0).abs() < f64::EPSILON);
303
304        let result = get_y_component(Angle::from_degrees(270.0), 1.0);
305        assert_eq!(result[0] as i32, -1);
306        assert!(result[1] < -100000.0);
307    }
308
309    #[test]
310    fn test_get_x_component() {
311        let mut expected = Vec::new();
312        let mut results = Vec::new();
313
314        for &(angle, expected_result) in EACH_QUAD.iter() {
315            let res = get_x_component(Angle::from_degrees(angle as f64), 1.0);
316            results.push([res[0].round() as i32, res[1].round() as i32]);
317            expected.push(expected_result);
318        }
319
320        assert_eq!(results, expected);
321
322        let result = get_x_component(Angle::zero(), 1.0);
323        assert!(result[0] > 100000.0);
324        assert_eq!(result[1] as i32, 1);
325
326        let result = get_x_component(Angle::from_degrees(90.0), 1.0);
327        assert!((result[0] - 0.0).abs() < f64::EPSILON);
328        assert_eq!(result[1] as i32, 1);
329
330        let result = get_x_component(Angle::from_degrees(180.0), 1.0);
331        assert!(result[0] < -100000.0);
332        assert_eq!(result[1] as i32, 1);
333
334        let result = get_x_component(Angle::from_degrees(270.0), 1.0);
335        assert!((result[0] - 0.0).abs() < f64::EPSILON);
336        assert_eq!(result[1] as i32, -1);
337    }
338
339    #[test]
340    fn test_arc_center_and_end() {
341        let (center, end) = super::arc_center_and_end([0.0, 0.0], Angle::zero(), Angle::from_degrees(90.0), 1.0);
342        assert_eq!(center[0].round(), -1.0);
343        assert_eq!(center[1], 0.0);
344        assert_eq!(end[0].round(), -1.0);
345        assert_eq!(end[1], 1.0);
346
347        let (center, end) = super::arc_center_and_end([0.0, 0.0], Angle::zero(), Angle::from_degrees(180.0), 1.0);
348        assert_eq!(center[0].round(), -1.0);
349        assert_eq!(center[1], 0.0);
350        assert_eq!(end[0].round(), -2.0);
351        assert_eq!(end[1].round(), 0.0);
352
353        let (center, end) = super::arc_center_and_end([0.0, 0.0], Angle::zero(), Angle::from_degrees(180.0), 10.0);
354        assert_eq!(center[0].round(), -10.0);
355        assert_eq!(center[1], 0.0);
356        assert_eq!(end[0].round(), -20.0);
357        assert_eq!(end[1].round(), 0.0);
358    }
359
360    #[test]
361    fn test_calculate_circle_center() {
362        const EPS: f64 = 1e-4;
363
364        // Test: circle center = (4.1, 1.9)
365        let p1 = [1.0, 2.0];
366        let p2 = [4.0, 5.0];
367        let p3 = [7.0, 3.0];
368        let center = calculate_circle_center(p1, p2, p3);
369        assert_relative_eq!(center[0], 4.1, epsilon = EPS);
370        assert_relative_eq!(center[1], 1.9, epsilon = EPS);
371
372        // Tests: Generate a few circles and test its points
373        let center = [3.2, 0.7];
374        let radius_array = [0.001, 0.01, 0.6, 1.0, 5.0, 60.0, 500.0, 2000.0, 400_000.0];
375        let points_array = [[0.0, 0.33, 0.66], [0.0, 0.1, 0.2], [0.0, -0.1, 0.1], [0.0, 0.5, 0.7]];
376
377        let get_point = |radius: f64, t: f64| {
378            let angle = t * TAU;
379            [
380                center[0] + radius * libm::cos(angle),
381                center[1] + radius * libm::sin(angle),
382            ]
383        };
384
385        for radius in radius_array {
386            for point in points_array {
387                let p1 = get_point(radius, point[0]);
388                let p2 = get_point(radius, point[1]);
389                let p3 = get_point(radius, point[2]);
390                let c = calculate_circle_center(p1, p2, p3);
391                assert_relative_eq!(c[0], center[0], epsilon = EPS);
392                assert_relative_eq!(c[1], center[1], epsilon = EPS);
393            }
394        }
395
396        // Test: Equilateral triangle
397        let p1 = [0.0, 0.0];
398        let p2 = [1.0, 0.0];
399        let p3 = [0.5, 3.0_f64.sqrt() / 2.0];
400        let center = calculate_circle_center(p1, p2, p3);
401        assert_relative_eq!(center[0], 0.5, epsilon = EPS);
402        assert_relative_eq!(center[1], 1.0 / (2.0 * 3.0_f64.sqrt()), epsilon = EPS);
403
404        // Test: Collinear points (should return the average of the points)
405        let p1 = [0.0, 0.0];
406        let p2 = [1.0, 0.0];
407        let p3 = [2.0, 0.0];
408        let center = calculate_circle_center(p1, p2, p3);
409        assert_relative_eq!(center[0], 1.0, epsilon = EPS);
410        assert_relative_eq!(center[1], 0.0, epsilon = EPS);
411
412        // Test: Points forming a circle with radius = 1
413        let p1 = [0.0, 0.0];
414        let p2 = [0.0, 2.0];
415        let p3 = [2.0, 0.0];
416        let center = calculate_circle_center(p1, p2, p3);
417        assert_relative_eq!(center[0], 1.0, epsilon = EPS);
418        assert_relative_eq!(center[1], 1.0, epsilon = EPS);
419
420        // Test: Integer coordinates
421        let p1 = [0.0, 0.0];
422        let p2 = [0.0, 6.0];
423        let p3 = [6.0, 0.0];
424        let center = calculate_circle_center(p1, p2, p3);
425        assert_relative_eq!(center[0], 3.0, epsilon = EPS);
426        assert_relative_eq!(center[1], 3.0, epsilon = EPS);
427        // Verify radius (should be 3 * sqrt(2))
428        let radius = ((center[0] - p1[0]).squared() + (center[1] - p1[1]).squared()).sqrt();
429        assert_relative_eq!(radius, 3.0 * 2.0_f64.sqrt(), epsilon = EPS);
430    }
431}
432
433pub(crate) type Coords2d = [f64; 2];
434
435pub fn is_points_ccw_wasm(points: &[f64]) -> i32 {
436    // CCW is positive as that the Math convention
437
438    let mut sum = 0.0;
439    for i in 0..(points.len() / 2) {
440        let point1 = [points[2 * i], points[2 * i + 1]];
441        let point2 = [points[(2 * i + 2) % points.len()], points[(2 * i + 3) % points.len()]];
442        sum += (point2[0] + point1[0]) * (point2[1] - point1[1]);
443    }
444    sum.signum() as i32
445}
446
447pub(crate) fn is_points_ccw(points: &[Coords2d]) -> i32 {
448    let flattened_points: Vec<f64> = points.iter().flat_map(|&p| vec![p[0], p[1]]).collect();
449    is_points_ccw_wasm(&flattened_points)
450}
451
452fn get_slope(start: Coords2d, end: Coords2d) -> (f64, f64) {
453    let slope = if start[0] - end[0] == 0.0 {
454        f64::INFINITY
455    } else {
456        (start[1] - end[1]) / (start[0] - end[0])
457    };
458
459    let perp_slope = if slope == f64::INFINITY { 0.0 } else { -1.0 / slope };
460
461    (slope, perp_slope)
462}
463
464fn get_angle(point1: Coords2d, point2: Coords2d) -> f64 {
465    let delta_x = point2[0] - point1[0];
466    let delta_y = point2[1] - point1[1];
467    let angle = libm::atan2(delta_y, delta_x);
468
469    let result = if angle < 0.0 { angle + TAU } else { angle };
470    result * (180.0 / PI)
471}
472
473fn delta_angle(from_angle: f64, to_angle: f64) -> f64 {
474    let norm_from_angle = normalize_rad(from_angle);
475    let norm_to_angle = normalize_rad(to_angle);
476    let provisional = norm_to_angle - norm_from_angle;
477
478    if provisional > -PI && provisional <= PI {
479        provisional
480    } else if provisional > PI {
481        provisional - TAU
482    } else if provisional < -PI {
483        provisional + TAU
484    } else {
485        provisional
486    }
487}
488
489fn deg2rad(deg: f64) -> f64 {
490    deg * (PI / 180.0)
491}
492
493fn get_mid_point(
494    center: Coords2d,
495    arc_start_point: Coords2d,
496    arc_end_point: Coords2d,
497    tan_previous_point: Coords2d,
498    radius: f64,
499    obtuse: bool,
500) -> Coords2d {
501    let angle_from_center_to_arc_start = get_angle(center, arc_start_point);
502    let angle_from_center_to_arc_end = get_angle(center, arc_end_point);
503    let delta_ang = delta_angle(
504        deg2rad(angle_from_center_to_arc_start),
505        deg2rad(angle_from_center_to_arc_end),
506    );
507    let delta_ang = delta_ang / 2.0 + deg2rad(angle_from_center_to_arc_start);
508    let shortest_arc_mid_point: Coords2d = [
509        libm::cos(delta_ang) * radius + center[0],
510        libm::sin(delta_ang) * radius + center[1],
511    ];
512    let opposite_delta = delta_ang + PI;
513    let longest_arc_mid_point: Coords2d = [
514        libm::cos(opposite_delta) * radius + center[0],
515        libm::sin(opposite_delta) * radius + center[1],
516    ];
517
518    let rotation_direction_original_points = is_points_ccw(&[tan_previous_point, arc_start_point, arc_end_point]);
519    let rotation_direction_points_on_arc = is_points_ccw(&[arc_start_point, shortest_arc_mid_point, arc_end_point]);
520    if rotation_direction_original_points != rotation_direction_points_on_arc && obtuse {
521        longest_arc_mid_point
522    } else {
523        shortest_arc_mid_point
524    }
525}
526
527fn intersect_point_n_slope(point1: Coords2d, slope1: f64, point2: Coords2d, slope2: f64) -> Coords2d {
528    let x = if slope1.abs() == f64::INFINITY {
529        point1[0]
530    } else if slope2.abs() == f64::INFINITY {
531        point2[0]
532    } else {
533        (point2[1] - slope2 * point2[0] - point1[1] + slope1 * point1[0]) / (slope1 - slope2)
534    };
535    let y = if slope1.abs() != f64::INFINITY {
536        slope1 * x - slope1 * point1[0] + point1[1]
537    } else {
538        slope2 * x - slope2 * point2[0] + point2[1]
539    };
540    [x, y]
541}
542
543/// Structure to hold input data for calculating tangential arc information.
544pub struct TangentialArcInfoInput {
545    /// The starting point of the arc.
546    pub arc_start_point: Coords2d,
547    /// The ending point of the arc.
548    pub arc_end_point: Coords2d,
549    /// The point from which the tangent is drawn.
550    pub tan_previous_point: Coords2d,
551    /// Flag to determine if the arc is obtuse. Obtuse means it flows smoothly from the previous segment.
552    pub obtuse: bool,
553}
554
555/// Structure to hold the output data from calculating tangential arc information.
556pub struct TangentialArcInfoOutput {
557    /// The center point of the arc.
558    pub center: Coords2d,
559    /// The midpoint on the arc.
560    pub arc_mid_point: Coords2d,
561    /// The radius of the arc.
562    pub radius: f64,
563    /// Start angle of the arc in radians.
564    pub start_angle: f64,
565    /// End angle of the arc in radians.
566    pub end_angle: f64,
567    /// If the arc is counter-clockwise.
568    pub ccw: i32,
569    /// The length of the arc.
570    pub arc_length: f64,
571}
572
573// tanPreviousPoint and arcStartPoint make up a straight segment leading into the arc (of which the arc should be tangential). The arc should start at arcStartPoint and end at, arcEndPoint
574// With this information we should everything we need to calculate the arc's center and radius. However there is two tangential arcs possible, that just varies on their direction
575// One is obtuse where the arc smoothly flows from the straight segment, and the other would be acute that immediately cuts back in the other direction. The obtuse boolean is there to control for this.
576pub fn get_tangential_arc_to_info(input: TangentialArcInfoInput) -> TangentialArcInfoOutput {
577    let (_, perp_slope) = get_slope(input.tan_previous_point, input.arc_start_point);
578    let tangential_line_perp_slope = perp_slope;
579
580    // Calculate the midpoint of the line segment between arcStartPoint and arcEndPoint
581    let mid_point: Coords2d = [
582        (input.arc_start_point[0] + input.arc_end_point[0]) / 2.0,
583        (input.arc_start_point[1] + input.arc_end_point[1]) / 2.0,
584    ];
585
586    let slope_mid_point_line = get_slope(input.arc_start_point, mid_point);
587
588    let center: Coords2d;
589    let radius: f64;
590
591    if tangential_line_perp_slope == slope_mid_point_line.0 {
592        // can't find the intersection of the two lines if they have the same gradient
593        // but in this case the center is the midpoint anyway
594        center = mid_point;
595        radius = ((input.arc_start_point[0] - center[0]).squared() + (input.arc_start_point[1] - center[1]).squared())
596            .sqrt();
597    } else {
598        center = intersect_point_n_slope(
599            mid_point,
600            slope_mid_point_line.1,
601            input.arc_start_point,
602            tangential_line_perp_slope,
603        );
604        radius = ((input.arc_start_point[0] - center[0]).squared() + (input.arc_start_point[1] - center[1]).squared())
605            .sqrt();
606    }
607
608    let arc_mid_point = get_mid_point(
609        center,
610        input.arc_start_point,
611        input.arc_end_point,
612        input.tan_previous_point,
613        radius,
614        input.obtuse,
615    );
616
617    let start_angle = libm::atan2(
618        input.arc_start_point[1] - center[1],
619        input.arc_start_point[0] - center[0],
620    );
621    let end_angle = libm::atan2(input.arc_end_point[1] - center[1], input.arc_end_point[0] - center[0]);
622    let ccw = is_points_ccw(&[input.arc_start_point, arc_mid_point, input.arc_end_point]);
623
624    let arc_mid_angle = libm::atan2(arc_mid_point[1] - center[1], arc_mid_point[0] - center[0]);
625    let start_to_mid_arc_length = radius
626        * delta(Angle::from_radians(start_angle), Angle::from_radians(arc_mid_angle))
627            .to_radians()
628            .abs();
629    let mid_to_end_arc_length = radius
630        * delta(Angle::from_radians(arc_mid_angle), Angle::from_radians(end_angle))
631            .to_radians()
632            .abs();
633    let arc_length = start_to_mid_arc_length + mid_to_end_arc_length;
634
635    TangentialArcInfoOutput {
636        center,
637        radius,
638        arc_mid_point,
639        start_angle,
640        end_angle,
641        ccw,
642        arc_length,
643    }
644}
645
646#[cfg(test)]
647mod get_tangential_arc_to_info_tests {
648    use approx::assert_relative_eq;
649
650    use super::*;
651
652    fn round_to_three_decimals(num: f64) -> f64 {
653        (num * 1000.0).round() / 1000.0
654    }
655
656    #[test]
657    fn test_basic_case() {
658        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
659            tan_previous_point: [0.0, -5.0],
660            arc_start_point: [0.0, 0.0],
661            arc_end_point: [4.0, 0.0],
662            obtuse: true,
663        });
664        assert_relative_eq!(result.center[0], 2.0);
665        assert_relative_eq!(result.center[1], 0.0);
666        assert_relative_eq!(result.arc_mid_point[0], 2.0);
667        assert_relative_eq!(result.arc_mid_point[1], 2.0);
668        assert_relative_eq!(result.radius, 2.0);
669        assert_relative_eq!(result.start_angle, PI);
670        assert_relative_eq!(result.end_angle, 0.0);
671        assert_eq!(result.ccw, -1);
672    }
673
674    #[test]
675    fn basic_case_with_arc_centered_at_0_0_and_the_tangential_line_being_45_degrees() {
676        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
677            tan_previous_point: [0.0, -4.0],
678            arc_start_point: [2.0, -2.0],
679            arc_end_point: [-2.0, 2.0],
680            obtuse: true,
681        });
682        assert_relative_eq!(result.center[0], 0.0);
683        assert_relative_eq!(result.center[1], 0.0);
684        assert_relative_eq!(round_to_three_decimals(result.arc_mid_point[0]), 2.0);
685        assert_relative_eq!(round_to_three_decimals(result.arc_mid_point[1]), 2.0);
686        assert_relative_eq!(result.radius, (2.0f64 * 2.0 + 2.0 * 2.0).sqrt());
687        assert_relative_eq!(result.start_angle, -PI / 4.0);
688        assert_relative_eq!(result.end_angle, 3.0 * PI / 4.0);
689        assert_eq!(result.ccw, 1);
690    }
691
692    #[test]
693    fn test_get_tangential_arc_to_info_moving_arc_end_point() {
694        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
695            tan_previous_point: [0.0, -4.0],
696            arc_start_point: [2.0, -2.0],
697            arc_end_point: [2.0, 2.0],
698            obtuse: true,
699        });
700        let expected_radius = (2.0f64 * 2.0 + 2.0 * 2.0).sqrt();
701        assert_relative_eq!(round_to_three_decimals(result.center[0]), 0.0);
702        assert_relative_eq!(result.center[1], 0.0);
703        assert_relative_eq!(result.arc_mid_point[0], expected_radius);
704        assert_relative_eq!(round_to_three_decimals(result.arc_mid_point[1]), -0.0);
705        assert_relative_eq!(result.radius, expected_radius);
706        assert_relative_eq!(result.start_angle, -PI / 4.0);
707        assert_relative_eq!(result.end_angle, PI / 4.0);
708        assert_eq!(result.ccw, 1);
709    }
710
711    #[test]
712    fn test_get_tangential_arc_to_info_moving_arc_end_point_again() {
713        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
714            tan_previous_point: [0.0, -4.0],
715            arc_start_point: [2.0, -2.0],
716            arc_end_point: [-2.0, -2.0],
717            obtuse: true,
718        });
719        let expected_radius = (2.0f64 * 2.0 + 2.0 * 2.0).sqrt();
720        assert_relative_eq!(result.center[0], 0.0);
721        assert_relative_eq!(result.center[1], 0.0);
722        assert_relative_eq!(result.radius, expected_radius);
723        assert_relative_eq!(round_to_three_decimals(result.arc_mid_point[0]), 0.0);
724        assert_relative_eq!(result.arc_mid_point[1], expected_radius);
725        assert_relative_eq!(result.start_angle, -PI / 4.0);
726        assert_relative_eq!(result.end_angle, -3.0 * PI / 4.0);
727        assert_eq!(result.ccw, 1);
728    }
729
730    #[test]
731    fn test_get_tangential_arc_to_info_acute_moving_arc_end_point() {
732        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
733            tan_previous_point: [0.0, -4.0],
734            arc_start_point: [2.0, -2.0],
735            arc_end_point: [-2.0, -2.0],
736            obtuse: false,
737        });
738        let expected_radius = (2.0f64 * 2.0 + 2.0 * 2.0).sqrt();
739        assert_relative_eq!(result.center[0], 0.0);
740        assert_relative_eq!(result.center[1], 0.0);
741        assert_relative_eq!(result.radius, expected_radius);
742        assert_relative_eq!(round_to_three_decimals(result.arc_mid_point[0]), -0.0);
743        assert_relative_eq!(result.arc_mid_point[1], -expected_radius);
744        assert_relative_eq!(result.start_angle, -PI / 4.0);
745        assert_relative_eq!(result.end_angle, -3.0 * PI / 4.0);
746        // would be cw if it was obtuse
747        assert_eq!(result.ccw, -1);
748    }
749
750    #[test]
751    fn test_get_tangential_arc_to_info_obtuse_with_wrap_around() {
752        let arc_end = libm::cos(std::f64::consts::PI / 4.0) * 2.0;
753        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
754            tan_previous_point: [2.0, -4.0],
755            arc_start_point: [2.0, 0.0],
756            arc_end_point: [0.0, -2.0],
757            obtuse: true,
758        });
759        assert_relative_eq!(result.center[0], -0.0);
760        assert_relative_eq!(result.center[1], 0.0);
761        assert_relative_eq!(result.radius, 2.0);
762        assert_relative_eq!(result.arc_mid_point[0], -arc_end);
763        assert_relative_eq!(result.arc_mid_point[1], arc_end);
764        assert_relative_eq!(result.start_angle, 0.0);
765        assert_relative_eq!(result.end_angle, -PI / 2.0);
766        assert_eq!(result.ccw, 1);
767    }
768
769    #[test]
770    fn test_arc_length_obtuse_cw() {
771        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
772            tan_previous_point: [-1.0, -1.0],
773            arc_start_point: [-1.0, 0.0],
774            arc_end_point: [0.0, -1.0],
775            obtuse: true,
776        });
777        let circumference = TAU * result.radius;
778        let expected_length = circumference * 3.0 / 4.0; // 3 quarters of a circle circle
779        assert_relative_eq!(result.arc_length, expected_length);
780    }
781
782    #[test]
783    fn test_arc_length_acute_cw() {
784        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
785            tan_previous_point: [-1.0, -1.0],
786            arc_start_point: [-1.0, 0.0],
787            arc_end_point: [0.0, 1.0],
788            obtuse: true,
789        });
790        let circumference = TAU * result.radius;
791        let expected_length = circumference / 4.0; // 1 quarters of a circle circle
792        assert_relative_eq!(result.arc_length, expected_length);
793    }
794
795    #[test]
796    fn test_arc_length_obtuse_ccw() {
797        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
798            tan_previous_point: [1.0, -1.0],
799            arc_start_point: [1.0, 0.0],
800            arc_end_point: [0.0, -1.0],
801            obtuse: true,
802        });
803        let circumference = TAU * result.radius;
804        let expected_length = circumference * 3.0 / 4.0; // 1 quarters of a circle circle
805        assert_relative_eq!(result.arc_length, expected_length);
806    }
807
808    #[test]
809    fn test_arc_length_acute_ccw() {
810        let result = get_tangential_arc_to_info(TangentialArcInfoInput {
811            tan_previous_point: [1.0, -1.0],
812            arc_start_point: [1.0, 0.0],
813            arc_end_point: [0.0, 1.0],
814            obtuse: true,
815        });
816        let circumference = TAU * result.radius;
817        let expected_length = circumference / 4.0; // 1 quarters of a circle circle
818        assert_relative_eq!(result.arc_length, expected_length);
819    }
820}
821
822pub(crate) fn get_tangent_point_from_previous_arc(
823    last_arc_center: Coords2d,
824    last_arc_ccw: bool,
825    last_arc_end: Coords2d,
826) -> Coords2d {
827    let angle_from_old_center_to_arc_start = get_angle(last_arc_center, last_arc_end);
828    let tangential_angle = angle_from_old_center_to_arc_start + if last_arc_ccw { -90.0 } else { 90.0 };
829    // What is the 10.0 constant doing???
830    [
831        libm::cos(tangential_angle.to_radians()) * 10.0 + last_arc_end[0],
832        libm::sin(tangential_angle.to_radians()) * 10.0 + last_arc_end[1],
833    ]
834}