GCC Code Coverage Report

Directory: ./
File: src/2geom/elliptical-arc.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 404 504 80.2%
Functions: 26 32 81.2%
Branches: 404 874 46.2%
Line Branch Exec Source
1 /*
2 * SVG Elliptical Arc Class
3 *
4 * Authors:
5 * Marco Cecchetti <mrcekets at gmail.com>
6 * Krzysztof Kosiński <tweenk.pl@gmail.com>
7 * Copyright 2008-2009 Authors
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it either under the terms of the GNU Lesser General Public
11 * License version 2.1 as published by the Free Software Foundation
12 * (the "LGPL") or, at your option, under the terms of the Mozilla
13 * Public License Version 1.1 (the "MPL"). If you do not alter this
14 * notice, a recipient may use your version of this file under either
15 * the MPL or the LGPL.
16 *
17 * You should have received a copy of the LGPL along with this library
18 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 * You should have received a copy of the MPL along with this library
21 * in the file COPYING-MPL-1.1
22 *
23 * The contents of this file are subject to the Mozilla Public License
24 * Version 1.1 (the "License"); you may not use this file except in
25 * compliance with the License. You may obtain a copy of the License at
26 * http://www.mozilla.org/MPL/
27 *
28 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
29 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
30 * the specific language governing rights and limitations.
31 */
32
33 #include <cfloat>
34 #include <limits>
35 #include <memory>
36
37 #include <2geom/bezier-curve.h>
38 #include <2geom/ellipse.h>
39 #include <2geom/elliptical-arc.h>
40 #include <2geom/path-sink.h>
41 #include <2geom/sbasis-geometric.h>
42 #include <2geom/transforms.h>
43 #include <2geom/utils.h>
44
45 #include <2geom/numeric/vector.h>
46 #include <2geom/numeric/fitting-tool.h>
47 #include <2geom/numeric/fitting-model.h>
48
49 namespace Geom
50 {
51
52 /**
53 * @class EllipticalArc
54 * @brief Elliptical arc curve
55 *
56 * Elliptical arc is a curve taking the shape of a section of an ellipse.
57 *
58 * The arc function has two forms: the regular one, mapping the unit interval to points
59 * on 2D plane (the linear domain), and a second form that maps some interval
60 * \f$A \subseteq [0,2\pi)\f$ to the same points (the angular domain). The interval \f$A\f$
61 * determines which part of the ellipse forms the arc. The arc is said to contain an angle
62 * if its angular domain includes that angle (and therefore it is defined for that angle).
63 *
64 * The angular domain considers each ellipse to be
65 * a rotated, scaled and translated unit circle: 0 corresponds to \f$(1,0)\f$ on the unit circle,
66 * \f$\pi/2\f$ corresponds to \f$(0,1)\f$, \f$\pi\f$ to \f$(-1,0)\f$ and \f$3\pi/2\f$
67 * to \f$(0,-1)\f$. After the angle is mapped to a point from a unit circle, the point is
68 * transformed using a matrix of this form
69 * \f[ M = \left[ \begin{array}{ccc}
70 r_X \cos(\theta) & -r_Y \sin(\theta) & 0 \\
71 r_X \sin(\theta) & r_Y \cos(\theta) & 0 \\
72 c_X & c_Y & 1 \end{array} \right] \f]
73 * where \f$r_X, r_Y\f$ are the X and Y rays of the ellipse, \f$\theta\f$ is its angle of rotation,
74 * and \f$c_X, c_Y\f$ the coordinates of the ellipse's center - thus mapping the angle
75 * to some point on the ellipse. Note that for example the point at angluar coordinate 0,
76 * the center and the point at angular coordinate \f$\pi/4\f$ do not necessarily
77 * create an angle of \f$\pi/4\f$ radians; it is only the case if both axes of the ellipse
78 * are of the same length (i.e. it is a circle).
79 *
80 * @image html ellipse-angular-coordinates.png "An illustration of the angular domain"
81 *
82 * Each arc is defined by five variables: The initial and final point, the ellipse's rays,
83 * and the ellipse's rotation. Each set of those parameters corresponds to four different arcs,
84 * with two of them larger than half an ellipse and two of them turning clockwise while traveling
85 * from initial to final point. The two flags disambiguate between them: "large arc flag" selects
86 * the bigger arc, while the "sweep flag" selects the arc going in the direction of positive
87 * angles. Angles always increase when going from the +X axis in the direction of the +Y axis,
88 * so if Y grows downwards, this means clockwise.
89 *
90 * @image html elliptical-arc-flags.png "Meaning of arc flags (Y grows downwards)"
91 *
92 * @ingroup Curves
93 */
94
95
96 /** @brief Compute bounds of an elliptical arc.
97 * The bounds computation works as follows. The extreme X and Y points
98 * are either the endpoints or local minima / maxima of the ellipse.
99 * We already have endpoints, and we compute the local extremes.
100 * The local extremes correspond to two angles separated by \f$\pi\f$.
101 * Once we compute these angles, we check whether they belong to the arc,
102 * and if they do, we evaluate the ellipse at these angles.
103 * The bounding box of the arc is equal to the bounding box of the endpoints
104 * and the local extrema that belong to the arc.
105 */
106 42468 Rect EllipticalArc::boundsExact() const
107 {
108
2/2
✓ Branch 1 taken 1200 times.
✓ Branch 2 taken 41268 times.
 42468 if (isChord()) {
109
1/2
✓ Branch 1 taken 1200 times.
✗ Branch 2 not taken.
 1200 return { _initial_point, _final_point };
110 }
111
112
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 41268 times.
 41268 if (_angles.isFull()) {
113 ✗ return _ellipse.boundsExact();
114 }
115
116
1/2
✓ Branch 2 taken 41268 times.
✗ Branch 3 not taken.
 41268 auto const trans = unitCircleTransform();
117
118 41268 auto proj_bounds = [&] (Dim2 d) {
119 // The dth coordinate function pulls back to trans[d] * x + trans[d + 2] * y + trans[d + 4]
120 // in the coordinate system where the ellipse is a unit circle. We compute its range of
121 // values on the unit circle arc.
122
1/2
✓ Branch 4 taken 82536 times.
✗ Branch 5 not taken.
 82536 auto result = Interval(_initial_point[d], _final_point[d]);
123
124 82536 auto const v = Point(trans[d], trans[d + 2]);
125 82536 auto const r = v.length();
126 82536 auto const mid = trans[d + 4];
127
1/2
✓ Branch 2 taken 82536 times.
✗ Branch 3 not taken.
 82536 auto const angle = Angle(v);
128
129
2/2
✓ Branch 1 taken 30878 times.
✓ Branch 2 taken 51658 times.
 82536 if (_angles.contains(angle)) {
130 30878 result.expandTo(mid + r);
131 }
132
3/4
✓ Branch 2 taken 82536 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 40868 times.
✓ Branch 7 taken 41668 times.
 82536 if (_angles.contains(angle + M_PI)) {
133 40868 result.expandTo(mid - r);
134 }
135
136 165072 return result;
137 41268 };
138
139
3/6
✓ Branch 2 taken 41268 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 41268 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 41268 times.
✗ Branch 10 not taken.
 82536 return { proj_bounds(X), proj_bounds(Y) };
140 }
141
142 2400 void EllipticalArc::expandToTransformed(Rect &bbox, Affine const &transform) const
143 {
144
1/2
✓ Branch 2 taken 2400 times.
✗ Branch 3 not taken.
 2400 bbox.expandTo(_final_point * transform);
145
146
2/2
✓ Branch 1 taken 1200 times.
✓ Branch 2 taken 1200 times.
 2400 if (isChord()) {
147 1200 return;
148 }
149
150
2/4
✓ Branch 3 taken 1200 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1200 times.
✗ Branch 7 not taken.
 1200 auto const trans = unitCircleTransform() * transform;
151
152
2/2
✓ Branch 4 taken 2400 times.
✓ Branch 5 taken 1200 times.
 3600 for (auto d : { X, Y }) {
153 // See boundsExact() for explanation.
154 2400 auto const v = Point(trans[d], trans[d + 2]);
155 2400 auto const r = v.length();
156 2400 auto const mid = trans[d + 4];
157
1/2
✓ Branch 2 taken 2400 times.
✗ Branch 3 not taken.
 2400 auto const interval = Interval(mid - r, mid + r);
158
159
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2400 times.
 2400 if (_angles.isFull()) {
160 ✗ bbox[d].unionWith(interval);
161
2/2
✓ Branch 2 taken 2349 times.
✓ Branch 3 taken 51 times.
 2400 } else if (!bbox[d].contains(interval)) {
162
1/2
✓ Branch 2 taken 2349 times.
✗ Branch 3 not taken.
 2349 auto const angle = Angle(v);
163
2/2
✓ Branch 1 taken 748 times.
✓ Branch 2 taken 1601 times.
 2349 if (_angles.contains(angle)) {
164 748 bbox[d].expandTo(mid + r);
165 }
166
3/4
✓ Branch 2 taken 2349 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 751 times.
✓ Branch 7 taken 1598 times.
 2349 if (_angles.contains(angle + M_PI)) {
167 751 bbox[d].expandTo(mid - r);
168 }
169 }
170 }
171 }
172
173 20000 Point EllipticalArc::pointAtAngle(Coord t) const
174 {
175
1/2
✓ Branch 2 taken 20000 times.
✗ Branch 3 not taken.
 20000 Point ret = _ellipse.pointAt(t);
176 20000 return ret;
177 }
178
179 ✗ Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const
180 {
181 ✗ return _ellipse.valueAt(t, d);
182 }
183
184 16 std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
185 {
186 16 std::vector<Coord> sol;
187
188
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
 16 if (isChord()) {
189 ✗ sol = chord().roots(v, d);
190 ✗ return sol;
191 }
192
193 16 Interval unit_interval(0, 1);
194
195 16 double rotx, roty;
196
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 10 times.
 16 if (d == X) {
197
1/2
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
 6 sincos(rotationAngle(), roty, rotx);
198 6 roty = -roty;
199 } else {
200
1/2
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 10 sincos(rotationAngle(), rotx, roty);
201 }
202
203 16 double rxrotx = ray(X) * rotx;
204 16 double c_v = center(d) - v;
205
206 16 double a = -rxrotx + c_v;
207 16 double b = ray(Y) * roty;
208 16 double c = rxrotx + c_v;
209 //std::cerr << "a = " << a << std::endl;
210 //std::cerr << "b = " << b << std::endl;
211 //std::cerr << "c = " << c << std::endl;
212
213
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 16 if (a == 0)
214 {
215
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 sol.push_back(M_PI);
216
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 8 if (b != 0)
217 {
218 8 double s = 2 * std::atan(-c/(2*b));
219
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 6 times.
 8 if ( s < 0 ) s += 2*M_PI;
220
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 sol.push_back(s);
221 }
222 }
223 else
224 {
225 8 double delta = b * b - a * c;
226 //std::cerr << "delta = " << delta << std::endl;
227
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if (delta == 0) {
228 ✗ double s = 2 * std::atan(-b/a);
229 ✗ if ( s < 0 ) s += 2*M_PI;
230 ✗ sol.push_back(s);
231 }
232
1/2
✓ Branch 0 taken 8 times.
✗ Branch 1 not taken.
 8 else if ( delta > 0 )
233 {
234 8 double sq = std::sqrt(delta);
235 8 double s = 2 * std::atan( (-b - sq) / a );
236
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if ( s < 0 ) s += 2*M_PI;
237
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 sol.push_back(s);
238 8 s = 2 * std::atan( (-b + sq) / a );
239
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
 8 if ( s < 0 ) s += 2*M_PI;
240
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 sol.push_back(s);
241 }
242 }
243
244 16 std::vector<double> arc_sol;
245
2/2
✓ Branch 7 taken 32 times.
✓ Branch 8 taken 16 times.
 48 for (double & i : sol) {
246 //std::cerr << "s = " << deg_from_rad(sol[i]);
247
2/4
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 32 times.
✗ Branch 6 not taken.
 32 i = timeAtAngle(i);
248 //std::cerr << " -> t: " << sol[i] << std::endl;
249
2/2
✓ Branch 1 taken 25 times.
✓ Branch 2 taken 7 times.
 32 if (unit_interval.contains(i)) {
250
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 arc_sol.push_back(i);
251 }
252 }
253 16 return arc_sol;
254 16 }
255
256
257 // D(E(t,C),t) = E(t+PI/2,O), where C is the ellipse center
258 // the derivative doesn't rotate the ellipse but there is a translation
259 // of the parameter t by an angle of PI/2 so the ellipse points are shifted
260 // of such an angle in the cw direction
261 8 Curve *EllipticalArc::derivative() const
262 {
263
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
 8 if (isChord()) {
264 ✗ return chord().derivative();
265 }
266
267 8 EllipticalArc *result = static_cast<EllipticalArc*>(duplicate());
268 8 result->_ellipse.setCenter(0, 0);
269
1/2
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 result->_angles.setInitial(result->_angles.initialAngle() + M_PI/2);
270
1/2
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 result->_angles.setFinal(result->_angles.finalAngle() + M_PI/2);
271
2/4
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
 8 result->_initial_point = result->pointAtAngle( result->initialAngle() );
272
2/4
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
 8 result->_final_point = result->pointAtAngle( result->finalAngle() );
273 8 return result;
274 }
275
276
277 std::vector<Point>
278 34 EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
279 {
280
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 34 times.
 34 if (isChord()) {
281 ✗ return chord().pointAndDerivatives(t, n);
282 }
283
284 34 unsigned int nn = n+1; // nn represents the size of the result vector.
285 34 std::vector<Point> result;
286
1/2
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
 34 result.reserve(nn);
287
2/4
✓ Branch 2 taken 34 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 34 times.
✗ Branch 6 not taken.
 34 double angle = angleAt(t);
288
1/2
✓ Branch 2 taken 34 times.
✗ Branch 3 not taken.
 34 std::unique_ptr<EllipticalArc> ea( static_cast<EllipticalArc*>(duplicate()) );
289 34 ea->_ellipse.setCenter(0, 0);
290 34 unsigned int m = std::min(nn, 4u);
291
2/2
✓ Branch 0 taken 136 times.
✓ Branch 1 taken 34 times.
 170 for ( unsigned int i = 0; i < m; ++i )
292 {
293
2/4
✓ Branch 3 taken 136 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 136 times.
✗ Branch 7 not taken.
 136 result.push_back( ea->pointAtAngle(angle) );
294
2/2
✓ Branch 1 taken 76 times.
✓ Branch 2 taken 60 times.
 136 angle += (sweep() ? M_PI/2 : -M_PI/2);
295
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 122 times.
 136 if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
296 }
297 34 m = nn / 4;
298
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 34 times.
 34 for ( unsigned int i = 1; i < m; ++i )
299 {
300 ✗ for ( unsigned int j = 0; j < 4; ++j )
301 ✗ result.push_back( result[j] );
302 }
303 34 m = nn - 4 * m;
304
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 34 times.
 34 for ( unsigned int i = 0; i < m; ++i )
305 {
306 ✗ result.push_back( result[i] );
307 }
308
1/2
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
 34 if ( !result.empty() ) // nn != 0
309
1/2
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
 34 result[0] = pointAtAngle(angle);
310 34 return result;
311 34 }
312
313 100664 Point EllipticalArc::pointAt(Coord t) const
314 {
315
2/2
✓ Branch 0 taken 93 times.
✓ Branch 1 taken 100571 times.
 100664 if (t == 0.0) {
316 93 return initialPoint();
317 }
318
2/2
✓ Branch 0 taken 78 times.
✓ Branch 1 taken 100493 times.
 100571 if (t == 1.0) {
319 78 return finalPoint();
320 }
321
2/2
✓ Branch 1 taken 80000 times.
✓ Branch 2 taken 20493 times.
 100493 if (isChord()) {
322
2/4
✓ Branch 2 taken 80000 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 80000 times.
✗ Branch 6 not taken.
 160000 return chord().pointAt(t);
323 }
324
3/6
✓ Branch 2 taken 20493 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 20493 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 20493 times.
✗ Branch 9 not taken.
 40986 return _ellipse.pointAt(angleAt(t));
325 }
326
327 ✗ Coord EllipticalArc::valueAt(Coord t, Dim2 d) const
328 {
329 ✗ if (isChord()) return chord().valueAt(t, d);
330 ✗ return valueAtAngle(angleAt(t), d);
331 }
332
333 63 Curve* EllipticalArc::portion(double f, double t) const
334 {
335 // fix input arguments
336 63 f = std::clamp(f, 0.0, 1.0);
337 63 t = std::clamp(t, 0.0, 1.0);
338
339
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 63 times.
 63 if (f == t) {
340 ✗ EllipticalArc *arc = new EllipticalArc();
341 ✗ arc->_initial_point = arc->_final_point = pointAt(f);
342 ✗ return arc;
343 }
344
4/4
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 19 times.
✓ Branch 3 taken 32 times.
 63 if (f == 0.0 && t == 1.0) {
345 19 return duplicate();
346 }
347
3/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 43 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 44 if (f == 1.0 && t == 0.0) {
348 ✗ return reverse();
349 }
350
351 44 EllipticalArc *arc = static_cast<EllipticalArc*>(duplicate());
352 44 arc->_initial_point = pointAt(f);
353 44 arc->_final_point = pointAt(t);
354 44 arc->_angles.setAngles(angleAt(f), angleAt(t));
355
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 41 times.
 44 if (f > t) arc->_angles.setSweep(!sweep());
356
5/6
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 21 times.
✓ Branch 3 taken 23 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 23 times.
✓ Branch 6 taken 21 times.
 44 if ( _large_arc && fabs(angularExtent() * (t-f)) <= M_PI) {
357 23 arc->_large_arc = false;
358 }
359 44 return arc;
360 }
361
362 // the arc is the same but traversed in the opposite direction
363 40016 Curve *EllipticalArc::reverse() const
364 {
365 using std::swap;
366 40016 EllipticalArc *rarc = static_cast<EllipticalArc*>(duplicate());
367 40016 rarc->_angles.reverse();
368 40016 swap(rarc->_initial_point, rarc->_final_point);
369 40016 return rarc;
370 }
371
372 #ifdef HAVE_GSL // GSL is required for function "solve_reals"
373 37 std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from, double to ) const
374 {
375 37 std::vector<double> result;
376
377
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 37 times.
 37 if ( from > to ) std::swap(from, to);
378
2/4
✓ Branch 0 taken 37 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 37 times.
 37 if ( from < 0 || to > 1 )
379 {
380 ✗ THROW_RANGEERROR("[from,to] interval out of range");
381 }
382
383
3/8
✗ Branch 2 not taken.
✓ Branch 3 taken 37 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 37 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 37 times.
 37 if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )
