GCC Code Coverage Report

Directory: ./
File: src/2geom/ellipse.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 316 392 80.6%
Functions: 28 34 82.4%
Branches: 264 534 49.4%
Line Branch Exec Source
1 /** @file
2 * @brief Ellipse shape
3 *//*
4 * Authors:
5 * Marco Cecchetti <mrcekets at gmail.com>
6 * Krzysztof KosiƄski <tweenk.pl@gmail.com>
7 *
8 * Copyright 2008-2014 Authors
9 *
10 * This library is free software; you can redistribute it and/or
11 * modify it either under the terms of the GNU Lesser General Public
12 * License version 2.1 as published by the Free Software Foundation
13 * (the "LGPL") or, at your option, under the terms of the Mozilla
14 * Public License Version 1.1 (the "MPL"). If you do not alter this
15 * notice, a recipient may use your version of this file under either
16 * the MPL or the LGPL.
17 *
18 * You should have received a copy of the LGPL along with this library
19 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 * You should have received a copy of the MPL along with this library
22 * in the file COPYING-MPL-1.1
23 *
24 * The contents of this file are subject to the Mozilla Public License
25 * Version 1.1 (the "License"); you may not use this file except in
26 * compliance with the License. You may obtain a copy of the License at
27 * http://www.mozilla.org/MPL/
28 *
29 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
30 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
31 * the specific language governing rights and limitations.
32 */
33
34 #include <2geom/conicsec.h>
35 #include <2geom/ellipse.h>
36 #include <2geom/elliptical-arc.h>
37 #include <2geom/numeric/fitting-tool.h>
38 #include <2geom/numeric/fitting-model.h>
39
40 namespace Geom {
41
42 ✗ Ellipse::Ellipse(Geom::Circle const &c)
43 ✗ : _center(c.center())
44 ✗ , _rays(c.radius(), c.radius())
45 ✗ , _angle(0)
46 ✗ {}
47
48 1203 void Ellipse::setCoefficients(double A, double B, double C, double D, double E, double F)
49 {
50 1203 double den = 4*A*C - B*B;
51
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1203 times.
 1203 if (den == 0) {
52 ✗ THROW_RANGEERROR("den == 0, while computing ellipse centre");
53 }
54 1203 _center[X] = (B*E - 2*C*D) / den;
55 1203 _center[Y] = (B*D - 2*A*E) / den;
56
57 // evaluate the a coefficient of the ellipse equation in normal form
58 // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1
59 // where b = a*B , c = a*C, (cx,cy) == centre
60 1203 double num = A * sqr(_center[X])
61 1203 + B * _center[X] * _center[Y]
62 1203 + C * sqr(_center[Y])
63 1203 - F;
64
65
66 //evaluate ellipse rotation angle
67
1/2
✓ Branch 2 taken 1203 times.
✗ Branch 3 not taken.
 1203 _angle = std::atan2( -B, -(A - C) )/2;
68
69 // evaluate the length of the ellipse rays
70 1203 double sinrot, cosrot;
71
1/2
✓ Branch 1 taken 1203 times.
✗ Branch 2 not taken.
 1203 sincos(_angle, sinrot, cosrot);
72 1203 double cos2 = cosrot * cosrot;
73 1203 double sin2 = sinrot * sinrot;
74 1203 double cossin = cosrot * sinrot;
75
76 1203 den = A * cos2 + B * cossin + C * sin2;
77
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1203 times.
 1203 if (den == 0) {
78 ✗ THROW_RANGEERROR("den == 0, while computing 'rx' coefficient");
79 }
80 1203 double rx2 = num / den;
81
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1203 times.
 1203 if (rx2 < 0) {
82 ✗ THROW_RANGEERROR("rx2 < 0, while computing 'rx' coefficient");
83 }
84 1203 _rays[X] = std::sqrt(rx2);
85
86 1203 den = C * cos2 - B * cossin + A * sin2;
87
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1203 times.
 1203 if (den == 0) {
88 ✗ THROW_RANGEERROR("den == 0, while computing 'ry' coefficient");
89 }
90 1203 double ry2 = num / den;
91
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1203 times.
 1203 if (ry2 < 0) {
92 ✗ THROW_RANGEERROR("ry2 < 0, while computing 'rx' coefficient");
93 }
94 1203 _rays[Y] = std::sqrt(ry2);
95
96 // the solution is not unique so we choose always the ellipse
97 // with a rotation angle between 0 and PI/2
98
1/2
✓ Branch 1 taken 1203 times.
✗ Branch 2 not taken.
 1203 makeCanonical();
99 1203 }
100
101 ✗ Point Ellipse::initialPoint() const
102 {
103 ✗ Coord sinrot, cosrot;
104 ✗ sincos(_angle, sinrot, cosrot);
105 ✗ Point p(ray(X) * cosrot + center(X), ray(X) * sinrot + center(Y));
106 ✗ return p;
107 }
108
109
110 364217 Affine Ellipse::unitCircleTransform() const
111 {
112
3/6
✓ Branch 2 taken 364217 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 364217 times.
✗ Branch 6 not taken.
✓ Branch 12 taken 364217 times.
✗ Branch 13 not taken.
 364217 Affine ret = Scale(ray(X), ray(Y)) * Rotate(_angle);
113
1/2
✓ Branch 3 taken 364217 times.
✗ Branch 4 not taken.
 364217 ret.setTranslation(center());
114 364217 return ret;
115 }
116
117 42608 Affine Ellipse::inverseUnitCircleTransform() const
118 {
119
3/6
✓ Branch 1 taken 42608 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 42608 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 42608 times.
 42608 if (ray(X) == 0 || ray(Y) == 0) {
120 ✗ THROW_RANGEERROR("a degenerate ellipse doesn't have an inverse unit circle transform");
121 }
122
4/8
✓ Branch 7 taken 42608 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 42608 times.
✗ Branch 11 not taken.
✓ Branch 19 taken 42608 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 42608 times.
✗ Branch 23 not taken.
 42608 Affine ret = Translate(-center()) * Rotate(-_angle) * Scale(1/ray(X), 1/ray(Y));
123 42608 return ret;
124 }
125
126
127 40000 LineSegment Ellipse::axis(Dim2 d) const
128 {
129 40000 Point a(0, 0), b(0, 0);
130 40000 a[d] = -1;
131 40000 b[d] = 1;
132
1/2
✓ Branch 1 taken 40000 times.
✗ Branch 2 not taken.
 40000 LineSegment ls(a, b);
133
2/4
✓ Branch 2 taken 40000 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 40000 times.
✗ Branch 6 not taken.
 40000 ls.transform(unitCircleTransform());
134 80000 return ls;
135 ✗ }
136
137 ✗ LineSegment Ellipse::semiaxis(Dim2 d, int sign) const
138 {
139 ✗ Point a(0, 0), b(0, 0);
140 ✗ b[d] = sgn(sign);
141 ✗ LineSegment ls(a, b);
142 ✗ ls.transform(unitCircleTransform());
143 ✗ return ls;
144 ✗ }
145
146 7 Rect Ellipse::boundsExact() const
147 {
148
1/2
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
 7 auto const trans = unitCircleTransform();
149
150 7 auto proj_bounds = [&] (Dim2 d) {
151 // The dth coordinate function pulls back to trans[d] * x + trans[d + 2] * y + trans[d + 4]
152 // in the coordinate system where the ellipse is a unit circle. We compute its range of
153 // values on the unit circle.
154 14 auto const r = std::hypot(trans[d], trans[d + 2]);
155 14 auto const mid = trans[d + 4];
156
1/2
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
 14 return Interval(mid - r, mid + r);
157 7 };
158
159
3/6
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 7 times.
✗ Branch 10 not taken.
 21 return { proj_bounds(X), proj_bounds(Y) };
160 }
161
162 30042 Rect Ellipse::boundsFast() const
163 {
164 // Every ellipse is contained in the circle with the same center and radius
165 // equal to the larger of the two rays. We return the bounding square
166 // of this circle (this is really fast but only exact for circles).
167
2/2
✓ Branch 2 taken 12030 times.
✓ Branch 3 taken 18012 times.
 30042 auto const larger_ray = (ray(X) > ray(Y) ? ray(X) : ray(Y));
168
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 30042 times.
 30042 assert(larger_ray >= 0.0);
169 30042 auto const rr = Point(larger_ray, larger_ray);
170
1/2
✓ Branch 5 taken 30042 times.
✗ Branch 6 not taken.
 90126 return Rect(_center - rr, _center + rr);
171 }
172
173 1199 std::vector<double> Ellipse::coefficients() const
174 {
175
1/2
✓ Branch 2 taken 1199 times.
✗ Branch 3 not taken.
 3597 std::vector<double> c(6);
176
1/2
✓ Branch 7 taken 1199 times.
✗ Branch 8 not taken.
 1199 coefficients(c[0], c[1], c[2], c[3], c[4], c[5]);
177 1199 return c;
178 ✗ }
179
180 21286 void Ellipse::coefficients(Coord &A, Coord &B, Coord &C, Coord &D, Coord &E, Coord &F) const
181 {
182
3/6
✓ Branch 1 taken 21286 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 21286 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 21286 times.
 21286 if (ray(X) == 0 || ray(Y) == 0) {
183 ✗ THROW_RANGEERROR("a degenerate ellipse doesn't have an implicit form");
184 }
185
186 21286 double cosrot, sinrot;
187
1/2
✓ Branch 1 taken 21286 times.
✗ Branch 2 not taken.
 21286 sincos(_angle, sinrot, cosrot);
188 21286 double cos2 = cosrot * cosrot;
189 21286 double sin2 = sinrot * sinrot;
190 21286 double cossin = cosrot * sinrot;
191 21286 double invrx2 = 1 / (ray(X) * ray(X));
192 21286 double invry2 = 1 / (ray(Y) * ray(Y));
193
194 21286 A = invrx2 * cos2 + invry2 * sin2;
195 21286 B = 2 * (invrx2 - invry2) * cossin;
196 21286 C = invrx2 * sin2 + invry2 * cos2;
197 21286 D = -2 * A * center(X) - B * center(Y);
198 21286 E = -2 * C * center(Y) - B * center(X);
199 21286 F = A * center(X) * center(X)
200 21286 + B * center(X) * center(Y)
201 21286 + C * center(Y) * center(Y)
202 21286 - 1;
203 21286 }
204
205
206 ✗ void Ellipse::fit(std::vector<Point> const &points)
207 {
208 ✗ size_t sz = points.size();
209 ✗ if (sz < 5) {
210 ✗ THROW_RANGEERROR("fitting error: too few points passed");
211 }
212 ✗ NL::LFMEllipse model;
213 ✗ NL::least_squeares_fitter<NL::LFMEllipse> fitter(model, sz);
214
215 ✗ for (size_t i = 0; i < sz; ++i) {
216 ✗ fitter.append(points[i]);
217 }
218 ✗ fitter.update();
219
220 ✗ NL::Vector z(sz, 0.0);
221 ✗ model.instance(*this, fitter.result(z));
222 ✗ }
223
224
225 EllipticalArc *
226 7 Ellipse::arc(Point const &ip, Point const &inner, Point const &fp)
227 {
228 // This is resistant to degenerate ellipses:
229 // both flags evaluate to false in that case.
230
231 7 bool large_arc_flag = false;
232 7 bool sweep_flag = false;
233
234 // Determination of large arc flag:
235 // large_arc is false when the inner point is on the same side
236 // of the center---initial point line as the final point, AND
237 // is on the same side of the center---final point line as the
238 // initial point.
239 // Additionally, large_arc is always false when we have exactly
240 // 1/2 of an arc, i.e. the cross product of the center -> initial point
241 // and center -> final point vectors is zero.
242 // Negating the above leads to the condition for large_arc being true.
243 7 Point fv = fp - _center;
244 7 Point iv = ip - _center;
245 7 Point innerv = inner - _center;
246 7 double ifcp = cross(fv, iv);
247
248
9/10
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4 times.
✓ Branch 12 taken 6 times.
✓ Branch 13 taken 1 times.
✓ Branch 15 taken 2 times.
✓ Branch 16 taken 5 times.
 22 if (ifcp != 0 && (sgn(cross(fv, innerv)) != sgn(ifcp) ||
249
4/4
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 3 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 3 times.
 15 sgn(cross(iv, innerv)) != sgn(-ifcp)))
250 {
251 2 large_arc_flag = true;
252 }
253
254 // Determination of sweep flag:
255 // For clarity, let's assume that Y grows up. Then the cross product
256 // is positive for points on the left side of a vector and negative
257 // on the right side of a vector.
258 //
259 // cross(?, v) > 0
260 // o------------------->
261 // cross(?, v) < 0
262 //
263 // If the arc is small (large_arc_flag is false) and the final point
264 // is on the right side of the vector initial point -> center,
265 // we have to go in the direction of increasing angles
266 // (counter-clockwise) and the sweep flag is true.
267 // If the arc is large, the opposite is true, since we have to reach
268 // the final point going the long way - in the other direction.
269 // We can express this observation as:
270 // cross(_center - ip, fp - _center) < 0 xor large_arc flag
271 // This is equal to:
272 // cross(-iv, fv) < 0 xor large_arc flag
273 // But cross(-iv, fv) is equal to cross(fv, iv) due to antisymmetry
274 // of the cross product, so we end up with the condition below.
275
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4 times.
 7 if ((ifcp < 0) ^ large_arc_flag) {
276 3 sweep_flag = true;
277 }
278
279 7 EllipticalArc *ret_arc = new EllipticalArc(ip, ray(X), ray(Y), rotationAngle(),
280
3/8
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 7 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 7 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 14 large_arc_flag, sweep_flag, fp);
281 14 return ret_arc;
282 }
283
284 10000 Ellipse &Ellipse::operator*=(Rotate const &r)
285 {
286 10000 _angle += r.angle();
287 10000 _center *= r;
288 10000 return *this;
289 }
290
291 1203 Ellipse &Ellipse::operator*=(Affine const& m)
292 {
293
3/6
✓ Branch 3 taken 1203 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1203 times.
✗ Branch 7 not taken.
✓ Branch 13 taken 1203 times.
✗ Branch 14 not taken.
 1203 Affine a = Scale(ray(X), ray(Y)) * Rotate(_angle);
294
1/2
✓ Branch 2 taken 1203 times.
✗ Branch 3 not taken.
 1203 Affine mwot = m.withoutTranslation();
295
1/2
✓ Branch 2 taken 1203 times.
✗ Branch 3 not taken.
 1203 Affine am = a * mwot;
296
1/2
✓ Branch 2 taken 1203 times.
✗ Branch 3 not taken.
 1203 Point new_center = _center * m;
297
298
2/4
✓ Branch 1 taken 1203 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1203 times.
 1203 if (are_near(am.descrim(), 0)) {
299 double angle;
300 ✗ if (am[0] != 0) {
301 ✗ angle = std::atan2(am[2], am[0]);
302 ✗ } else if (am[1] != 0) {
303 ✗ angle = std::atan2(am[3], am[1]);
304 } else {
305 ✗ angle = M_PI/2;
306 }
307 ✗ Point v = Point::polar(angle) * am;
308 ✗ _center = new_center;
309 ✗ _rays[X] = L2(v);
310 ✗ _rays[Y] = 0;
311 ✗ _angle = atan2(v);
312 ✗ return *this;
313
7/8
✓ Branch 1 taken 1203 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 14 times.
✓ Branch 4 taken 1189 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 4 times.
✓ Branch 9 taken 1199 times.
 1203 } else if (mwot.isScale(0) && _angle.radians() == 0) {
314 4 _rays[X] *= std::abs(mwot[0]);
315 4 _rays[Y] *= std::abs(mwot[3]);
316 4 _center = new_center;
317 4 return *this;
318 }
319
320
1/2
✓ Branch 2 taken 1199 times.
✗ Branch 3 not taken.
 1199 std::vector<double> coeff = coefficients();
321 1199 Affine q( coeff[0], coeff[1]/2,
322 1199 coeff[1]/2, coeff[2],
323 1199 0, 0 );
324
325
1/2
✓ Branch 2 taken 1199 times.
✗ Branch 3 not taken.
 1199 Affine invm = mwot.inverse();
326
1/2
✓ Branch 1 taken 1199 times.
✗ Branch 2 not taken.
 1199 q = invm * q ;
327 1199 std::swap(invm[1], invm[2]);
328
1/2
✓ Branch 1 taken 1199 times.
✗ Branch 2 not taken.
 1199 q *= invm;
329
1/2
✓ Branch 4 taken 1199 times.
✗ Branch 5 not taken.
 1199 setCoefficients(q[0], 2*q[1], q[3], 0, 0, -1);
330 1199 _center = new_center;
331
332 1199 return *this;
333 1199 }
334
335 36 Ellipse Ellipse::canonicalForm() const
336 {
337 36 Ellipse result(*this);
338 36 result.makeCanonical();
339 36 return result;
340 }
341
342 1239 void Ellipse::makeCanonical()
343 {
344
2/2
✓ Branch 2 taken 32 times.
✓ Branch 3 taken 1207 times.
 1239 if (_rays[X] == _rays[Y]) {
345
1/2
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
 32 _angle = 0;
346 32 return;
347 }
348
349
2/2
✓ Branch 1 taken 601 times.
✓ Branch 2 taken 606 times.
 1207 if (_angle < 0) {
350 601 _angle += M_PI;
351 }
352
2/2
✓ Branch 1 taken 599 times.
✓ Branch 2 taken 608 times.
 1207 if (_angle >= M_PI/2) {
353 599 std::swap(_rays[X], _rays[Y]);
354 599 _angle -= M_PI/2;
355 }
356 }
357
358 281378 Point Ellipse::pointAt(Coord t) const
359 {
360
1/2
✓ Branch 2 taken 281378 times.
✗ Branch 3 not taken.
 281378 Point p = Point::polar(t);
361
2/4
✓ Branch 2 taken 281378 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 281378 times.
✗ Branch 6 not taken.
 281378 p *= unitCircleTransform();
362 281378 return p;
363 }
364
365 ✗ Coord Ellipse::valueAt(Coord t, Dim2 d) const
366 {
367 ✗ Coord sinrot, cosrot, cost, sint;
368 ✗ sincos(rotationAngle(), sinrot, cosrot);
369 ✗ sincos(t, sint, cost);
370
371 ✗ if ( d == X ) {
372 ✗ return ray(X) * cosrot * cost
373 ✗ - ray(Y) * sinrot * sint
374 ✗ + center(X);
375 } else {
376 ✗ return ray(X) * sinrot * cost
377 ✗ + ray(Y) * cosrot * sint
378 ✗ + center(Y);
379 }
380 }
381
382 2512 Coord Ellipse::timeAt(Point const &p) const
383 {
384 // degenerate ellipse is basically a reparametrized line segment
385
3/6
✓ Branch 1 taken 2512 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2512 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2512 times.
 2512 if (ray(X) == 0 || ray(Y) == 0) {
386 ✗ if (ray(X) != 0) {
387 ✗ return asin(Line(axis(X)).timeAt(p));
388 ✗ } else if (ray(Y) != 0) {
389 ✗ return acos(Line(axis(Y)).timeAt(p));
390 } else {
391 ✗ return 0;
392 }
393 }
394
1/2
✓ Branch 2 taken 2512 times.
✗ Branch 3 not taken.
 2512 Affine iuct = inverseUnitCircleTransform();
395
3/6
✓ Branch 3 taken 2512 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2512 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2512 times.
✗ Branch 10 not taken.
 2512 return Angle(atan2(p * iuct)).radians0(); // return a value in [0, 2pi)
396 }
397
398 8 Point Ellipse::unitTangentAt(Coord t) const
399 {
400
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 Point p = Point::polar(t + M_PI/2);
401
3/6
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
 8 p *= unitCircleTransform().withoutTranslation();
402
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 p.normalize();
403 8 return p;
404 }
405
406 19863 bool Ellipse::contains(Point const &p) const
407 {
408
2/4
✓ Branch 3 taken 19863 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 19863 times.
✗ Branch 7 not taken.
 19863 Point tp = p * inverseUnitCircleTransform();
409 39726 return tp.length() <= 1;
410 }
411
412 /** @brief Convert curve time on the major axis to the corresponding angle
413 * parameters on a degenerate ellipse collapsed onto that axis.
414 * @param t The curve time on the major axis of an ellipse.
415 * @param vertical If true, the major axis goes from angle -π/2 to +π/2;
416 * otherwise, the major axis connects angles π and 0.
417 * @return The two angles at which the collapsed ellipse passes through the
418 * major axis point corresponding to the given time \f$t \in [0, 1]\f$.
419 */
420 60000 static std::array<Coord, 2> axis_time_to_angles(Coord t, bool vertical)
421 {
422 60000 Coord const to_unit = std::clamp(2.0 * t - 1.0, -1.0, 1.0);
423
2/2
✓ Branch 0 taken 10000 times.
✓ Branch 1 taken 50000 times.
 60000 if (vertical) {
424 10000 double const arcsin = std::asin(to_unit);
425 10000 return {arcsin, M_PI - arcsin};
426 } else {
427 50000 double const arccos = std::acos(to_unit);
428 50000 return {arccos, -arccos};
429 }
430 }
431
432 /** @brief For each intersection of some shape with the major axis of an ellipse, produce one or two
433 * intersections of a degenerate ellipse (collapsed onto that axis) with the same shape.
434 *
435 * @param axis_intersections The intersections of some shape with the major axis.
436 * @param vertical Whether this is the vertical major axis (in the ellipse's natural coordinates).
437 * @return A vector with doubled intersections (corresponding to the two passages of the squashed
438 * ellipse through that point) and swapped order of the intersected shapes.
439 */
440 40000 static std::vector<ShapeIntersection> double_axis_intersections(std::vector<ShapeIntersection> &&axis_intersections,
441 bool vertical)
442 {
443
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 40000 times.
 40000 if (axis_intersections.empty()) {
444 ✗ return {};
445 }
446 40000 std::vector<ShapeIntersection> result;
447
1/2
✓ Branch 2 taken 40000 times.
✗ Branch 3 not taken.
 40000 result.reserve(2 * axis_intersections.size());
448
449
2/2
✓ Branch 7 taken 60000 times.
✓ Branch 8 taken 40000 times.
 100000 for (auto const &x : axis_intersections) {
450
3/4
✓ Branch 2 taken 60000 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 120000 times.
✓ Branch 6 taken 60000 times.
 180000 for (auto a : axis_time_to_angles(x.second, vertical)) {
451
1/2
✓ Branch 3 taken 120000 times.
✗ Branch 4 not taken.
 120000 result.emplace_back(a, x.first, x.point()); // Swap first <-> converted second.
452
2/4
✓ Branch 0 taken 120000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 120000 times.
✗ Branch 3 not taken.
 120000 if (x.second == 0.0 || x.second == 1.0) {
453 break; // Do not double up endpoint intersections.
454 }
455 }
456 }
457 40000 return result;
458 40000 }
459
460 40053 std::vector<ShapeIntersection> Ellipse::intersect(Line const &line) const
461 {
462 40053 std::vector<ShapeIntersection> result;
463
464
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 40053 times.
 40053 if (line.isDegenerate()) {
465 ✗ return result;
466 }
467
6/6
✓ Branch 1 taken 30053 times.
✓ Branch 2 taken 10000 times.
✓ Branch 4 taken 10000 times.
✓ Branch 5 taken 20053 times.
✓ Branch 6 taken 20000 times.
✓ Branch 7 taken 20053 times.
 40053 if (ray(X) == 0 || ray(Y) == 0) {
468
3/6
✓ Branch 4 taken 20000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20000 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 20000 times.
✗ Branch 11 not taken.
 20000 return double_axis_intersections(line.intersect(majorAxis()), ray(X) == 0);
469 }
470
471 // Ax^2 + Bxy + Cy^2 + Dx + Ey + F
472 20053 std::array<Coord, 6> coeffs;
473
1/2
✓ Branch 7 taken 20053 times.
✗ Branch 8 not taken.
 20053 coefficients(coeffs[0], coeffs[1], coeffs[2], coeffs[3], coeffs[4], coeffs[5]);
474 20053 rescale_homogenous(coeffs);
475 20053 auto [A, B, C, D, E, F] = coeffs;
476
1/2
✓ Branch 2 taken 20053 times.
✗ Branch 3 not taken.
 20053 Affine iuct = inverseUnitCircleTransform();
477
478 // generic case
479 20053 std::array<Coord, 3> line_coeffs;
480
1/2
✓ Branch 4 taken 20053 times.
✗ Branch 5 not taken.
 20053 line.coefficients(line_coeffs[0], line_coeffs[1], line_coeffs[2]);
481 20053 rescale_homogenous(line_coeffs);
482 20053 auto [a, b, c] = line_coeffs;
483
1/2
✓ Branch 2 taken 20053 times.
✗ Branch 3 not taken.
 20053 Point lv = line.vector();
484
485
2/2
✓ Branch 2 taken 19128 times.
✓ Branch 3 taken 925 times.
 20053 if (fabs(lv[X]) > fabs(lv[Y])) {
486 // y = -a/b x - c/b
487 19128 Coord q = -a/b;
488 19128 Coord r = -c/b;
489
490 // substitute that into the ellipse equation, making it quadratic in x
491 19128 Coord I = A + B*q + C*q*q; // x^2 terms
492 19128 Coord J = B*r + C*2*q*r + D + E*q; // x^1 terms
493 19128 Coord K = C*r*r + E*r + F; // x^0 terms
494
1/2
✓ Branch 2 taken 19128 times.
✗ Branch 3 not taken.
 19128 std::vector<Coord> xs = solve_quadratic(I, J, K);
495
496
2/2
✓ Branch 7 taken 38255 times.
✓ Branch 8 taken 19128 times.
 57383 for (double x : xs) {
497 38255 Point p(x, q*x + r);
498
4/8
✓ Branch 2 taken 38255 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 38255 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 38255 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 38255 times.
✗ Branch 14 not taken.
 38255 result.emplace_back(atan2(p * iuct), line.timeAt(p), p);
499 }
500 19128 } else {
501 925 Coord q = -b/a;
502 925 Coord r = -c/a;
503
504 925 Coord I = A*q*q + B*q + C;
505 925 Coord J = A*2*q*r + B*r + D*q + E;
506 925 Coord K = A*r*r + D*r + F;
507
1/2
✓ Branch 2 taken 925 times.
✗ Branch 3 not taken.
 925 std::vector<Coord> xs = solve_quadratic(I, J, K);
508
509
2/2
✓ Branch 7 taken 1812 times.
✓ Branch 8 taken 925 times.
 2737 for (double x : xs) {
510 1812 Point p(q*x + r, x);
511
4/8
✓ Branch 2 taken 1812 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 1812 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1812 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1812 times.
✗ Branch 14 not taken.
 1812 result.emplace_back(atan2(p * iuct), line.timeAt(p), p);
512 }
513 925 }
514 20053 return result;
515 40053 }
516
517 30010 std::vector<ShapeIntersection> Ellipse::intersect(LineSegment const &seg) const
518 {
519
4/8
✓ Branch 2 taken 30010 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 30010 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 30010 times.
✗ Branch 10 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 30010 times.
 30010 if (!boundsFast().intersects(seg.boundsFast())) {
520 ✗ return {};
521 }
522
523 // We simply reuse the procedure for lines and filter out
524 // results where the line time value is outside of the unit interval,
525 // but we apply a small tolerance to account for numerical errors.
526
1/2
✓ Branch 1 taken 30010 times.
✗ Branch 2 not taken.
 30010 double const param_prec = EPSILON / seg.length(0.0);
527 // TODO: accept a precision setting instead of always using EPSILON
528 // (requires an ABI break).
529
530
2/4
✓ Branch 3 taken 30010 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 30010 times.
✗ Branch 7 not taken.
 30010 auto xings = intersect(Line(seg));
531
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 30010 times.
 30010 if (xings.empty()) {
532 ✗ return xings;
533 }
534 30010 decltype(xings) result;
535
1/2
✓ Branch 2 taken 30010 times.
✗ Branch 3 not taken.
 30010 result.reserve(xings.size());
536
537
2/2
✓ Branch 7 taken 60015 times.
✓ Branch 8 taken 30010 times.
 90025 for (auto const &x : xings) {
538
3/4
✓ Branch 0 taken 60013 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 60013 times.
 60015 if (x.second < -param_prec || x.second > 1.0 + param_prec) {
539 2 continue;
540 }
541
1/2
✓ Branch 6 taken 60013 times.
✗ Branch 7 not taken.
 60013 result.emplace_back(x.first, std::clamp(x.second, 0.0, 1.0), x.point());
542 }
543 30010 return result;
544 30010 }
545
546 20022 std::vector<ShapeIntersection> Ellipse::intersect(Ellipse const &other) const
547 {
548 // Handle degenerate cases first.
549
5/6
✓ Branch 1 taken 20022 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20000 times.
✓ Branch 5 taken 22 times.
✓ Branch 6 taken 20000 times.
✓ Branch 7 taken 22 times.
 20022 if (ray(X) == 0 || ray(Y) == 0) { // Degenerate ellipse, collapsed to the major axis.
550
3/6
✓ Branch 4 taken 20000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20000 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 20000 times.
✗ Branch 11 not taken.
 20000 return double_axis_intersections(other.intersect(majorAxis()), ray(X) == 0);
551 }
552
3/4
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 13 times.
 22 if (*this == other) { // Two identical ellipses.
553
1/2
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
 9 THROW_INFINITELY_MANY_SOLUTIONS("The two ellipses are identical.");
554 }
555
4/8
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 13 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 13 times.
✗ Branch 10 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 13 times.
 13 if (!boundsFast().intersects(other.boundsFast())) {
556 ✗ return {};
557 }
558
559 // Find coefficients of the implicit equations of the two ellipses and rescale
560 // them (losslessly) for better numerical conditioning.
561 13 std::array<double, 6> coeffs;
562
1/2
✓ Branch 7 taken 13 times.
✗ Branch 8 not taken.
 13 coefficients(coeffs[0], coeffs[1], coeffs[2], coeffs[3], coeffs[4], coeffs[5]);
563 13 rescale_homogenous(coeffs);
564 13 auto [A, B, C, D, E, F] = coeffs;
565
566 13 std::array<double, 6> otheffs;
567
1/2
✓ Branch 7 taken 13 times.
✗ Branch 8 not taken.
 13 other.coefficients(otheffs[0], otheffs[1], otheffs[2], otheffs[3], otheffs[4], otheffs[5]);
568 13 rescale_homogenous(otheffs);
569 13 auto [a, b, c, d, e, f] = otheffs;
570
571 // Assume that Q(x, y) = 0 is the ellipse equation given by uppercase letters
572 // and R(x, y) = 0 is the equation given by lowercase ones.
573 // In other words, Q is the quadratic function describing this ellipse and
574 // R is the quadratic function for the other ellipse.
575 //
576 // A point (x, y) is common to both ellipses if and only if it solves the system
577 // { Q(x, y) = 0,
578 // { R(x, y) = 0.
579 //
580 // If ” is any real number, we can multiply the first equation by ” and add that
581 // to the first equation, obtaining the new system of equations:
582 // { Q(x, y) = 0,
583 // { ”Q(x, y) + R(x, y) = 0.
584 //
585 // The first equation still says that (x, y) is a point on this ellipse, but the
586 // second equation uses the new expression (”Q + R) instead of the original R.
587 //
588 // Why do we do this? The reason is that the set of functions {”Q + R : ” real}
589 // is a "real system of conics" and there's a theorem which guarantees that such a system
590 // always contains a "degenerate conic" [proof below].
591 // In practice, the degenerate conic will describe a line or a pair of lines, and intersecting
592 // a line with an ellipse is much easier than intersecting two ellipses directly.
593 //
594 // But in order to be able to do this, we must find a value of ” for which ”Q + R is degenerate.
595 // We can write the expression (”Q + R)(x, y) in the following way:
596 //
597 // | aa bb/2 dd/2 | |x|
598 // (”Q + R)(x, y) = [x y 1] | bb/2 cc ee/2 | |y|
599 // | dd/2 ee/2 ff | |1|
600 //
601 // where aa = ”A + a and so on. The determinant can be explicitly written out,
602 // giving an equation which is cubic in ” and can be solved analytically.
603 // The conic ”Q + R is degenerate if and only if this determinant is 0.
604 //
605 // Proof that there's always a degenerate conic: a cubic real polynomial always has a root,
606 // and if the polynomial in ” isn't cubic (coefficient of ”^3 is zero), then the starting
607 // conic is already degenerate.
608
609 Coord I, J, K, L; // Coefficients of ” in the expression for the determinant.
610 13 I = (-B*B*F + 4*A*C*F + D*E*B - A*E*E - C*D*D) / 4;
611 13 J = -((B*B - 4*A*C) * f + (2*B*F - D*E) * b + (2*A*E - D*B) * e +
612 13 (2*C*D - E*B) * d + (D*D - 4*A*F) * c + (E*E - 4*C*F) * a) / 4;
613 13 K = -((b*b - 4*a*c) * F + (2*b*f - d*e) * B + (2*a*e - d*b) * E +
614 13 (2*c*d - e*b) * D + (d*d - 4*a*f) * C + (e*e - 4*c*f) * A) / 4;
615 13 L = (-b*b*f + 4*a*c*f + d*e*b - a*e*e - c*d*d) / 4;
616
617
1/2
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
 13 std::vector<Coord> mus = solve_cubic(I, J, K, L);
618 13 Coord mu = infinity();
619
620 // Now that we have solved for ”, we need to check whether the conic
621 // determined by ”Q + R is reducible to a product of two lines. If it's not,
622 // it means that there are no intersections. If it is, the intersections of these
623 // lines with the original ellipses (if there are any) give the coordinates
624 // of intersections.
625
626 // Prefer middle root if there are three.
627 // Out of three possible pairs of lines that go through four points of intersection
628 // of two ellipses, this corresponds to cross-lines. These intersect the ellipses
629 // at less shallow angles than the other two options.
630
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 7 times.
 13 if (mus.size() == 3) {
631 6 std::swap(mus[1], mus[0]);
632 }
633 /// Discriminant within this radius of 0 will be considered zero.
634 static Coord const discriminant_precision = 1e-10;
635
636
1/2
✓ Branch 7 taken 17 times.
✗ Branch 8 not taken.
 17 for (Coord candidate_mu : mus) {
637 17 Coord const aa = std::fma(candidate_mu, A, a);
638 17 Coord const bb = std::fma(candidate_mu, B, b);
639 17 Coord const cc = std::fma(candidate_mu, C, c);
640 17 Coord const delta = sqr(bb) - 4*aa*cc;
641
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 13 times.
 17 if (delta < -discriminant_precision) {
642 4 continue;
643 }
644 13 mu = candidate_mu;
645 13 break;
646 }
647
648 // if no suitable mu was found, there are no intersections
649
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 13 times.
 13 if (mu == infinity()) {
650 ✗ return {};
651 }
652
653 // Create the degenerate conic and decompose it into lines.
654 13 std::array<double, 6> degen = {std::fma(mu, A, a), std::fma(mu, B, b), std::fma(mu, C, c),
655 13 std::fma(mu, D, d), std::fma(mu, E, e), std::fma(mu, F, f)};
656 13 rescale_homogenous(degen);
657
1/2
✓ Branch 5 taken 13 times.
✗ Branch 6 not taken.
 39 auto const lines = xAx(degen[0], degen[1], degen[2],
658
1/2
✓ Branch 5 taken 13 times.
✗ Branch 6 not taken.
 26 degen[3], degen[4], degen[5]).decompose_df(discriminant_precision);
659
660 // intersect with the obtained lines and report intersections
661 13 std::vector<ShapeIntersection> result;
662
2/2
✓ Branch 2 taken 26 times.
✓ Branch 3 taken 13 times.
 39 for (auto const &line : lines) {
663
2/2
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 21 times.
 26 if (line.isDegenerate()) {
664 5 continue;
665 }
666
1/2
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
 21 auto as = intersect(line);
667 // NOTE: If we only cared about the intersection points, we could simply
668 // intersect this ellipse with the lines and ignore the other ellipse.
669 // But we need the time coordinates on the other ellipse as well.
670
1/2
✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
 21 auto bs = other.intersect(line);
671
5/6
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 15 times.
✓ Branch 6 taken 6 times.
✓ Branch 7 taken 15 times.
 21 if (as.empty() || bs.empty()) {
672 6 continue;
673 }
674 // Due to numerical errors, a tangency may sometimes be found as 1 intersection
675 // on one ellipse and 2 intersections on the other. If this happens, we average
676 // the points of the two intersections.
677 15 auto const intersection_average = [](ShapeIntersection const &i,
678 ShapeIntersection const &j) -> ShapeIntersection
679 {
680
1/2
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 2 return ShapeIntersection(i.first, j.first, middle_point(i.point(), j.point()));
681 };
682 15 auto const synthesize_intersection = [&](ShapeIntersection const &i,
683 ShapeIntersection const &j) -> void
684 {
685
2/4
✓ Branch 6 taken 24 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 24 times.
✗ Branch 10 not taken.
 24 result.emplace_back(i.first, j.first, middle_point(i.point(), j.point()));
686 24 };
687
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 4 times.
 15 if (as.size() == 2) {
688
2/2
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 2 times.
 11 if (bs.size() == 2) {
689
1/2
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
 9 synthesize_intersection(as[0], bs[0]);
690
1/2
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
 9 synthesize_intersection(as[1], bs[1]);
691
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 } else if (bs.size() == 1) {
692
2/4
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 2 synthesize_intersection(intersection_average(as[0], as[1]), bs[0]);
693 }
694
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 } else if (as.size() == 1) {
695
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if (bs.size() == 2) {
696 ✗ synthesize_intersection(as[0], intersection_average(bs[0], bs[1]));
697
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 } else if (bs.size() == 1) {
698
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 synthesize_intersection(as[0], bs[0]);
699 }
700 }
701
4/4
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 6 times.
✓ Branch 4 taken 15 times.
✓ Branch 5 taken 6 times.
 27 }
702 13 return result;
703 13 }
704
705 4 std::vector<ShapeIntersection> Ellipse::intersect(D2<Bezier> const &b) const
706 {
707 4 Coord A, B, C, D, E, F;
708
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 coefficients(A, B, C, D, E, F);
709
710 // We plug the X and Y curves into the implicit equation and solve for t.
711
13/26
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
✓ Branch 19 taken 4 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 4 times.
✗ Branch 23 not taken.
✓ Branch 30 taken 4 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 4 times.
✗ Branch 34 not taken.
✓ Branch 40 taken 4 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 4 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 4 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 4 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 4 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 4 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 4 times.
✗ Branch 59 not taken.
 4 Bezier x = A*b[X]*b[X] + B*b[X]*b[Y] + C*b[Y]*b[Y] + D*b[X] + E*b[Y] + F;
712
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector<Coord> r = x.roots();
713
714 4 std::vector<ShapeIntersection> result;
715
2/2
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 4 times.
 14 for (double & i : r) {
716
1/2
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 10 Point p = b.valueAt(i);
717
2/4
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
 10 result.emplace_back(timeAt(p), i, p);
718 }
719 4 return result;
720 4 }
721
722 31 bool Ellipse::operator==(Ellipse const &other) const
723 {
724
2/2
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 18 times.
 31 if (_center != other._center) return false;
725
726
1/2
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
 18 Ellipse a = this->canonicalForm();
727
1/2
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
 18 Ellipse b = other.canonicalForm();
728
729
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
 18 if (a._rays != b._rays) return false;
730
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
 18 if (a._angle != b._angle) return false;
731
732 18 return true;
733 }
734
735
736 28 bool are_near(Ellipse const &a, Ellipse const &b, Coord precision)
737 {
738 // We want to know whether no point on ellipse a is further than precision
739 // from the corresponding point on ellipse b. To check this, we compute
740 // the four extreme points at the end of each ray for each ellipse
741 // and check whether they are sufficiently close.
742
743 // First, we need to correct the angles on the ellipses, so that they are
744 // no further than M_PI/4 apart. This can always be done by rotating
745 // and exchanging axes.
746 28 Ellipse ac = a, bc = b;
747
3/4
✓ Branch 6 taken 28 times.
✗ Branch 7 not taken.
✓ Branch 12 taken 12 times.
✓ Branch 13 taken 16 times.
 28 if (distance(ac.rotationAngle(), bc.rotationAngle()).radians0() >= M_PI/2) {
748
1/2
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
 12 ac.setRotationAngle(ac.rotationAngle() + M_PI);
749 }
750
4/6
✓ Branch 6 taken 28 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 28 times.
✗ Branch 10 not taken.
✓ Branch 14 taken 10 times.
✓ Branch 15 taken 18 times.
 28 if (distance(ac.rotationAngle(), bc.rotationAngle()) >= M_PI/4) {
751
2/4
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 10 times.
✗ Branch 12 not taken.
 10 Angle d1 = distance(ac.rotationAngle() + M_PI/2, bc.rotationAngle());
752
2/4
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 10 times.
✗ Branch 12 not taken.
 10 Angle d2 = distance(ac.rotationAngle() - M_PI/2, bc.rotationAngle());
753
2/2
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 6 times.
 10 Coord adj = d1.radians0() < d2.radians0() ? M_PI/2 : -M_PI/2;
754
1/2
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 10 ac.setRotationAngle(ac.rotationAngle() + adj);
755 10 ac.setRays(ac.ray(Y), ac.ray(X));
756 }
757
758 // Do the actual comparison by computing four points on each ellipse.
759 28 Point tps[] = {Point(1,0), Point(0,1), Point(-1,0), Point(0,-1)};
760
2/2
✓ Branch 0 taken 88 times.
✓ Branch 1 taken 20 times.
 108 for (auto & tp : tps) {
761
5/8
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 88 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 88 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 8 times.
✓ Branch 15 taken 80 times.
 264 if (!are_near(tp * ac.unitCircleTransform(),
762
2/4
✓ Branch 2 taken 88 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 88 times.
✗ Branch 6 not taken.
 176 tp * bc.unitCircleTransform(),
763 precision))
764 8 return false;
765 }
766 20 return true;
767 }
768
769 ✗ std::ostream &operator<<(std::ostream &out, Ellipse const &e)
770 {
771 ✗ out << "Ellipse(" << e.center() << ", " << e.rays()
772 ✗ << ", " << format_coord_nice(e.rotationAngle()) << ")";
773 ✗ return out;
774 }
775
776 } // end namespace Geom
777
778
779 /*
780 Local Variables:
781 mode:c++
782 c-file-style:"stroustrup"
783 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
784 indent-tabs-mode:nil
785 fill-column:99
786 End:
787 */
788 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
789
790
791