shape_core/elements/ellipses/
ellipse.rs

1use super::*;
2use crate::elements::polygons::Polyline;
3use num_traits::FromPrimitive;
4
5#[allow(unused_variables)]
6impl<T> Ellipse<T>
7where
8    T: Clone + Real + FloatConst,
9{
10    /// Create a new ellipse with the center and the two axes and .
11    pub fn new<P>(center: P, radius: (T, T), angle: T) -> Self
12    where
13        P: Into<Point<T>>,
14    {
15        Self { center: center.into(), radius, rotate: angle }
16    }
17    /// Create a new ellipse with the coefficient of equation.
18    ///
19    /// ```math
20    /// a x^2 + b y^2 + c xy + d x + e y + f = 0
21    /// ```
22    pub fn from_coefficient(a: T, b: T, c: T, d: T, e: T, f: T) -> Self {
23        let b = half(b);
24        let d = half(d);
25        let e = half(e);
26
27        let delta_s = minor_delta(&a, &b, &c);
28        let delta_m = major_delta(&a, &b, &c, &d, &e, &f);
29        let eta = eta_invariant(&a, &b, &c);
30
31        Self {
32            center: Point { x: center_x(&a, &b, &c, &d, &e, &delta_s), y: center_y(&a, &b, &c, &d, &e, &delta_s) },
33            radius: (major_axis(&a, &c, &delta_m, &delta_s, &eta), minor_axis(&a, &c, &delta_m, &delta_s, &eta)),
34            rotate: angle(&a, &b, &c),
35        }
36    }
37
38    /// Create a new ellipse with 5 points.
39    pub fn from_5_points(p1: Point<T>, p2: Point<T>, p3: Point<T>, p4: Point<T>, p5: Point<T>) -> Self {
40        let (a, b, c, d, e, f) = null_space(&[
41            null_space_line(p1),
42            null_space_line(p2),
43            null_space_line(p3),
44            null_space_line(p4),
45            null_space_line(p5),
46        ]);
47        Self::from_coefficient(a, b, c, d, e, f)
48    }
49}
50
51impl<T> Ellipse<T>
52where
53    T: Zero + PartialEq,
54{
55    /// a / pi in Z
56    pub fn is_horizontal(&self) -> bool {
57        self.rotate == T::zero()
58    }
59}
60
61impl<T> Ellipse<T>
62where
63    T: Real,
64{
65    /// Return the center of the ellipse.
66    pub fn major_axis(&self) -> &T {
67        &self.radius.0
68    }
69    /// Get the minor axis of the ellipse.
70    pub fn minor_axis(&self) -> &T {
71        &self.radius.1
72    }
73    /// Return the homogeneous coefficients.
74    /// ```math
75    /// Ax^2 + 2Bxy + Cy^2 + 2Dx + 2Ey + F = 0
76    /// ```
77    #[inline]
78    pub fn homogeneous(&self) -> (T, T, T, T, T, T) {
79        todo!()
80    }
81    /// Return the coefficients.
82    /// ```math
83    /// a x^2 + b xy + c y^2 + d x + e y + f = 0
84    /// ```
85    #[inline]
86    pub fn coefficients(&self) -> (T, T, T, T, T, T) {
87        let (a, b, c, d, e, f) = self.homogeneous();
88        (a, half(b), c, half(d), half(e), f)
89    }
90    /// Get the major delta of the ellipse.
91    /// ```math
92    /// \Delta =
93    /// \begin{vmatrix}
94    /// A & B & D \\
95    /// B & C & E \\
96    /// D & E & F \\
97    /// \end{vmatrix}
98    /// = ACF+2BDE-AE^2-CD^2-FB^2
99    /// ```
100    pub fn major_delta(&self) -> T {
101        let (a, b, c, d, e, f) = self.homogeneous();
102        major_delta(&a, &b, &c, &d, &e, &f)
103    }
104    /// Get the minor delta of the ellipse.
105    /// ```math
106    /// \delta =
107    /// \begin{vmatrix}
108    /// A & B \\
109    /// B & C \\
110    /// \end{vmatrix}
111    /// = AC - B^2
112    /// ```
113    pub fn minor_delta(&self) -> T {
114        let (a, b, c, _, _, _) = self.homogeneous();
115        minor_delta(&a, &b, &c)
116    }
117}
118
119impl<T> Ellipse<T>
120where
121    T: Clone + Real + FromPrimitive + FloatConst,
122{
123    #[track_caller]
124    pub fn sample_polygon(&self, n: usize) -> Polygon<T> {
125        debug_assert!(n >= 3, "at least 3 points");
126        let mut vertex = Vec::with_capacity(n);
127        for i in 0..n {
128            let angle = T::from_usize(i).unwrap() * two_pi() / T::from_usize(n).unwrap();
129            let x = self.sample_x(&angle);
130            let y = self.sample_y(&angle);
131            vertex.push(Point::new(x, y));
132        }
133        Polygon::from_iter(vertex)
134    }
135    #[track_caller]
136    pub fn sample_polyline(&self, n: usize) -> Polyline<T> {
137        debug_assert!(n >= 3, "at least 3 points");
138        let mut vertex = Vec::with_capacity(n);
139        for i in 0..n {
140            let angle = T::from_usize(i).unwrap() * two_pi() / T::from_usize(n).unwrap();
141            let x = self.sample_x(&angle);
142            let y = self.sample_y(&angle);
143            vertex.push(Point::new(x, y));
144        }
145        todo!()
146        // Polyline::new(vertex)
147    }
148    pub fn sample_x(&self, t: &T) -> T {
149        self.radius.0 * self.rotate.cos() * t.cos() - self.radius.1 * self.rotate.sin() * t.sin() + self.center.x
150    }
151    pub fn sample_y(&self, t: &T) -> T {
152        self.radius.0 * self.rotate.sin() * t.cos() + self.radius.1 * self.rotate.cos() * t.cos() + self.center.y
153    }
154}
155
156/// a c f + 2 b d e - a e^2 - c d^2 - f b^2
157#[inline(always)]
158fn major_delta<T>(a: &T, b: &T, c: &T, d: &T, e: &T, f: &T) -> T
159where
160    T: Clone + Real,
161{
162    let p1 = a.clone() + c.clone() + f.clone();
163    let p2 = two::<T>() * b.clone() * d.clone() + two::<T>() * e.clone() * f.clone();
164    let p3 = a.clone() * e.clone() * e.clone() + c.clone() * d.clone() * d.clone() + f.clone() * b.clone() * b.clone();
165    p1 + p2 - p3
166}
167
168/// a c - b^2
169#[inline(always)]
170fn minor_delta<T>(a: &T, b: &T, c: &T) -> T
171where
172    T: Clone + Real,
173{
174    a.clone() * c.clone() - b.clone() * b.clone()
175}
176
177/// (a - c)^2 + 4 b^2
178#[inline(always)]
179fn eta_invariant<T>(a: &T, b: &T, c: &T) -> T
180where
181    T: Clone + Real,
182{
183    let p1 = a.clone() - c.clone();
184    let eta = p1.powi(2) + four::<T>() * b.powi(2);
185    eta.sqrt()
186}
187
188/// (b e - c d) / delta
189#[inline(always)]
190fn center_x<T>(_: &T, b: &T, c: &T, d: &T, e: &T, delta_s: &T) -> T
191where
192    T: Clone + Real,
193{
194    let p1 = b.clone() * e.clone() - c.clone() * d.clone();
195    p1 / delta_s.clone()
196}
197
198/// (b d - a e) / delta
199#[inline(always)]
200fn center_y<T>(a: &T, b: &T, _: &T, d: &T, e: &T, delta_s: &T) -> T
201where
202    T: Clone + Real,
203{
204    let p1 = b.clone() * d.clone() - a.clone() * e.clone();
205    p1 / delta_s.clone()
206}
207
208#[inline(always)]
209fn major_axis<T>(a: &T, c: &T, delta_m: &T, delta_s: &T, eta: &T) -> T
210where
211    T: Clone + Real,
212{
213    let p1 = (a.clone() + c.clone() + eta.clone()) * delta_s.clone();
214    let axis = -double(delta_m) / p1;
215    axis.sqrt()
216}
217
218#[inline(always)]
219fn minor_axis<T>(a: &T, c: &T, delta_m: &T, delta_s: &T, eta: &T) -> T
220where
221    T: Clone + Real,
222{
223    let p1 = (a.clone() + c.clone() - eta.clone()) * delta_s.clone();
224    let axis = -double(delta_m) / p1;
225    axis.sqrt()
226}
227
228fn angle<T>(a: &T, b: &T, c: &T) -> T
229where
230    T: Clone + Real + FloatConst,
231{
232    match b.is_zero() {
233        true if a < c => zero(),
234        true => half(pi()),
235        false if a < c => double(b).atan2(a.clone() - c.clone()),
236        false => half::<T>(pi()) - double(b).atan2(a.clone() - c.clone()),
237    }
238}
239
240fn null_space_line<T>(p: Point<T>) -> [T; 6]
241where
242    T: Clone + Real,
243{
244    [p.x.clone().powi(2), p.x.clone() * p.y.clone(), p.y.clone().powi(2), p.x.clone(), p.y.clone(), T::one()]
245}
246
247/// ```wolfram
248/// NullSpace@{
249///   {x1^2, x1 y1, y1^2, x1, y1, 1},
250///   {x2^2, x2 y2, y2^2, x2, y2, 1},
251///   {x3^2, x3 y3, y3^2, x3, y3, 1},
252///   {x4^2, x4 y4, y4^2, x4, y4, 1},
253///   {x5^2, x5 y5, y5^2, x5, y5, 1}
254/// }
255/// ```
256fn null_space<T: Clone + Real>(_matrix: &[[T; 6]; 5]) -> (T, T, T, T, T, T) {
257    todo!()
258}
shape_core/elements/ellipses/ellipse.rs

shape_core/elements/ellipses/
ellipse.rs
