Source code

Revision control

Other Tools

1
/* This Source Code Form is subject to the terms of the Mozilla Public
2
* License, v. 2.0. If a copy of the MPL was not distributed with this
3
* file, You can obtain one at http://mozilla.org/MPL/2.0/. */
4
5
use api::units::*;
6
use euclid::Size2D;
7
use std::f32::consts::FRAC_PI_2;
8
9
10
/// Number of steps to integrate arc length over.
11
const STEP_COUNT: usize = 20;
12
13
/// Represents an ellipse centred at a local space origin.
14
#[derive(Debug, Clone)]
15
pub struct Ellipse<U> {
16
pub radius: Size2D<f32, U>,
17
pub total_arc_length: f32,
18
}
19
20
impl<U> Ellipse<U> {
21
pub fn new(radius: Size2D<f32, U>) -> Ellipse<U> {
22
// Approximate the total length of the first quadrant of this ellipse.
23
let total_arc_length = get_simpson_length(FRAC_PI_2, radius.width, radius.height);
24
25
Ellipse {
26
radius,
27
total_arc_length,
28
}
29
}
30
31
/// Binary search to estimate the angle of an ellipse
32
/// for a given arc length. This only searches over the
33
/// first quadrant of an ellipse.
34
pub fn find_angle_for_arc_length(&self, arc_length: f32) -> f32 {
35
// Clamp arc length to [0, pi].
36
let arc_length = arc_length.max(0.0).min(self.total_arc_length);
37
38
let epsilon = 0.01;
39
let mut low = 0.0;
40
let mut high = FRAC_PI_2;
41
let mut theta = 0.0;
42
let mut new_low = 0.0;
43
let mut new_high = FRAC_PI_2;
44
45
while low <= high {
46
theta = 0.5 * (low + high);
47
let length = get_simpson_length(theta, self.radius.width, self.radius.height);
48
49
if (length - arc_length).abs() < epsilon {
50
break;
51
} else if length < arc_length {
52
new_low = theta;
53
} else {
54
new_high = theta;
55
}
56
57
// If we have stopped moving down the arc, the answer that we have is as good as
58
// it is going to get. We break to avoid going into an infinite loop.
59
if new_low == low && new_high == high {
60
break;
61
}
62
63
high = new_high;
64
low = new_low;
65
}
66
67
theta
68
}
69
70
/// Get a point and tangent on this ellipse from a given angle.
71
/// This only works for the first quadrant of the ellipse.
72
pub fn get_point_and_tangent(&self, theta: f32) -> (LayoutPoint, LayoutPoint) {
73
let (sin_theta, cos_theta) = theta.sin_cos();
74
let point = LayoutPoint::new(
75
self.radius.width * cos_theta,
76
self.radius.height * sin_theta,
77
);
78
let tangent = LayoutPoint::new(
79
-self.radius.width * sin_theta,
80
self.radius.height * cos_theta,
81
);
82
(point, tangent)
83
}
84
85
pub fn contains(&self, point: LayoutPoint) -> bool {
86
self.signed_distance(point.to_vector()) <= 0.0
87
}
88
89
/// Find the signed distance from this ellipse given a point.
91
fn signed_distance(&self, point: LayoutVector2D) -> f32 {
92
// This algorithm fails for circles, so we handle them here.
93
if self.radius.width == self.radius.height {
94
return point.length() - self.radius.width;
95
}
96
97
let mut p = LayoutVector2D::new(point.x.abs(), point.y.abs());
98
let mut ab = self.radius.to_vector();
99
if p.x > p.y {
100
p = p.yx();
101
ab = ab.yx();
102
}
103
104
let l = ab.y * ab.y - ab.x * ab.x;
105
106
let m = ab.x * p.x / l;
107
let n = ab.y * p.y / l;
108
let m2 = m * m;
109
let n2 = n * n;
110
111
let c = (m2 + n2 - 1.0) / 3.0;
112
let c3 = c * c * c;
113
114
let q = c3 + m2 * n2 * 2.0;
115
let d = c3 + m2 * n2;
116
let g = m + m * n2;
117
118
let co = if d < 0.0 {
119
let p = (q / c3).acos() / 3.0;
120
let s = p.cos();
121
let t = p.sin() * (3.0_f32).sqrt();
122
let rx = (-c * (s + t + 2.0) + m2).sqrt();
123
let ry = (-c * (s - t + 2.0) + m2).sqrt();
124
(ry + l.signum() * rx + g.abs() / (rx * ry) - m) / 2.0
125
} else {
126
let h = 2.0 * m * n * d.sqrt();
127
let s = (q + h).signum() * (q + h).abs().powf(1.0 / 3.0);
128
let u = (q - h).signum() * (q - h).abs().powf(1.0 / 3.0);
129
let rx = -s - u - c * 4.0 + 2.0 * m2;
130
let ry = (s - u) * (3.0_f32).sqrt();
131
let rm = (rx * rx + ry * ry).sqrt();
132
let p = ry / (rm - rx).sqrt();
133
(p + 2.0 * g / rm - m) / 2.0
134
};
135
136
let si = (1.0 - co * co).sqrt();
137
let r = LayoutVector2D::new(ab.x * co, ab.y * si);
138
(r - p).length() * (p.y - r.y).signum()
139
}
140
}
141
142
/// Use Simpsons rule to approximate the arc length of
143
/// part of an ellipse. Note that this only works over
144
/// the range of [0, pi/2].
145
// TODO(gw): This is a simplistic way to estimate the
146
// arc length of an ellipse segment. We can probably use
147
// a faster / more accurate method!
148
fn get_simpson_length(theta: f32, rx: f32, ry: f32) -> f32 {
149
let df = theta / STEP_COUNT as f32;
150
let mut sum = 0.0;
151
152
for i in 0 .. (STEP_COUNT + 1) {
153
let (sin_theta, cos_theta) = (i as f32 * df).sin_cos();
154
let a = rx * sin_theta;
155
let b = ry * cos_theta;
156
let y = (a * a + b * b).sqrt();
157
let q = if i == 0 || i == STEP_COUNT {
158
1.0
159
} else if i % 2 == 0 {
160
2.0
161
} else {
162
4.0
163
};
164
165
sum += q * y;
166
}
167
168
(df / 3.0) * sum
169
}
170
171
#[cfg(test)]
172
pub mod test {
173
use super::*;
174
175
#[test]
176
fn find_angle_for_arc_length_for_long_eclipse() {
177
// Ensure that finding the angle on giant ellipses produces and answer and
178
// doesn't send us into an infinite loop.
179
let ellipse = Ellipse::new(LayoutSize::new(57500.0, 25.0));
180
let _ = ellipse.find_angle_for_arc_length(55674.53);
181
assert!(true);
182
183
let ellipse = Ellipse::new(LayoutSize::new(25.0, 57500.0));
184
let _ = ellipse.find_angle_for_arc_length(55674.53);
185
assert!(true);
186
}
187
}