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 can find the local extremes
100 * by computing a partial derivative with respect to the angle
101 * and equating that to zero:
102 * \f{align*}{
103 x &= r_x \cos \varphi \cos \theta - r_y \sin \varphi \sin \theta + c_x \\
104 \frac{\partial x}{\partial \theta} &= -r_x \cos \varphi \sin \theta - r_y \sin \varphi \cos \theta = 0 \\
105 \frac{\sin \theta}{\cos \theta} &= \tan\theta = -\frac{r_y \sin \varphi}{r_x \cos \varphi} \\
106 \theta &= \tan^{-1} \frac{-r_y \sin \varphi}{r_x \cos \varphi}
107 \f}
108 * The local extremes correspond to two angles separated by \f$\pi\f$.
109 * Once we compute these angles, we check whether they belong to the arc,
110 * and if they do, we evaluate the ellipse at these angles.
111 * The bounding box of the arc is equal to the bounding box of the endpoints
112 * and the local extrema that belong to the arc.
113 */
boundsExact() const114 Rect EllipticalArc::boundsExact() const
115 {
116 if (isChord()) {
117 return chord().boundsExact();
118 }
119
120 Coord coord[2][4] = {
121 { _initial_point[X], _final_point[X], 0, 0 },
122 { _initial_point[Y], _final_point[Y], 0, 0 }
123 };
124 int ncoord[2] = { 2, 2 };
125
126 Angle extremes[2][2];
127 double sinrot, cosrot;
128 sincos(rotationAngle(), sinrot, cosrot);
129
130 extremes[X][0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
131 extremes[X][1] = extremes[X][0] + M_PI;
132 extremes[Y][0] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
133 extremes[Y][1] = extremes[Y][0] + M_PI;
134
135 for (unsigned d = 0; d < 2; ++d) {
136 for (unsigned i = 0; i < 2; ++i) {
137 if (containsAngle(extremes[d][i])) {
138 coord[d][ncoord[d]++] = valueAtAngle(extremes[d][i], d ? Y : X);
139 }
140 }
141 }
142
143 Interval xival = Interval::from_range(coord[X], coord[X] + ncoord[X]);
144 Interval yival = Interval::from_range(coord[Y], coord[Y] + ncoord[Y]);
145 Rect result(xival, yival);
146 return result;
147 }
148
149
pointAtAngle(Coord t) const150 Point EllipticalArc::pointAtAngle(Coord t) const
151 {
152 Point ret = _ellipse.pointAt(t);
153 return ret;
154 }
155
valueAtAngle(Coord t,Dim2 d) const156 Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const
157 {
158 return _ellipse.valueAt(t, d);
159 }
160
roots(Coord v,Dim2 d) const161 std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const
162 {
163 std::vector<Coord> sol;
164
165 if (isChord()) {
166 sol = chord().roots(v, d);
167 return sol;
168 }
169
170 Interval unit_interval(0, 1);
171
172 double rotx, roty;
173 if (d == X) {
174 sincos(rotationAngle(), roty, rotx);
175 roty = -roty;
176 } else {
177 sincos(rotationAngle(), rotx, roty);
178 }
179
180 double rxrotx = ray(X) * rotx;
181 double c_v = center(d) - v;
182
183 double a = -rxrotx + c_v;
184 double b = ray(Y) * roty;
185 double c = rxrotx + c_v;
186 //std::cerr << "a = " << a << std::endl;
187 //std::cerr << "b = " << b << std::endl;
188 //std::cerr << "c = " << c << std::endl;
189
190 if (a == 0)
191 {
192 sol.push_back(M_PI);
193 if (b != 0)
194 {
195 double s = 2 * std::atan(-c/(2*b));
196 if ( s < 0 ) s += 2*M_PI;
197 sol.push_back(s);
198 }
199 }
200 else
201 {
202 double delta = b * b - a * c;
203 //std::cerr << "delta = " << delta << std::endl;
204 if (delta == 0) {
205 double s = 2 * std::atan(-b/a);
206 if ( s < 0 ) s += 2*M_PI;
207 sol.push_back(s);
208 }
209 else if ( delta > 0 )
210 {
211 double sq = std::sqrt(delta);
212 double s = 2 * std::atan( (-b - sq) / a );
213 if ( s < 0 ) s += 2*M_PI;
214 sol.push_back(s);
215 s = 2 * std::atan( (-b + sq) / a );
216 if ( s < 0 ) s += 2*M_PI;
217 sol.push_back(s);
218 }
219 }
220
221 std::vector<double> arc_sol;
222 for (double & i : sol) {
223 //std::cerr << "s = " << deg_from_rad(sol[i]);
224 i = timeAtAngle(i);
225 //std::cerr << " -> t: " << sol[i] << std::endl;
226 if (unit_interval.contains(i)) {
227 arc_sol.push_back(i);
228 }
229 }
230 return arc_sol;
231 }
232
233
234 // D(E(t,C),t) = E(t+PI/2,O), where C is the ellipse center
235 // the derivative doesn't rotate the ellipse but there is a translation
236 // of the parameter t by an angle of PI/2 so the ellipse points are shifted
237 // of such an angle in the cw direction
derivative() const238 Curve *EllipticalArc::derivative() const
239 {
240 if (isChord()) {
241 return chord().derivative();
242 }
243
244 EllipticalArc *result = static_cast<EllipticalArc*>(duplicate());
245 result->_ellipse.setCenter(0, 0);
246 result->_angles.setInitial(result->_angles.initialAngle() + M_PI/2);
247 result->_angles.setFinal(result->_angles.finalAngle() + M_PI/2);
248 result->_initial_point = result->pointAtAngle( result->initialAngle() );
249 result->_final_point = result->pointAtAngle( result->finalAngle() );
250 return result;
251 }
252
253
254 std::vector<Point>
pointAndDerivatives(Coord t,unsigned int n) const255 EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
256 {
257 if (isChord()) {
258 return chord().pointAndDerivatives(t, n);
259 }
260
261 unsigned int nn = n+1; // nn represents the size of the result vector.
262 std::vector<Point> result;
263 result.reserve(nn);
264 double angle = angleAt(t);
265 std::unique_ptr<EllipticalArc> ea( static_cast<EllipticalArc*>(duplicate()) );
266 ea->_ellipse.setCenter(0, 0);
267 unsigned int m = std::min(nn, 4u);
268 for ( unsigned int i = 0; i < m; ++i )
269 {
270 result.push_back( ea->pointAtAngle(angle) );
271 angle += (sweep() ? M_PI/2 : -M_PI/2);
272 if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
273 }
274 m = nn / 4;
275 for ( unsigned int i = 1; i < m; ++i )
276 {
277 for ( unsigned int j = 0; j < 4; ++j )
278 result.push_back( result[j] );
279 }
280 m = nn - 4 * m;
281 for ( unsigned int i = 0; i < m; ++i )
282 {
283 result.push_back( result[i] );
284 }
285 if ( !result.empty() ) // nn != 0
286 result[0] = pointAtAngle(angle);
287 return result;
288 }
289
pointAt(Coord t) const290 Point EllipticalArc::pointAt(Coord t) const
291 {
292 if (isChord()) return chord().pointAt(t);
293 return _ellipse.pointAt(angleAt(t));
294 }
295
valueAt(Coord t,Dim2 d) const296 Coord EllipticalArc::valueAt(Coord t, Dim2 d) const
297 {
298 if (isChord()) return chord().valueAt(t, d);
299 return valueAtAngle(angleAt(t), d);
300 }
301
portion(double f,double t) const302 Curve* EllipticalArc::portion(double f, double t) const
303 {
304 // fix input arguments
305 if (f < 0) f = 0;
306 if (f > 1) f = 1;
307 if (t < 0) t = 0;
308 if (t > 1) t = 1;
309
310 if (f == t) {
311 EllipticalArc *arc = new EllipticalArc();
312 arc->_initial_point = arc->_final_point = pointAt(f);
313 return arc;
314 }
315
316 EllipticalArc *arc = static_cast<EllipticalArc*>(duplicate());
317 arc->_initial_point = pointAt(f);
318 arc->_final_point = pointAt(t);
319 arc->_angles.setAngles(angleAt(f), angleAt(t));
320 if (f > t) arc->_angles.setSweep(!sweep());
321 if ( _large_arc && fabs(angularExtent() * (t-f)) <= M_PI) {
322 arc->_large_arc = false;
323 }
324 return arc;
325 }
326
327 // the arc is the same but traversed in the opposite direction
reverse() const328 Curve *EllipticalArc::reverse() const
329 {
330 using std::swap;
331 EllipticalArc *rarc = static_cast<EllipticalArc*>(duplicate());
332 rarc->_angles.reverse();
333 swap(rarc->_initial_point, rarc->_final_point);
334 return rarc;
335 }
336
337 #ifdef HAVE_GSL // GSL is required for function "solve_reals"
allNearestTimes(Point const & p,double from,double to) const338 std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from, double to ) const
339 {
340 std::vector<double> result;
341
342 if ( from > to ) std::swap(from, to);
343 if ( from < 0 || to > 1 )
344 {
345 THROW_RANGEERROR("[from,to] interval out of range");
346 }
347
348 if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )
349 {
350 result.push_back(from);
351 return result;
352 }
353 else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
354 {
355 LineSegment seg(pointAt(from), pointAt(to));
356 Point np = seg.pointAt( seg.nearestTime(p) );
357 if ( are_near(ray(Y), 0) )
358 {
359 if ( are_near(rotationAngle(), M_PI/2)
360 || are_near(rotationAngle(), 3*M_PI/2) )
361 {
362 result = roots(np[Y], Y);
363 }
364 else
365 {
366 result = roots(np[X], X);
367 }
368 }
369 else
370 {
371 if ( are_near(rotationAngle(), M_PI/2)
372 || are_near(rotationAngle(), 3*M_PI/2) )
373 {
374 result = roots(np[X], X);
375 }
376 else
377 {
378 result = roots(np[Y], Y);
379 }
380 }
381 return result;
382 }
383 else if ( are_near(ray(X), ray(Y)) )
384 {
385 Point r = p - center();
386 if ( are_near(r, Point(0,0)) )
387 {
388 THROW_INFINITESOLUTIONS(0);
389 }
390 // TODO: implement case r != 0
391 // Point np = ray(X) * unit_vector(r);
392 // std::vector<double> solX = roots(np[X],X);
393 // std::vector<double> solY = roots(np[Y],Y);
394 // double t;
395 // if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
396 // {
397 // t = solX[0];
398 // }
399 // else
400 // {
401 // t = solX[1];
402 // }
403 // if ( !(t < from || t > to) )
404 // {
405 // result.push_back(t);
406 // }
407 // else
408 // {
409 //
410 // }
411 }
412
413 // solve the equation <D(E(t),t)|E(t)-p> == 0
414 // that provides min and max distance points
415 // on the ellipse E wrt the point p
416 // after the substitutions:
417 // cos(t) = (1 - s^2) / (1 + s^2)
418 // sin(t) = 2t / (1 + s^2)
419 // where s = tan(t/2)
420 // we get a 4th degree equation in s
421 /*
422 * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
423 * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
424 * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
425 * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
426 */
427
428 Point p_c = p - center();
429 double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
430 double sinrot, cosrot;
431 sincos(rotationAngle(), sinrot, cosrot);
432 double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
433 Poly coeff;
434 coeff.resize(5);
435 coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
436 coeff[3] = 2 * ( rx2_ry2 + expr1 );
437 coeff[2] = 0;
438 coeff[1] = 2 * ( -rx2_ry2 + expr1 );
439 coeff[0] = -coeff[4];
440
441 // for ( unsigned int i = 0; i < 5; ++i )
442 // std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
443
444 std::vector<double> real_sol;
445 // gsl_poly_complex_solve raises an error
446 // if the leading coefficient is zero
447 if ( are_near(coeff[4], 0) )
448 {
449 real_sol.push_back(0);
450 if ( !are_near(coeff[3], 0) )
451 {
452 double sq = -coeff[1] / coeff[3];
453 if ( sq > 0 )
454 {
455 double s = std::sqrt(sq);
456 real_sol.push_back(s);
457 real_sol.push_back(-s);
458 }
459 }
460 }
461 else
462 {
463 real_sol = solve_reals(coeff);
464 }
465
466 for (double & i : real_sol)
467 {
468 i = 2 * std::atan(i);
469 if ( i < 0 ) i += 2*M_PI;
470 }
471 // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
472 // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
473 if ( (real_sol.size() % 2) != 0 )
474 {
475 real_sol.push_back(M_PI);
476 }
477
478 double mindistsq1 = std::numeric_limits<double>::max();
479 double mindistsq2 = std::numeric_limits<double>::max();
480 double dsq = 0;
481 unsigned int mi1 = 0, mi2 = 0;
482 for ( unsigned int i = 0; i < real_sol.size(); ++i )
483 {
484 dsq = distanceSq(p, pointAtAngle(real_sol[i]));
485 if ( mindistsq1 > dsq )
486 {
487 mindistsq2 = mindistsq1;
488 mi2 = mi1;
489 mindistsq1 = dsq;
490 mi1 = i;
491 }
492 else if ( mindistsq2 > dsq )
493 {
494 mindistsq2 = dsq;
495 mi2 = i;
496 }
497 }
498
499 double t = timeAtAngle(real_sol[mi1]);
500 if ( !(t < from || t > to) )
501 {
502 result.push_back(t);
503 }
504
505 bool second_sol = false;
506 t = timeAtAngle(real_sol[mi2]);
507 if ( real_sol.size() == 4 && !(t < from || t > to) )
508 {
509 if ( result.empty() || are_near(mindistsq1, mindistsq2) )
510 {
511 result.push_back(t);
512 second_sol = true;
513 }
514 }
515
516 // we need to test extreme points too
517 double dsq1 = distanceSq(p, pointAt(from));
518 double dsq2 = distanceSq(p, pointAt(to));
519 if ( second_sol )
520 {
521 if ( mindistsq2 > dsq1 )
522 {
523 result.clear();
524 result.push_back(from);
525 mindistsq2 = dsq1;
526 }
527 else if ( are_near(mindistsq2, dsq) )
528 {
529 result.push_back(from);
530 }
531 if ( mindistsq2 > dsq2 )
532 {
533 result.clear();
534 result.push_back(to);
535 }
536 else if ( are_near(mindistsq2, dsq2) )
537 {
538 result.push_back(to);
539 }
540
541 }
542 else
543 {
544 if ( result.empty() )
545 {
546 if ( are_near(dsq1, dsq2) )
547 {
548 result.push_back(from);
549 result.push_back(to);
550 }
551 else if ( dsq2 > dsq1 )
552 {
553 result.push_back(from);
554 }
555 else
556 {
557 result.push_back(to);
558 }
559 }
560 }
561
562 return result;
563 }
564 #endif
565
566
_filterIntersections(std::vector<ShapeIntersection> & xs,bool is_first) const567 void EllipticalArc::_filterIntersections(std::vector<ShapeIntersection> &xs, bool is_first) const
568 {
569 Interval unit(0, 1);
570 std::vector<ShapeIntersection>::reverse_iterator i = xs.rbegin(), last = xs.rend();
571 while (i != last) {
572 Coord &t = is_first ? i->first : i->second;
573 constexpr auto eps = 1e-4;
574 assert(are_near_rel(_ellipse.pointAt(t), i->point(), eps));
575 t = timeAtAngle(t);
576 if (!unit.contains(t)) {
577 xs.erase((++i).base());
578 continue;
579 } else {
580 assert(are_near_rel(pointAt(t), i->point(), eps));
581 ++i;
582 }
583 }
584 }
585
intersect(Curve const & other,Coord eps) const586 std::vector<CurveIntersection> EllipticalArc::intersect(Curve const &other, Coord eps) const
587 {
588 if (isLineSegment()) {
589 LineSegment ls(_initial_point, _final_point);
590 return ls.intersect(other, eps);
591 }
592
593 std::vector<CurveIntersection> result;
594
595 if (other.isLineSegment()) {
596 LineSegment ls(other.initialPoint(), other.finalPoint());
597 result = _ellipse.intersect(ls);
598 _filterIntersections(result, true);
599 return result;
600 }
601
602 BezierCurve const *bez = dynamic_cast<BezierCurve const *>(&other);
603 if (bez) {
604 result = _ellipse.intersect(bez->fragment());
605 _filterIntersections(result, true);
606 return result;
607 }
608
609 EllipticalArc const *arc = dynamic_cast<EllipticalArc const *>(&other);
610 if (arc) {
611 result = _ellipse.intersect(arc->_ellipse);
612 _filterIntersections(result, true);
613 arc->_filterIntersections(result, false);
614 return result;
615 }
616
617 // in case someone wants to make a custom curve type
618 result = other.intersect(*this, eps);
619 transpose_in_place(result);
620 return result;
621 }
622
623
_updateCenterAndAngles()624 void EllipticalArc::_updateCenterAndAngles()
625 {
626 // See: http://www.w3.org/TR/SVG/implnote.html#ArcImplementationNotes
627 Point d = initialPoint() - finalPoint();
628 Point mid = middle_point(initialPoint(), finalPoint());
629
630 // if ip = sp, the arc contains no other points
631 if (initialPoint() == finalPoint()) {
632 _ellipse = Ellipse();
633 _ellipse.setCenter(initialPoint());
634 _angles = AngleInterval();
635 _large_arc = false;
636 return;
637 }
638
639 // rays should be positive
640 _ellipse.setRays(std::fabs(ray(X)), std::fabs(ray(Y)));
641
642 if (isChord()) {
643 _ellipse.setRays(L2(d) / 2, 0);
644 _ellipse.setRotationAngle(atan2(d));
645 _ellipse.setCenter(mid);
646 _angles.setAngles(0, M_PI);
647 _angles.setSweep(false);
648 _large_arc = false;
649 return;
650 }
651
652 Rotate rot(rotationAngle()); // the matrix in F.6.5.3
653 Rotate invrot = rot.inverse(); // the matrix in F.6.5.1
654
655 Point r = rays();
656 Point p = (initialPoint() - mid) * invrot; // x', y' in F.6.5.1
657 Point c(0,0); // cx', cy' in F.6.5.2
658
659 // Correct out-of-range radii
660 Coord lambda = hypot(p[X]/r[X], p[Y]/r[Y]);
661 if (lambda > 1) {
662 r *= lambda;
663 _ellipse.setRays(r);
664 _ellipse.setCenter(mid);
665 } else {
666 // evaluate F.6.5.2
667 Coord rxry = r[X]*r[X] * r[Y]*r[Y];
668 Coord pxry = p[X]*p[X] * r[Y]*r[Y];
669 Coord rxpy = r[X]*r[X] * p[Y]*p[Y];
670 Coord rad = (rxry - pxry - rxpy)/(rxpy + pxry);
671 // normally rad should never be negative, but numerical inaccuracy may cause this
672 if (rad > 0) {
673 rad = std::sqrt(rad);
674 if (sweep() == _large_arc) {
675 rad = -rad;
676 }
677 c = rad * Point(r[X]*p[Y]/r[Y], -r[Y]*p[X]/r[X]);
678 _ellipse.setCenter(c * rot + mid);
679 } else {
680 _ellipse.setCenter(mid);
681 }
682 }
683
684 // Compute start and end angles.
685 // If the ellipse was enlarged, c will be zero - this is correct.
686 Point sp((p[X] - c[X]) / r[X], (p[Y] - c[Y]) / r[Y]);
687 Point ep((-p[X] - c[X]) / r[X], (-p[Y] - c[Y]) / r[Y]);
688 Point v(1, 0);
689
690 _angles.setInitial(angle_between(v, sp));
691 _angles.setFinal(angle_between(v, ep));
692
693 /*double sweep_angle = angle_between(sp, ep);
694 if (!sweep() && sweep_angle > 0) sweep_angle -= 2*M_PI;
695 if (sweep() && sweep_angle < 0) sweep_angle += 2*M_PI;*/
696 }
697
toSBasis() const698 D2<SBasis> EllipticalArc::toSBasis() const
699 {
700 if (isChord()) {
701 return chord().toSBasis();
702 }
703
704 D2<SBasis> arc;
705 // the interval of parametrization has to be [0,1]
706 Coord et = initialAngle().radians() + sweepAngle();
707 Linear param(initialAngle().radians(), et);
708 Coord cosrot, sinrot;
709 sincos(rotationAngle(), sinrot, cosrot);
710
711 // order = 4 seems to be enough to get a perfect looking elliptical arc
712 SBasis arc_x = ray(X) * cos(param,4);
713 SBasis arc_y = ray(Y) * sin(param,4);
714 arc[0] = arc_x * cosrot - arc_y * sinrot + Linear(center(X), center(X));
715 arc[1] = arc_x * sinrot + arc_y * cosrot + Linear(center(Y), center(Y));
716
717 // ensure that endpoints remain exact
718 for ( int d = 0 ; d < 2 ; d++ ) {
719 arc[d][0][0] = initialPoint()[d];
720 arc[d][0][1] = finalPoint()[d];
721 }
722
723 return arc;
724 }
725
726 // All operations that do not contain skew can be evaluated
727 // without passing through the implicit form of the ellipse,
728 // which preserves precision.
729
operator *=(Translate const & tr)730 void EllipticalArc::operator*=(Translate const &tr)
731 {
732 _initial_point *= tr;
733 _final_point *= tr;
734 _ellipse *= tr;
735 }
736
operator *=(Scale const & s)737 void EllipticalArc::operator*=(Scale const &s)
738 {
739 _initial_point *= s;
740 _final_point *= s;
741 _ellipse *= s;
742 }
743
operator *=(Rotate const & r)744 void EllipticalArc::operator*=(Rotate const &r)
745 {
746 _initial_point *= r;
747 _final_point *= r;
748 _ellipse *= r;
749 }
750
operator *=(Zoom const & z)751 void EllipticalArc::operator*=(Zoom const &z)
752 {
753 _initial_point *= z;
754 _final_point *= z;
755 _ellipse *= z;
756 }
757
operator *=(Affine const & m)758 void EllipticalArc::operator*=(Affine const& m)
759 {
760 if (isChord()) {
761 _initial_point *= m;
762 _final_point *= m;
763 _ellipse.setCenter(middle_point(_initial_point, _final_point));
764 _ellipse.setRays(0, 0);
765 _ellipse.setRotationAngle(0);
766 return;
767 }
768
769 _initial_point *= m;
770 _final_point *= m;
771 _ellipse *= m;
772 if (m.det() < 0) {
773 _angles.setSweep(!sweep());
774 }
775
776 // ellipse transformation does not preserve its functional form,
777 // i.e. e.pointAt(0.5)*m and (e*m).pointAt(0.5) can be different.
778 // We need to recompute start / end angles.
779 _angles.setInitial(_ellipse.timeAt(_initial_point));
780 _angles.setFinal(_ellipse.timeAt(_final_point));
781 }
782
operator ==(Curve const & c) const783 bool EllipticalArc::operator==(Curve const &c) const
784 {
785 EllipticalArc const *other = dynamic_cast<EllipticalArc const *>(&c);
786 if (!other) return false;
787 if (_initial_point != other->_initial_point) return false;
788 if (_final_point != other->_final_point) return false;
789 // TODO: all arcs with ellipse rays which are too small
790 // and fall back to a line should probably be equal
791 if (rays() != other->rays()) return false;
792 if (rotationAngle() != other->rotationAngle()) return false;
793 if (_large_arc != other->_large_arc) return false;
794 if (sweep() != other->sweep()) return false;
795 return true;
796 }
797
isNear(Curve const & c,Coord precision) const798 bool EllipticalArc::isNear(Curve const &c, Coord precision) const
799 {
800 EllipticalArc const *other = dynamic_cast<EllipticalArc const *>(&c);
801 if (!other) {
802 if (isChord()) {
803 return c.isNear(chord(), precision);
804 }
805 return false;
806 }
807
808 if (!are_near(_initial_point, other->_initial_point, precision)) return false;
809 if (!are_near(_final_point, other->_final_point, precision)) return false;
810 if (isChord() && other->isChord()) return true;
811
812 if (sweep() != other->sweep()) return false;
813 if (!are_near(_ellipse, other->_ellipse, precision)) return false;
814 return true;
815 }
816
feed(PathSink & sink,bool moveto_initial) const817 void EllipticalArc::feed(PathSink &sink, bool moveto_initial) const
818 {
819 if (moveto_initial) {
820 sink.moveTo(_initial_point);
821 }
822 sink.arcTo(ray(X), ray(Y), rotationAngle(), _large_arc, sweep(), _final_point);
823 }
824
winding(Point const & p) const825 int EllipticalArc::winding(Point const &p) const
826 {
827 using std::swap;
828
829 double sinrot, cosrot;
830 sincos(rotationAngle(), sinrot, cosrot);
831
832 Angle ymin_a = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
833 Angle ymax_a = ymin_a + M_PI;
834
835 Point ymin = pointAtAngle(ymin_a);
836 Point ymax = pointAtAngle(ymax_a);
837 if (ymin[Y] > ymax[Y]) {
838 swap(ymin, ymax);
839 swap(ymin_a, ymax_a);
840 }
841
842 Interval yspan(ymin[Y], ymax[Y]);
843 if (!yspan.lowerContains(p[Y])) return 0;
844
845 bool left = cross(ymax - ymin, p - ymin) > 0;
846 bool inside = _ellipse.contains(p);
847 bool includes_ymin = _angles.contains(ymin_a);
848 bool includes_ymax = _angles.contains(ymax_a);
849
850 AngleInterval rarc(ymin_a, ymax_a, true),
851 larc(ymax_a, ymin_a, true);
852
853 // we'll compute the result for an arc in the direction of increasing angles
854 // and then negate if necessary
855 Angle ia = initialAngle(), fa = finalAngle();
856 Point ip = _initial_point, fp = _final_point;
857 if (!sweep()) {
858 swap(ia, fa);
859 swap(ip, fp);
860 }
861
862 bool initial_left = larc.contains(ia);
863 bool initial_right = !initial_left; // rarc.contains(ia);
864 bool final_left = larc.contains(fa);
865 bool final_right = !final_left; // rarc.contains(fa);
866
867 int result = 0;
868 if (inside || left) {
869 if (includes_ymin && final_right) {
870 Interval ival(ymin[Y], fp[Y]);
871 if (ival.lowerContains(p[Y])) {
872 ++result;
873 }
874 }
875 if (initial_right && final_right && !largeArc()) {
876 Interval ival(ip[Y], fp[Y]);
877 if (ival.lowerContains(p[Y])) {
878 ++result;
879 }
880 }
881 if (initial_right && includes_ymax) {
882 Interval ival(ip[Y], ymax[Y]);
883 if (ival.lowerContains(p[Y])) {
884 ++result;
885 }
886 }
887 if (!initial_right && !final_right && includes_ymin && includes_ymax) {
888 Interval ival(ymax[Y], ymin[Y]);
889 if (ival.lowerContains(p[Y])) {
890 ++result;
891 }
892 }
893 }
894 if (left && !inside) {
895 if (includes_ymin && initial_left) {
896 Interval ival(ymin[Y], ip[Y]);
897 if (ival.lowerContains(p[Y])) {
898 --result;
899 }
900 }
901 if (initial_left && final_left && !largeArc()) {
902 Interval ival(ip[Y], fp[Y]);
903 if (ival.lowerContains(p[Y])) {
904 --result;
905 }
906 }
907 if (final_left && includes_ymax) {
908 Interval ival(fp[Y], ymax[Y]);
909 if (ival.lowerContains(p[Y])) {
910 --result;
911 }
912 }
913 if (!initial_left && !final_left && includes_ymin && includes_ymax) {
914 Interval ival(ymax[Y], ymin[Y]);
915 if (ival.lowerContains(p[Y])) {
916 --result;
917 }
918 }
919 }
920 return sweep() ? result : -result;
921 }
922
operator <<(std::ostream & out,EllipticalArc const & ea)923 std::ostream &operator<<(std::ostream &out, EllipticalArc const &ea)
924 {
925 out << "EllipticalArc("
926 << ea.initialPoint() << ", "
927 << format_coord_nice(ea.ray(X)) << ", " << format_coord_nice(ea.ray(Y)) << ", "
928 << format_coord_nice(ea.rotationAngle()) << ", "
929 << "large_arc=" << (ea.largeArc() ? "true" : "false") << ", "
930 << "sweep=" << (ea.sweep() ? "true" : "false") << ", "
931 << ea.finalPoint() << ")";
932 return out;
933 }
934
935 } // end namespace Geom
936
937 /*
938 Local Variables:
939 mode:c++
940 c-file-style:"stroustrup"
941 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
942 indent-tabs-mode:nil
943 fill-column:99
944 End:
945 */
946 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
947
948