384 {
385 ✗ result.push_back(from);
386 ✗ return result;
387 }
388
3/6
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 37 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 37 times.
 37 else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
389 {
390 ✗ LineSegment seg(pointAt(from), pointAt(to));
391 ✗ Point np = seg.pointAt( seg.nearestTime(p) );
392 ✗ if ( are_near(ray(Y), 0) )
393 {
394 ✗ if ( are_near(rotationAngle(), M_PI/2)
395 ✗ || are_near(rotationAngle(), 3*M_PI/2) )
396 {
397 ✗ result = roots(np[Y], Y);
398 }
399 else
400 {
401 ✗ result = roots(np[X], X);
402 }
403 }
404 else
405 {
406 ✗ if ( are_near(rotationAngle(), M_PI/2)
407 ✗ || are_near(rotationAngle(), 3*M_PI/2) )
408 {
409 ✗ result = roots(np[X], X);
410 }
411 else
412 {
413 ✗ result = roots(np[Y], Y);
414 }
415 }
416 ✗ return result;
417 ✗ }
418
2/2
✓ Branch 3 taken 25 times.
✓ Branch 4 taken 12 times.
 37 else if ( are_near(ray(X), ray(Y)) )
419 {
420 25 Point r = p - center();
421
2/4
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 25 times.
 25 if ( are_near(r, Point(0,0)) )
422 {
423 ✗ THROW_INFINITESOLUTIONS(0);
424 }
425 // TODO: implement case r != 0
426 // Point np = ray(X) * unit_vector(r);
427 // std::vector<double> solX = roots(np[X],X);
428 // std::vector<double> solY = roots(np[Y],Y);
429 // double t;
430 // if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
431 // {
432 // t = solX[0];
433 // }
434 // else
435 // {
436 // t = solX[1];
437 // }
438 // if ( !(t < from || t > to) )
439 // {
440 // result.push_back(t);
441 // }
442 // else
443 // {
444 //
445 // }
446 }
447
448 // solve the equation <D(E(t),t)|E(t)-p> == 0
449 // that provides min and max distance points
450 // on the ellipse E wrt the point p
451 // after the substitutions:
452 // cos(t) = (1 - s^2) / (1 + s^2)
453 // sin(t) = 2t / (1 + s^2)
454 // where s = tan(t/2)
455 // we get a 4th degree equation in s
456 /*
457 * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
458 * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
459 * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
460 * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
461 */
462
463 37 Point p_c = p - center();
464 37 double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
465 37 double sinrot, cosrot;
466
1/2
✓ Branch 3 taken 37 times.
✗ Branch 4 not taken.
 37 sincos(rotationAngle(), sinrot, cosrot);
