diff --git a/crates/fj-math/src/arc.rs b/crates/fj-math/src/arc.rs index d5039731d..0cd991056 100644 --- a/crates/fj-math/src/arc.rs +++ b/crates/fj-math/src/arc.rs @@ -1,23 +1,25 @@ -use crate::{Point, Scalar}; +use num_traits::Float; + +use crate::{Point, Scalar, Vector}; /// Calculated geometry that is useful when dealing with an arc pub struct Arc { - /// Start point of the arc - pub start: Point<2>, - /// End point of the arc - pub end: Point<2>, /// Center of the circle the arc is constructed on pub center: Point<2>, + /// Radius of the circle the arc is constructed on pub radius: Scalar, + /// Angle of `start` relative to `center`, in radians /// /// Guaranteed to be less than `end_angle`. pub start_angle: Scalar, + /// Angle of `end` relative to `center`, in radians /// /// Guaranteed to be greater than `end_angle`. pub end_angle: Scalar, + /// True if `start` and `end` were switched to ensure `end_angle` > `start_angle` pub flipped_construction: bool, } @@ -27,14 +29,22 @@ impl Arc { pub fn from_endpoints_and_angle( p0: impl Into>, p1: impl Into>, - angle: Scalar, + angle_rad: Scalar, ) -> Self { - use num_traits::Float; + let p0 = p0.into(); + let p1 = p1.into(); - let (p0, p1) = (p0.into(), p1.into()); + // This is an implementation of this solution: + // https://math.stackexchange.com/a/87374 - let flipped_construction = angle <= Scalar::ZERO; - let angle_rad = angle.abs(); + let distance_between_endpoints = (p1 - p0).magnitude(); + let radius = distance_between_endpoints + / (2. * (angle_rad.abs().into_f64() / 2.).sin()); + let distance_center_to_midpoint = + (radius.powi(2) - (distance_between_endpoints.powi(2) / 4.)).sqrt(); + + let flipped_construction = angle_rad <= Scalar::ZERO; + let angle_rad = angle_rad.abs(); let [p0, p1] = if flipped_construction { [p1, p0] @@ -47,27 +57,25 @@ impl Arc { } else { (Scalar::ONE, Scalar::ZERO) }; - let [[x0, y0], [x1, y1]] = [p0, p1].map(|p| p.coords.components); - // https://math.stackexchange.com/questions/27535/how-to-find-center-of-an-arc-given-start-point-end-point-radius-and-arc-direc - // distance between endpoints - let d = ((x1 - x0).powi(2) + (y1 - y0).powi(2)).sqrt(); - // radius - let r = d / (2. * (angle_rad.into_f64() / 2.).sin()); - // distance from center to midpoint between endpoints - let h = (r.powi(2) - (d.powi(2) / 4.)).sqrt(); - // (u, v) is the unit normal in the direction of p1 - p0 - let u = (x1 - x0) / d * uv_factor; - let v = (y1 - y0) / d * uv_factor; - // (cx, cy) is the center of the circle - let cx = ((x0 + x1) / 2.) - h * v; - let cy = ((y0 + y1) / 2.) + h * u; - let start_angle = (y0 - cy).atan2(x0 - cx); - let end_angle = (y1 - cy).atan2(x1 - cx) + end_angle_offset; + let unit_vector_p0_to_p1 = + (p1 - p0) / distance_between_endpoints * uv_factor; + let unit_vector_midpoint_to_center = + Vector::from([-unit_vector_p0_to_p1.v, unit_vector_p0_to_p1.u]); + let center = Point { + coords: (p0.coords + p1.coords) / 2. + + unit_vector_midpoint_to_center * distance_center_to_midpoint, + }; + let start_angle = { + let center_to_start = p0 - center; + center_to_start.v.atan2(center_to_start.u) + }; + let end_angle = { + let center_to_end = p1 - center; + center_to_end.v.atan2(center_to_end.u) + end_angle_offset + }; Self { - start: p0, - end: p1, - center: Point::from([cx, cy]), - radius: r, + center, + radius, start_angle, end_angle, flipped_construction, @@ -77,39 +85,11 @@ impl Arc { #[cfg(test)] mod tests { - use crate::{Point, Scalar}; + use crate::{Point, Scalar, Vector}; use super::Arc; - use approx::AbsDiffEq; - - fn check_arc_calculation(center: [f64; 2], radius: f64, a0: f64, a1: f64) { - let angle = a1 - a0; - - let p0 = [center[0] + radius * a0.cos(), center[1] + radius * a0.sin()]; - let p1 = [center[0] + radius * a1.cos(), center[1] + radius * a1.sin()]; - - let arc = Arc::from_endpoints_and_angle(p0, p1, Scalar::from(angle)); - - let epsilon = Scalar::default_epsilon() * 10.; - - dbg!(center, arc.center); - dbg!(arc.start_angle); - dbg!(arc.end_angle); - dbg!(arc.flipped_construction); - assert!(arc.center.abs_diff_eq(&Point::from(center), epsilon)); - assert!(arc.radius.abs_diff_eq(&Scalar::from(radius), epsilon)); - - if a0 < a1 { - assert!(!arc.flipped_construction); - assert!(arc.start_angle.abs_diff_eq(&Scalar::from(a0), epsilon)); - assert!(arc.end_angle.abs_diff_eq(&Scalar::from(a1), epsilon)); - } else { - assert!(arc.flipped_construction); - assert!(arc.end_angle.abs_diff_eq(&Scalar::from(a0), epsilon)); - assert!(arc.start_angle.abs_diff_eq(&Scalar::from(a1), epsilon)); - } - } + use approx::{assert_abs_diff_eq, AbsDiffEq}; #[test] fn arc_construction() { @@ -144,4 +124,57 @@ mod tests { 270_f64.to_radians(), ); } + + fn check_arc_calculation( + center: impl Into>, + radius: f64, + a0: f64, + a1: f64, + ) { + let center = center.into(); + let angle = a1 - a0; + + let p0 = center + Vector::from([a0.cos(), a0.sin()]) * radius; + let p1 = center + Vector::from([a1.cos(), a1.sin()]) * radius; + + let arc = Arc::from_endpoints_and_angle(p0, p1, Scalar::from(angle)); + + let epsilon = Scalar::default_epsilon() * 10.; + + dbg!(arc.start_angle); + dbg!(arc.end_angle); + dbg!(arc.flipped_construction); + assert_abs_diff_eq!(arc.center, center, epsilon = epsilon); + assert_abs_diff_eq!( + arc.radius, + Scalar::from(radius), + epsilon = epsilon + ); + + if a0 < a1 { + assert!(!arc.flipped_construction); + assert_abs_diff_eq!( + arc.start_angle, + Scalar::from(a0), + epsilon = epsilon + ); + assert_abs_diff_eq!( + arc.end_angle, + Scalar::from(a1), + epsilon = epsilon + ); + } else { + assert!(arc.flipped_construction); + assert_abs_diff_eq!( + arc.end_angle, + Scalar::from(a0), + epsilon = epsilon + ); + assert_abs_diff_eq!( + arc.start_angle, + Scalar::from(a1), + epsilon = epsilon + ); + } + } }