Repository navigation
Expand file tree
/
Copy pathellipse.go
More file actions
372 lines (328 loc) · 10.8 KB
/
Copy pathellipse.go
File metadata and controls
372 lines (328 loc) · 10.8 KB
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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
// SPDX-FileCopyrightText: 2018 Raph Levien
// SPDX-FileCopyrightText: 2024 Dominik Honnef and contributors
//
// SPDX-License-Identifier: MIT
// SPDX-FileAttributionText: https://github.com/linebender/kurbo
package curve
import (
"math"
)
type Ellipse struct {
inner Affine
}
var _ ClosedShape = Ellipse{}
// NewEllipse returns a new ellipse with a given center, radii, and rotation.
//
// The returned ellipse will be the result of taking a circle, stretching
// it by the radii along the x and y axes, then rotating it from the
// x axis by xRotation radians, before finally translating the center
// to center.
//
// Rotation is clockwise in a y-down coordinate system. For more on
// rotation, see [Rotate].
func NewEllipse(center Point, radii Vec2, xRotation Angle) Ellipse {
rx, ry := radii.Splat()
return newEllipse(Vec2(center), rx, ry, xRotation)
}
// NewEllipseFromRect returns the largest ellipse that can be bounded by the
// provided rectangle.
//
// This uses the absolute width and height of the rectangle.
//
// This ellipse is always axis-aligned; to apply rotation you can call
// [Ellipse.WithRotation] on the result.
func NewEllipseFromRect(rect Rect) Ellipse {
center := Vec2(rect.Center())
width, height := rect.Size().Scale(1.0 / 2.0).Splat()
return newEllipse(center, width, height, 0.0)
}
// NewEllipseFromAffine returns an ellipse from an affine transformation of the unit
// circle.
func NewEllipseFromAffine(aff Affine) Ellipse {
return Ellipse{inner: aff}
}
func NewEllipseFromCircle(c Circle) Ellipse {
return NewEllipse(c.Center, Vec(c.Radius, c.Radius), 0)
}
// WithCenter returns a new ellipse centered on the provided point.
func (e Ellipse) WithCenter(center Point) Ellipse {
return Ellipse{inner: e.inner.WithTranslation(Vec2(center))}
}
// WithRadii returns a new ellipse, with the radii replaced by the argument.
func (e Ellipse) WithRadii(radii Vec2) Ellipse {
rotation := e.inner.svd1()
translation := e.inner.Translation()
return newEllipse(translation, radii.X, radii.Y, rotation)
}
// WithRotation returns a new ellipse, with the rotation replaced by the
// argument.
//
// The rotation is clockwise, for a y-down coordinate system. For more
// on rotation, See [Rotate].
func (e Ellipse) WithRotation(rotation Angle) Ellipse {
scale := e.inner.svd0()
translation := e.inner.Translation()
return newEllipse(translation, scale.X, scale.Y, rotation)
}
// Since the circle is symmetric about the x and y axes, using absolute values for the
// radii results in the same ellipse. For simplicity we make this change here.
return Ellipse{
inner: Translate(center).
Mul(Rotate(xRotation)).
Mul(Scale(math.Abs(scaleX), math.Abs(scaleY))),
}
}
// Contains implements ClosedShape.
func (e Ellipse) Contains(pt Point) bool {
return e.Winding(pt) != 0
}
func (e Ellipse) IsInf() bool {
return e.inner.IsInf()
}
func (e Ellipse) IsNaN() bool {
return e.inner.IsNaN()
}
// Area implements ClosedShape.
func (e Ellipse) Area() float64 {
// An ellipse is a unit circle transformed by an affine transformation. The
// unsigned area of a transformed region is area × abs(det(affine)). The
// area of the unit circle is π.
return math.Pi * math.Abs(e.inner.Determinant())
}
// BoundingBox implements Shape.
func (e Ellipse) BoundingBox() Rect {
// Compute a tight bounding box of the ellipse.
//
// See https://www.iquilezles.org/articles/ellipses/. We can get the two
// radius vectors by applying the affine map to the two impulses (1, 0) and (0, 1) which gives
// (a, b) and (c, d) if the affine map is
//
// a | c | e
// -----------
// b | d | f
//
// We can then use the method in the link with the translation to get the bounding box.
aff := e.inner.Coefficients()
a2 := aff[0] * aff[0]
b2 := aff[1] * aff[1]
c2 := aff[2] * aff[2]
d2 := aff[3] * aff[3]
cx := aff[4]
cy := aff[5]
rangeX := math.Sqrt(a2 + c2)
rangeY := math.Sqrt(b2 + d2)
return Rect{
X0: cx - rangeX,
Y0: cy - rangeY,
X1: cx + rangeX,
Y1: cy + rangeY,
}
}
// Path implements Shape.
func (e Ellipse) Path(tolerance float64, out BezPath) BezPath {
radii, xRotation := e.inner.svd()
return Arc{
Center: e.Center(),
Radii: radii,
StartAngle: 0.0,
SweepAngle: 2 * math.Pi,
XRotation: xRotation,
}.Path(tolerance, out)
}
// PathLength returns the approximated ellipse perimeter.
//
// This uses a numerical approximation. The absolute error between the calculated perimeter
// and the true perimeter is bounded by `accuracy` (modulo floating point rounding errors).
//
// For circular ellipses (equal horizontal and vertical radii), the calculated perimeter is
// exact.
func (e Ellipse) PathLength(accuracy float64) float64 {
radii := e.Radii()
if radii.IsInf() {
return math.Inf(1)
}
// Check for the trivial case where the ellipse has one of its radii
// equal to 0, i.e., where it describes a line, as the numerical method
// used breaks down with this extreme.
if radii.X == 0 || radii.Y == 0 {
return 4 * max(radii.X, radii.Y)
}
// Evaluate an approximation based on a truncated infinite series. If it
// returns a good enough value, we do not need to iterate.
if kummerEllipticPerimeterRange(radii) <= accuracy {
return kummerEllipticPerimeter(radii)
}
return agmEllipticPerimeter(accuracy, radii)
}
// kummerEllipticPerimeter calculates circumference C of an ellipse with radii
// (x, y) as the infinite series
//
// C = π (x+y) · ∑ binom(1/2, n)^2 * h^n from n = 0 to ∞
//
// with h = (x - y)^2 / (x + y)^2
// and binom(.,.) the binomial coefficient
//
// as described by Kummer ("Über die Hypergeometrische Reihe", 1837) and
// rediscovered by Linderholm and Segal
// ("An Overlooked Series for the Elliptic Perimeter", 1995).
//
// The series converges very quickly for ellipses with only moderate
// eccentricity (h not close to 1).
//
// The series is truncated to the sixth power, meaning a lower bound on the true
// value is returned. Adding the value of [kummerEllipticPerimeterRange] to the
// value returned by this function calculates an upper bound on the true value.
func kummerEllipticPerimeter(radii Vec2) float64 {
x, y := radii.Splat()
h := pow2((x - y) / (x + y))
h2 := h * h
h3 := h2 * h
h4 := h3 * h
h5 := h4 * h
h6 := h5 * h
lower := math.Pi +
h*(math.Pi/4) +
h2*(math.Pi/64) +
h3*(math.Pi/256) +
h4*(math.Pi*25/16384) +
h5*(math.Pi*49/65536) +
h6*(math.Pi*441/1048576)
return (x + y) * lower
}
// kummerEllipticPerimeterRange this calculates the error range of
// [kummerEllipticPerimeter]. That function returns a lower bound on the true
// value, and though we do not know what the remainder of the infinite series
// sums to, we can calculate an upper bound:
//
// ∑ binom(1/2, n)^2 for n = 0 to inf
//
// = 1 + (1 / 2‼)² + (1‼ / 4‼)² + (3‼ / 6‼)² + (5‼ / 8‼)² + …
// = 4 / π
// with ‼ the [double factorial]
//
// (equation 274 in "Summation of Series", L. B. W. Jolley, 1961).
//
// This means the remainder of the infinite series for C, assuming the series was truncated to the
// mᵗʰ term and h = 1, sums to
// 4 / π - ∑ binom(1/2, n)² for n = 0 to m-1
//
// As 0 ≤ h ≤ 1, this is an upper bound.
//
// [double factorial]: https://en.wikipedia.org/wiki/Double_factorial
func kummerEllipticPerimeterRange(radii Vec2) float64 {
x, y := radii.Splat()
h := pow2((x - y) / (x + y))
const binomSquaredRemainder = 4.0/math.Pi -
(1.0 +
1.0/4.0 +
1.0/64.0 +
1.0/256.0 +
25.0/16384.0 +
49.0/65536.0 +
441.0/1048576.0)
return math.Pi * binomSquaredRemainder * pow7(h) * (x + y)
}
// agmEllipticPerimeter calculates circumference C of an ellipse with radii
// (x, y) using the arithmetic-geometric mean, as described in equation 19.8.6 of
// https://web.archive.org/web/20240926233336/https://dlmf.nist.gov/19.8#i.
func agmEllipticPerimeter(accuracy float64, radii Vec2) float64 {
var x, y float64
if radii.X >= radii.Y {
x, y = radii.Splat()
} else {
y, x = radii.Splat()
}
accuracy = accuracy / (2 * math.Pi * radii.X)
sum := 1.0
a := 1.0
g := y / x
c := math.Sqrt(1 - pow2(g))
mul := 0.5
for {
c2 := pow2(c)
// term = 2^(n-1) c_n^2
term := mul * c2
sum -= term
// We have c_(n+1) ≤ 1/2 c_n
// (for a derivation, see e.g. section 2.1 of "Elliptic integrals, the
// arithmetic-geometric mean and the Brent-Salamin algorithm for π" by G.J.O. Jameson:
// https://web.archive.org/web/20241002140956/https://www.maths.lancs.ac.uk/jameson/ellagm.pdf)
//
// Therefore
// ∑ 2^(i-1) c_i^2 from i = 1 ≤ ∑ 2^(i-1) ((1/2)^i c_0)^2 from i = 1
// = ∑ 2^-(i+1) c_0^2 from i = 1
// = 1/2 c_0^2
//
// or, for arbitrary starting point i = n+1:
// ∑ 2^(i-1) c_i^2 from i = n+1 ≤ ∑ 2^(i-1) ((1/2)^(i-n) c_n)^2 from i = n+1
// = ∑ 2^(2n - i - 1) c_n^2 from i = n+1
// = 2^(2n) ∑ 2^(-(i+1)) c_n^2 from i = n+1
// = 2^(2n) 2^(-(n+1)) c_n^2
// = 2^(n-1) c_n^2
//
// Therefore, the remainder of the series sums to less than or equal to 2^(n-1) c_n^2,
// which is exactly the value of the nth term.
//
// Furthermore, a_m ≥ g_n, and g_n ≤ 1, for all m, n.
if term <= accuracy*g {
// sum currently overestimates the true value - subtract the upper
// bound of the remaining series. We will then underestimate the
// true value, but by no more than 'accuracy'.
sum -= term
break
}
mul *= 2
// This is equal to c_next = c^2 / (4 * a_next)
c = (a - g) / 2
aNext := (a + g) / 2
g = math.Sqrt(a * g)
a = aNext
}
return 2 * math.Pi * radii.X / a * sum
}
// Winding implements ClosedShape.
func (e Ellipse) Winding(pt Point) int {
// Strategy here is to apply the inverse map to the point and see if it is in the unit
// circle.
inv := e.inner.Invert()
if Vec2(pt.Transform(inv)).Hypot2() < 1.0 {
return 1
} else {
return 0
}
}
// Center returns the center of the ellipse.
func (e Ellipse) Center() Point {
return Point(e.inner.Translation())
}
// Radii returns the two radii of the ellipse.
//
// The first number is the horizontal radius and the second is the
// vertical radius, before rotation.
func (e Ellipse) Radii() Vec2 {
radii := e.inner.svd0()
return radii
}
// Rotation returns the ellipse's rotation, in radians.
func (e Ellipse) Rotation() Angle {
rot := e.inner.svd1()
return rot
}
// RadiiRotation returns the radii and the rotation of this ellipse.
//
// This is equivalent to, but more efficiant than, using [Ellipse.Radii] and
// [Ellipse.Rotation].
func (e Ellipse) RadiiRotation() (Vec2, Angle) {
return e.inner.svd()
}
func (e Ellipse) Translate(v Vec2) Ellipse {
return Ellipse{
inner: Translate(v).Mul(e.inner),
}
}
func (e Ellipse) Transform(aff Affine) Ellipse {
return Ellipse{
inner: e.inner.Mul(aff),
}
}