467 37 double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
468 37 Poly coeff;
469
1/2
✓ Branch 1 taken 37 times.
✗ Branch 2 not taken.
 37 coeff.resize(5);
470 37 coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
471 37 coeff[3] = 2 * ( rx2_ry2 + expr1 );
472 37 coeff[2] = 0;
473 37 coeff[1] = 2 * ( -rx2_ry2 + expr1 );
474 37 coeff[0] = -coeff[4];
475
476 // for ( unsigned int i = 0; i < 5; ++i )
477 // std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
478
479 37 std::vector<double> real_sol;
480 // gsl_poly_complex_solve raises an error
481 // if the leading coefficient is zero
482
2/2
✓ Branch 2 taken 7 times.
✓ Branch 3 taken 30 times.
 37 if ( are_near(coeff[4], 0) )
483 {
484
1/2
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
 7 real_sol.push_back(0);
485
1/2
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
 7 if ( !are_near(coeff[3], 0) )
486 {
487 7 double sq = -coeff[1] / coeff[3];
488
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 if ( sq > 0 )
489 {
490 ✗ double s = std::sqrt(sq);
491 ✗ real_sol.push_back(s);
492 ✗ real_sol.push_back(-s);
493 }
494 }
495 }
496 else
497 {
498
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 real_sol = solve_reals(coeff);
499 }
500
501
2/2
✓ Branch 7 taken 77 times.
✓ Branch 8 taken 37 times.
 114 for (double & i : real_sol)
502 {
503 77 i = 2 * std::atan(i);
504
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 37 times.
 77 if ( i < 0 ) i += 2*M_PI;
505 }
506 // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
507 // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
508
2/2
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 30 times.
 37 if ( (real_sol.size() % 2) != 0 )
509 {
510
1/2
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
 7 real_sol.push_back(M_PI);
511 }
512
513 37 double mindistsq1 = std::numeric_limits<double>::max();
514 37 double mindistsq2 = std::numeric_limits<double>::max();
515 37 double dsq = 0;
516 37 unsigned int mi1 = 0, mi2 = 0;
517
2/2
✓ Branch 1 taken 84 times.
✓ Branch 2 taken 37 times.
 121 for ( unsigned int i = 0; i < real_sol.size(); ++i )
518 {
519
2/4
✓ Branch 3 taken 84 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 84 times.
✗ Branch 7 not taken.
 84 dsq = distanceSq(p, pointAtAngle(real_sol[i]));
520
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 33 times.
 84 if ( mindistsq1 > dsq )
521 {
522 51 mindistsq2 = mindistsq1;
523 51 mi2 = mi1;
524 51 mindistsq1 = dsq;
525 51 mi1 = i;
526 }
527
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 10 times.
 33 else if ( mindistsq2 > dsq )
528 {
529 23 mindistsq2 = dsq;
530 23 mi2 = i;
531 }
532 }
533
534
2/4
✓ Branch 4 taken 37 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 37 times.
✗ Branch 8 not taken.
 37 double t = timeAtAngle(real_sol[mi1]);
535
3/4
✓ Branch 0 taken 37 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 33 times.
✓ Branch 3 taken 4 times.
 37 if ( !(t < from || t > to) )
536 {
537
1/2
✓ Branch 1 taken 33 times.
✗ Branch 2 not taken.
 33 result.push_back(t);
538 }
539
540 37 bool second_sol = false;
541
2/4
✓ Branch 3 taken 37 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 37 times.
✗ Branch 7 not taken.
 37 t = timeAtAngle(real_sol[mi2]);
542
5/8
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 32 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 37 times.
 37 if ( real_sol.size() == 4 && !(t < from || t > to) )
543 {
544 ✗ if ( result.empty() || are_near(mindistsq1, mindistsq2) )
545 {
546 ✗ result.push_back(t);
547 ✗ second_sol = true;
548 }
549 }
550
551 // we need to test extreme points too
552
2/4
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 37 times.
✗ Branch 6 not taken.
 37 double dsq1 = distanceSq(p, pointAt(from));
553
2/4
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 37 times.
✗ Branch 6 not taken.
 37 double dsq2 = distanceSq(p, pointAt(to));
554
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 37 times.
 37 if ( second_sol )
555 {
556 ✗ if ( mindistsq2 > dsq1 )
557 {
558 ✗ result.clear();
559 ✗ result.push_back(from);
560 ✗ mindistsq2 = dsq1;
561 }
562 ✗ else if ( are_near(mindistsq2, dsq) )
563 {
564 ✗ result.push_back(from);
565 }
566 ✗ if ( mindistsq2 > dsq2 )
567 {
568 ✗ result.clear();
569 ✗ result.push_back(to);
570 }
571 ✗ else if ( are_near(mindistsq2, dsq2) )
572 {
573 ✗ result.push_back(to);
574 }
575
576 }
577 else
578 {
579
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 33 times.
 37 if ( result.empty() )
580 {
581
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( are_near(dsq1, dsq2) )
582 {
583 ✗ result.push_back(from);
584 ✗ result.push_back(to);
585 }
586
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 else if ( dsq2 > dsq1 )
587 {
588 ✗ result.push_back(from);
589 }
590 else
591 {
592
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 result.push_back(to);
593 }
594 }
595 }
596
597 37 return result;
598 37 }
599 #endif
600
601 /** @brief Convert the passed intersections to curve time parametrization
602 * and filter out any invalid intersections.
603 */
604 10021 std::vector<ShapeIntersection> EllipticalArc::_filterIntersections(std::vector<ShapeIntersection> &&xs,
605 bool is_first) const
606 {
607 10021 std::vector<ShapeIntersection> result;
608
1/2
✓ Branch 2 taken 10021 times.
✗ Branch 3 not taken.
 10021 result.reserve(xs.size());
609
2/2
✓ Branch 7 taken 20038 times.
✓ Branch 8 taken 10021 times.
 30059 for (auto &xing : xs) {
610
3/4
✓ Branch 1 taken 20038 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10036 times.
✓ Branch 4 taken 10002 times.
 20038 if (_validateIntersection(xing, is_first)) {
611
1/2
✓ Branch 2 taken 10036 times.
✗ Branch 3 not taken.
 10036 result.emplace_back(std::move(xing));
612 }
613 }
614 10021 return result;
615 ✗ }
616
617 /** @brief Convert the passed intersection to curve time and check whether the intersection
618 * is numerically sane.
619 *
620 * @param xing The intersection to convert to curve time and to validate.
621 * @param is_first If true, this arc is the first of the intersected curves; if false, it's second.
622 * @return Whether the intersection is valid.
623 *
624 * Note that the intersection is guaranteed to be converted only if the return value is true.
625 */
626 20038 bool EllipticalArc::_validateIntersection(ShapeIntersection &xing, bool is_first) const
627 {
628 static auto const UNIT_INTERVAL = Interval(0, 1);
629 20038 constexpr auto EPS = 1e-4;
630
631
2/2
✓ Branch 0 taken 20030 times.
✓ Branch 1 taken 8 times.
 20038 Coord &t = is_first ? xing.first : xing.second;
632
4/6
✓ Branch 4 taken 20038 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20038 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 20037 times.
 20038 if (!are_near_rel(_ellipse.pointAt(t), xing.point(), EPS)) {
633 1 return false;
634 }
635
636
2/4
✓ Branch 2 taken 20037 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 20037 times.
✗ Branch 6 not taken.
 20037 t = timeAtAngle(t);
637
2/2
✓ Branch 1 taken 10001 times.
✓ Branch 2 taken 10036 times.
 20037 if (!UNIT_INTERVAL.contains(t)) {
638 10001 return false;
639 }
640
3/6
✓ Branch 4 taken 10036 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10036 times.
✗ Branch 8 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 10036 times.
 10036 if (!are_near_rel(pointAt(t), xing.point(), EPS)) {
641 ✗ return false;
642 }
643 10036 return true;
644 }
645
646 90026 std::vector<CurveIntersection> EllipticalArc::intersect(Curve const &other, Coord eps) const
647 {
648
3/4
✓ Branch 1 taken 90026 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 80000 times.
✓ Branch 4 taken 10026 times.
 90026 if (isLineSegment()) {
649
1/2
✓ Branch 2 taken 80000 times.
✗ Branch 3 not taken.
 80000 LineSegment ls(_initial_point, _final_point);
650
1/2
✓ Branch 1 taken 80000 times.
✗ Branch 2 not taken.
 80000 return ls.intersect(other, eps);
651 80000 }
652
653
3/4
✓ Branch 1 taken 10026 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10010 times.
✓ Branch 4 taken 16 times.
 10026 if (other.isLineSegment()) {
654
3/6
✓ Branch 2 taken 10010 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 10010 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 10010 times.
✗ Branch 9 not taken.
 10010 LineSegment ls(other.initialPoint(), other.finalPoint());
655
2/4
✓ Branch 2 taken 10010 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 10010 times.
✗ Branch 6 not taken.
 20020 return _filterIntersections(_ellipse.intersect(ls), true);
656 10010 }
657
658
3/4
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 13 times.
 16 if (auto bez = dynamic_cast<BezierCurve const *>(&other)) {
659
2/4
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
 6 return _filterIntersections(_ellipse.intersect(bez->fragment()), true);
660 }
661
662
2/4
✓ Branch 0 taken 13 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
 13 if (auto arc = dynamic_cast<EllipticalArc const *>(&other)) {
663 13 std::vector<CurveIntersection> crossings;
664 try {
665
2/2
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 9 times.
 13 crossings = _ellipse.intersect(arc->_ellipse);
666
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 } catch (InfinitelyManySolutions &) {
667 // This could happen if the two arcs come from the same ellipse.
668
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 return _intersectSameEllipse(arc);
669
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 }
670
2/4
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 8 return arc->_filterIntersections(_filterIntersections(std::move(crossings), true), false);
671 13 }
672
673 // in case someone wants to make a custom curve type
674 ✗ auto result = other.intersect(*this, eps);
675 ✗ transpose_in_place(result);
676 ✗ return result;
677 ✗ }
678
679 /** @brief Check if two arcs on the same ellipse intersect/overlap.
680 *
681 * @param other Another elliptical arc on the same ellipse as this one.
682 * @return If the arcs overlap, the returned vector contains synthesized intersections
683 * at the start and end of the overlap.
684 * If the arcs do not overlap, an empty vector is returned.
685 */
686 9 std::vector<ShapeIntersection> EllipticalArc::_intersectSameEllipse(EllipticalArc const *other) const
687 {
688
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
 9 assert(_ellipse == other->_ellipse);
689 9 auto const &other_angles = other->angularInterval();
690 9 std::vector<ShapeIntersection> result;
691
692 /// A closure to create an "intersection" at the prescribed angle.
693 9 auto const synthesize_intersection = [&](Angle angle) {
694
1/2
✓ Branch 2 taken 22 times.
✗ Branch 3 not taken.
 22 auto const time = timeAtAngle(angle);
695
3/4
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 11 taken 13 times.
✓ Branch 12 taken 9 times.
 22 if (result.end() == std::find_if(result.begin(), result.end(),
696 15 [=](ShapeIntersection const &xing) -> bool {
697 15 return xing.first == time;
698 }))
699 {
700
4/8
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 13 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 13 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 13 times.
✗ Branch 13 not taken.
 13 result.emplace_back(time, other->timeAtAngle(angle), _ellipse.pointAt(angle));
701 }
702 22 };
703
704
2/2
✓ Branch 7 taken 18 times.
✓ Branch 8 taken 9 times.
 27 for (auto a : {_angles.initialAngle(), _angles.finalAngle()}) {
705
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 7 times.
 18 if (other_angles.contains(a)) {
706
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 synthesize_intersection(a);
707 }
708 }
709
2/2
✓ Branch 7 taken 18 times.
✓ Branch 8 taken 9 times.
 27 for (auto a : {other_angles.initialAngle(), other_angles.finalAngle()}) {
710
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 7 times.
 18 if (_angles.contains(a)) {
711
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 synthesize_intersection(a);
712 }
713 }
714 9 return result;
715 ✗ }
716
717 10389 void EllipticalArc::_updateCenterAndAngles()
718 {
719 // See: http://www.w3.org/TR/SVG/implnote.html#ArcImplementationNotes
720
2/4
✓ Branch 3 taken 10389 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 10389 times.
✗ Branch 8 not taken.
 10389 Point d = initialPoint() - finalPoint();
721
3/6
✓ Branch 3 taken 10389 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 10389 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10389 times.
✗ Branch 11 not taken.
 10389 Point mid = middle_point(initialPoint(), finalPoint());
722
723 10389 auto degenerate_ellipse = [&] {
724 2 _ellipse = Ellipse();
725
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 _ellipse.setCenter(initialPoint());
726 2 _angles = AngleInterval();
727 2 _large_arc = false;
728 2 };
729
730 // if ip = fp, the arc contains no other points
731
4/6
✓ Branch 2 taken 10389 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 10389 times.
✗ Branch 7 not taken.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 10387 times.
 10389 if (initialPoint() == finalPoint()) {
732
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 degenerate_ellipse();
733 2 return;
734 }
735
736 // rays should be positive
737 10387 _ellipse.setRays(std::fabs(ray(X)), std::fabs(ray(Y)));
738
739
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 10383 times.
 10387 if (isChord()) {
740 4 _ellipse.setRays(L2(d) / 2, 0);
741
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 _ellipse.setRotationAngle(atan2(d));
742 4 _ellipse.setCenter(mid);
743
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 4 _angles.setAngles(0, M_PI);
744 4 _angles.setSweep(false);
745 4 _large_arc = false;
746 4 return;
747 }
748
749
2/4
✓ Branch 4 taken 10383 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10383 times.
✗ Branch 8 not taken.
 10383 Rotate rot(rotationAngle()); // the matrix in F.6.5.3
750
1/2
✓ Branch 2 taken 10383 times.
✗ Branch 3 not taken.
 10383 Rotate invrot = rot.inverse(); // the matrix in F.6.5.1
751
752 10383 Point r = rays();
753
1/2
✓ Branch 5 taken 10383 times.
✗ Branch 6 not taken.
 10383 Point p = d / 2 * invrot; // x', y' in F.6.5.1
754 10383 Point c(0,0); // cx', cy' in F.6.5.2
755
756 // Correct out-of-range radii
757 10383 Coord lambda = hypot(p[X]/r[X], p[Y]/r[Y]);
758
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 10331 times.
 10383 if (lambda > 1) {
759 52 r *= lambda;
760 52 _ellipse.setRays(r);
761 52 _ellipse.setCenter(mid);
762 } else {
763 // evaluate F.6.5.2
764 10331 Coord rxry = r[X]*r[X] * r[Y]*r[Y];
765 10331 Coord pxry = p[X]*p[X] * r[Y]*r[Y];
766 10331 Coord rxpy = r[X]*r[X] * p[Y]*p[Y];
767 10331 Coord const denominator = rxpy + pxry;
768
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10331 times.
 10331 if (denominator == 0.0) {
769 ✗ degenerate_ellipse();
770 ✗ return;
771 }
772 10331 Coord rad = (rxry - pxry - rxpy) / denominator;
773 // normally rad should never be negative, but numerical inaccuracy may cause this
774
2/2
✓ Branch 0 taken 10044 times.
✓ Branch 1 taken 287 times.
 10331 if (rad > 0) {
775 10044 rad = std::sqrt(rad);
776
2/2
✓ Branch 1 taken 4998 times.
✓ Branch 2 taken 5046 times.
 10044 if (sweep() == _large_arc) {
777 4998 rad = -rad;
778 }
779 10044 c = rad * Point(r[X]*p[Y]/r[Y], -r[Y]*p[X]/r[X]);
780
1/2
✓ Branch 3 taken 10044 times.
✗ Branch 4 not taken.
 10044 _ellipse.setCenter(c * rot + mid);
781 } else {
782 287 _ellipse.setCenter(mid);
783 }
784 }
785
786 // Compute start and end angles.
787 // If the ellipse was enlarged, c will be zero - this is correct.
788 10383 Point sp((p[X] - c[X]) / r[X], (p[Y] - c[Y]) / r[Y]);
789 10383 Point ep((-p[X] - c[X]) / r[X], (-p[Y] - c[Y]) / r[Y]);
790 10383 Point v(1, 0);
791
792
2/4
✓ Branch 2 taken 10383 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 10383 times.
✗ Branch 6 not taken.
 10383 _angles.setInitial(angle_between(v, sp));
793
2/4
✓ Branch 2 taken 10383 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 10383 times.
✗ Branch 6 not taken.
 10383 _angles.setFinal(angle_between(v, ep));
794
795 /*double sweep_angle = angle_between(sp, ep);
796 if (!sweep() && sweep_angle > 0) sweep_angle -= 2*M_PI;
797 if (sweep() && sweep_angle < 0) sweep_angle += 2*M_PI;*/
798 }
799
800 23 D2<SBasis> EllipticalArc::toSBasis() const
801 {
802
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 23 times.
 23 if (isChord()) {
803 ✗ return chord().toSBasis();
804 }
805
806
1/2
✓ Branch 2 taken 23 times.
✗ Branch 3 not taken.
 23 D2<SBasis> arc;
807 // the interval of parametrization has to be [0,1]
808 23 Coord et = initialAngle().radians() + sweepAngle();
809 23 Linear param(initialAngle().radians(), et);
810 23 Coord cosrot, sinrot;
811
1/2
✓ Branch 3 taken 23 times.
✗ Branch 4 not taken.
 23 sincos(rotationAngle(), sinrot, cosrot);
812
813 // order = 4 seems to be enough to get a perfect looking elliptical arc
814
2/4
✓ Branch 3 taken 23 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 23 times.
✗ Branch 8 not taken.
 23 SBasis arc_x = ray(X) * cos(param,4);
815
2/4
✓ Branch 3 taken 23 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 23 times.
✗ Branch 8 not taken.
 23 SBasis arc_y = ray(Y) * sin(param,4);
816
6/12
✓ Branch 7 taken 23 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 23 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 23 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 23 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 23 times.
✗ Branch 23 not taken.
✓ Branch 26 taken 23 times.
✗ Branch 27 not taken.
 23 arc[0] = arc_x * cosrot - arc_y * sinrot + Linear(center(X), center(X));
817
6/12
✓ Branch 7 taken 23 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 23 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 23 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 23 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 23 times.
✗ Branch 23 not taken.
✓ Branch 26 taken 23 times.
✗ Branch 27 not taken.
 23 arc[1] = arc_x * sinrot + arc_y * cosrot + Linear(center(Y), center(Y));
818
819 // ensure that endpoints remain exact
820
2/2
✓ Branch 0 taken 46 times.
✓ Branch 1 taken 23 times.
 69 for ( int d = 0 ; d < 2 ; d++ ) {
821
2/4
✓ Branch 2 taken 46 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 46 times.
✗ Branch 8 not taken.
 46 arc[d][0][0] = initialPoint()[d];
822
2/4
✓ Branch 2 taken 46 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 46 times.
✗ Branch 8 not taken.
 46 arc[d][0][1] = finalPoint()[d];
823 }
824
825
1/2
✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
 23 return arc;
826 23 }
827
828 // All operations that do not contain skew can be evaluated
829 // without passing through the implicit form of the ellipse,
830 // which preserves precision.
831
832 ✗ void EllipticalArc::operator*=(Translate const &tr)
833 {
834 ✗ _initial_point *= tr;
835 ✗ _final_point *= tr;
836 ✗ _ellipse *= tr;
837 ✗ }
838
839 4 void EllipticalArc::operator*=(Scale const &s)
840 {
841 4 _initial_point *= s;
842 4 _final_point *= s;
843 4 _ellipse *= s;
844 4 }
845
846 ✗ void EllipticalArc::operator*=(Rotate const &r)
847 {
848 ✗ _initial_point *= r;
849 ✗ _final_point *= r;
850 ✗ _ellipse *= r;
851 ✗ }
852
853 ✗ void EllipticalArc::operator*=(Zoom const &z)
854 {
855 ✗ _initial_point *= z;
856 ✗ _final_point *= z;
857 ✗ _ellipse *= z;
858 ✗ }
859
860 2401 void EllipticalArc::operator*=(Affine const& m)
861 {
862
2/2
✓ Branch 1 taken 1200 times.
✓ Branch 2 taken 1201 times.
 2401 if (isChord()) {
863 1200 _initial_point *= m;
864 1200 _final_point *= m;
865
1/2
✓ Branch 2 taken 1200 times.
✗ Branch 3 not taken.
 1200 _ellipse.setCenter(middle_point(_initial_point, _final_point));
866 1200 _ellipse.setRays(0, 0);
867
1/2
✓ Branch 2 taken 1200 times.
✗ Branch 3 not taken.
 1200 _ellipse.setRotationAngle(0);
868 1200 return;
869 }
870
871 1201 _initial_point *= m;
872 1201 _final_point *= m;
873 1201 _ellipse *= m;
874
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1201 times.
 1201 if (m.det() < 0) {
875 ✗ _angles.setSweep(!sweep());
876 }
877
878 // ellipse transformation does not preserve its functional form,
879 // i.e. e.pointAt(0.5)*m and (e*m).pointAt(0.5) can be different.
880 // We need to recompute start / end angles.
881
2/4
✓ Branch 2 taken 1201 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1201 times.
✗ Branch 6 not taken.
 1201 _angles.setInitial(_ellipse.timeAt(_initial_point));
882
2/4
✓ Branch 2 taken 1201 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1201 times.
✗ Branch 6 not taken.
 1201 _angles.setFinal(_ellipse.timeAt(_final_point));
883 }
884
885 15 bool EllipticalArc::_equalTo(Curve const &c) const
886 {
887
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 if (this == &c) return true;
888
1/2
✓ Branch 0 taken 15 times.
✗ Branch 1 not taken.
 15 auto other = dynamic_cast<EllipticalArc const *>(&c);
889
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 if (!other) return false;
890
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 15 times.
 15 if (_initial_point != other->_initial_point) return false;
891
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 15 times.
 15 if (_final_point != other->_final_point) return false;
892 // TODO: all arcs with ellipse rays which are too small
893 // and fall back to a line should probably be equal
894
1/2
✗ Branch 7 not taken.
✓ Branch 8 taken 15 times.
 15 if (rays() != other->rays()) return false;
895
1/2
✗ Branch 7 not taken.
✓ Branch 8 taken 15 times.
 15 if (rotationAngle() != other->rotationAngle()) return false;
896
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 if (_large_arc != other->_large_arc) return false;
897
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 15 times.
 15 if (sweep() != other->sweep()) return false;
898 15 return true;
899 }
900
901 13 bool EllipticalArc::isNear(Curve const &c, Coord precision) const
902 {
903
1/2
✓ Branch 0 taken 13 times.
✗ Branch 1 not taken.
 13 EllipticalArc const *other = dynamic_cast<EllipticalArc const *>(&c);
904
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 11 times.
 13 if (!other) {
905
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if (isChord()) {
906 ✗ return c.isNear(chord(), precision);
907 }
908 2 return false;
909 }
910
911
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 11 times.
 11 if (!are_near(_initial_point, other->_initial_point, precision)) return false;
912
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 10 times.
 11 if (!are_near(_final_point, other->_final_point, precision)) return false;
913
2/6
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 times.
 10 if (isChord() && other->isChord()) return true;
914
915
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
 10 if (sweep() != other->sweep()) return false;
916
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
 10 if (!are_near(_ellipse, other->_ellipse, precision)) return false;
917 10 return true;
918 }
919
920 169 void EllipticalArc::feed(PathSink &sink, bool moveto_initial) const
921 {
922
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 169 times.
 169 if (moveto_initial) {
923 ✗ sink.moveTo(_initial_point);
924 }
925
2/4
✓ Branch 4 taken 169 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 169 times.
✗ Branch 10 not taken.
 169 sink.arcTo(ray(X), ray(Y), rotationAngle(), _large_arc, sweep(), _final_point);
926 169 }
927
928 9863 int EllipticalArc::winding(Point const &p) const
929 {
930 using std::swap;
931
932 9863 double sinrot, cosrot;
933
1/2
✓ Branch 3 taken 9863 times.
✗ Branch 4 not taken.
 9863 sincos(rotationAngle(), sinrot, cosrot);
934
935
1/2
✓ Branch 4 taken 9863 times.
✗ Branch 5 not taken.
 9863 Angle ymin_a = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
936
1/2
✓ Branch 3 taken 9863 times.
✗ Branch 4 not taken.
 9863 Angle ymax_a = ymin_a + M_PI;
937
938
2/4
✓ Branch 2 taken 9863 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9863 times.
✗ Branch 6 not taken.
 9863 Point ymin = pointAtAngle(ymin_a);
939
2/4
✓ Branch 2 taken 9863 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9863 times.
✗ Branch 6 not taken.
 9863 Point ymax = pointAtAngle(ymax_a);
940
1/2
✓ Branch 2 taken 9863 times.
✗ Branch 3 not taken.
 9863 if (ymin[Y] > ymax[Y]) {
941 9863 swap(ymin, ymax);
942 9863 swap(ymin_a, ymax_a);
943 }
944
945
3/6
✓ Branch 4 taken 9863 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 9863 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 9863 times.
 9863 if (!Interval(ymin[Y], ymax[Y]).lowerContains(p[Y])) {
946 ✗ return 0;
947 }
948
949 9863 bool const left = cross(ymax - ymin, p - ymin) > 0;
950
1/2
✓ Branch 1 taken 9863 times.
✗ Branch 2 not taken.
 9863 bool const inside = _ellipse.contains(p);
951
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9863 times.
 9863 if (_angles.isFull()) {
952 ✗ if (inside) {
953 ✗ return sweep() ? 1 : -1;
954 }
955 ✗ return 0;
956 }
957 9863 bool const includes_ymin = _angles.contains(ymin_a);
958 9863 bool const includes_ymax = _angles.contains(ymax_a);
959
960 9863 AngleInterval rarc(ymin_a, ymax_a, true),
961 9863 larc(ymax_a, ymin_a, true);
962
963 // we'll compute the result for an arc in the direction of increasing angles
964 // and then negate if necessary
965 9863 Angle ia = initialAngle(), fa = finalAngle();
966 9863 Point ip = _initial_point, fp = _final_point;
967
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 9857 times.
 9863 if (!sweep()) {
968 6 swap(ia, fa);
969 6 swap(ip, fp);
970 }
971
972 9863 bool const initial_left = larc.contains(ia);
973 9863 bool const final_left = larc.contains(fa);
974
975 9863 bool intersects_left = false, intersects_right = false;
976
4/4
✓ Branch 0 taken 2153 times.
✓ Branch 1 taken 7710 times.
✓ Branch 2 taken 1060 times.
✓ Branch 3 taken 1093 times.
 9863 if (inside || left) {
977 // The point is inside the ellipse or to the left of it, so the rightwards horizontal ray
978 // may intersect the part of the arc contained in the right half of the ellipse.
979 // There are four ways in which this can happen.
980
981 1060 intersects_right =
982 // Possiblity 1: the arc extends into the right half through the min-Y point
983 // and the ray intersects this extension:
984
7/12
✓ Branch 0 taken 1967 times.
✓ Branch 1 taken 2435 times.
✓ Branch 6 taken 1967 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 1967 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1967 times.
✓ Branch 14 taken 6836 times.
✓ Branch 15 taken 1967 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 15172 (includes_ymin && !final_left && Interval(ymin[Y], fp[Y]).lowerContains(p[Y]))
985 6836 ||
986 // Possibility 2: the arc starts and ends within the right half (hence, it cannot be the
987 // "large arc") and the ray's Y-coordinate is within the Y-coordinate range of the arc:
988
4/16
✓ Branch 0 taken 1938 times.
✓ Branch 1 taken 4898 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1938 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 8803 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 15639 (!initial_left && !final_left && !largeArc() && Interval(ip[Y], fp[Y]).lowerContains(p[Y]))
989 6836 ||
990 // Possibility 3: the arc starts in the right half and continues through the max-Y
991 // point into the left half:
992
9/14
✓ Branch 0 taken 1938 times.
✓ Branch 1 taken 4898 times.
✓ Branch 2 taken 1938 times.
✗ Branch 3 not taken.
✓ Branch 8 taken 1938 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 1938 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✓ Branch 15 taken 1936 times.
✓ Branch 16 taken 1938 times.
✓ Branch 17 taken 6865 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 15639 (!initial_left && includes_ymax && Interval(ip[Y], ymax[Y]).lowerContains(p[Y]))
993
4/4
✓ Branch 0 taken 4402 times.
✓ Branch 1 taken 4401 times.
✓ Branch 2 taken 4898 times.
✓ Branch 3 taken 2 times.
 22504 ||
994 // Possibility 4: the entire right half of the ellipse is contained in the arc.
995
4/6
✓ Branch 0 taken 4898 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2433 times.
✓ Branch 3 taken 2465 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2433 times.
 4898 (initial_left && final_left && includes_ymin && includes_ymax);
996 }
997
4/4
✓ Branch 0 taken 4899 times.
✓ Branch 1 taken 4964 times.
✓ Branch 2 taken 1093 times.
✓ Branch 3 taken 3806 times.
 9863 if (left && !inside) {
998 // The point is to the left of the ellipse, so the rightwards horizontal ray
999 // may intersect the part of the arc contained in the left half of the ellipse.
1000 // There are four ways in which this can happen.
1001
1002 1093 intersects_left =
1003 // Possibility 1: the arc starts in the left half and continues through the min-Y
1004 // point into the right half:
1005
7/12
✓ Branch 0 taken 527 times.
✓ Branch 1 taken 1 times.
✓ Branch 6 taken 527 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 527 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 527 times.
✓ Branch 14 taken 566 times.
✓ Branch 15 taken 527 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 2148 (includes_ymin && initial_left && Interval(ymin[Y], ip[Y]).lowerContains(p[Y]))
1006
2/2
✓ Branch 0 taken 565 times.
✓ Branch 1 taken 1 times.
 566 ||
1007 // Possibility 2: the arc starts and ends within the left half (hence, it cannot be the
1008 // "large arc") and the ray's Y-coordinate is within the Y-coordinate range of the arc:
1009
7/14
✓ Branch 0 taken 565 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 565 times.
✗ Branch 4 not taken.
✓ Branch 9 taken 565 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 565 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 565 times.
✓ Branch 17 taken 565 times.
✓ Branch 18 taken 528 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 1658 (initial_left && final_left && !largeArc() && Interval(ip[Y], fp[Y]).lowerContains(p[Y]))
1010
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 ||
1011 // Possibility 3: the arc extends into the left half through the max-Y point
1012 // and the ray intersects this extension:
1013
6/12
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 1 times.
✓ Branch 15 taken 1092 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 1094 (final_left && includes_ymax && Interval(fp[Y], ymax[Y]).lowerContains(p[Y]))
1014
2/2
✓ Branch 0 taken 528 times.
✓ Branch 1 taken 565 times.
 2186 ||
1015 // Possibility 4: the entire left half of the ellipse is contained in the arc.
1016 ✗ (!initial_left && !final_left && includes_ymin && includes_ymax);
1017
1018 }
1019 9863 int const winding_assuming_increasing_angles = (int)intersects_right - (int)intersects_left;
1020
2/2
✓ Branch 1 taken 9857 times.
✓ Branch 2 taken 6 times.
 9863 return sweep() ? winding_assuming_increasing_angles : -winding_assuming_increasing_angles;
1021 }
1022
1023 ✗ std::ostream &operator<<(std::ostream &out, EllipticalArc const &ea)
1024 {
1025 ✗ out << "EllipticalArc("
1026 ✗ << ea.initialPoint() << ", "
1027 ✗ << format_coord_nice(ea.ray(X)) << ", " << format_coord_nice(ea.ray(Y)) << ", "
1028 ✗ << format_coord_nice(ea.rotationAngle()) << ", "
1029 ✗ << "large_arc=" << (ea.largeArc() ? "true" : "false") << ", "
1030 ✗ << "sweep=" << (ea.sweep() ? "true" : "false") << ", "
1031 ✗ << ea.finalPoint() << ")";
1032 ✗ return out;
1033 }
1034
1035 } // end namespace Geom
1036
1037 /*
1038 Local Variables:
1039 mode:c++
1040 c-file-style:"stroustrup"
1041 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1042 indent-tabs-mode:nil
1043 fill-column:99
1044 End:
1045 */
1046 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
1047
1048