1 /*
2 * SVG Elliptical Arc Class
3 *
4 * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it either under the terms of the GNU Lesser General Public
8 * License version 2.1 as published by the Free Software Foundation
9 * (the "LGPL") or, at your option, under the terms of the Mozilla
10 * Public License Version 1.1 (the "MPL"). If you do not alter this
11 * notice, a recipient may use your version of this file under either
12 * the MPL or the LGPL.
13 *
14 * You should have received a copy of the LGPL along with this library
15 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 * You should have received a copy of the MPL along with this library
18 * in the file COPYING-MPL-1.1
19 *
20 * The contents of this file are subject to the Mozilla Public License
21 * Version 1.1 (the "License"); you may not use this file except in
22 * compliance with the License. You may obtain a copy of the License at
23 * http://www.mozilla.org/MPL/
24 *
25 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
26 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
27 * the specific language governing rights and limitations.
28 */
31 #include <2geom/elliptical-arc.h>
32 #include <2geom/bezier-curve.h>
33 #include <2geom/poly.h>
34 #include <2geom/sbasis-math.h>
36 #include <cfloat>
37 #include <limits>
42 namespace Geom
43 {
46 OptRect EllipticalArc::boundsExact() const
47 {
48 std::vector<double> extremes(4);
49 double cosrot = std::cos(rotation_angle());
50 double sinrot = std::sin(rotation_angle());
51 extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
52 extremes[1] = extremes[0] + M_PI;
53 if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;
54 extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
55 extremes[3] = extremes[2] + M_PI;
56 if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;
59 std::vector<double>arc_extremes(4);
60 arc_extremes[0] = initialPoint()[X];
61 arc_extremes[1] = finalPoint()[X];
62 if ( arc_extremes[0] < arc_extremes[1] )
63 std::swap(arc_extremes[0], arc_extremes[1]);
64 arc_extremes[2] = initialPoint()[Y];
65 arc_extremes[3] = finalPoint()[Y];
66 if ( arc_extremes[2] < arc_extremes[3] )
67 std::swap(arc_extremes[2], arc_extremes[3]);
70 if ( start_angle() < end_angle() )
71 {
72 if ( sweep_flag() )
73 {
74 for ( unsigned int i = 0; i < extremes.size(); ++i )
75 {
76 if ( start_angle() < extremes[i] && extremes[i] < end_angle() )
77 {
78 arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
79 }
80 }
81 }
82 else
83 {
84 for ( unsigned int i = 0; i < extremes.size(); ++i )
85 {
86 if ( start_angle() > extremes[i] || extremes[i] > end_angle() )
87 {
88 arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
89 }
90 }
91 }
92 }
93 else
94 {
95 if ( sweep_flag() )
96 {
97 for ( unsigned int i = 0; i < extremes.size(); ++i )
98 {
99 if ( start_angle() < extremes[i] || extremes[i] < end_angle() )
100 {
101 arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
102 }
103 }
104 }
105 else
106 {
107 for ( unsigned int i = 0; i < extremes.size(); ++i )
108 {
109 if ( start_angle() > extremes[i] && extremes[i] > end_angle() )
110 {
111 arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
112 }
113 }
114 }
115 }
117 return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
118 Point(arc_extremes[0], arc_extremes[2]) );
120 }
123 std::vector<double>
124 EllipticalArc::roots(double v, Dim2 d) const
125 {
126 if ( d > Y )
127 {
128 THROW_RANGEERROR("dimention out of range");
129 }
131 std::vector<double> sol;
132 if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
133 {
134 if ( center(d) == v )
135 sol.push_back(0);
136 return sol;
137 }
139 const char* msg[2][2] =
140 {
141 { "d == X; ray(X) == 0; "
142 "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "
143 "s should be contained in [-1,1]",
144 "d == X; ray(Y) == 0; "
145 "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "
146 "s should be contained in [-1,1]"
147 },
148 { "d == Y; ray(X) == 0; "
149 "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "
150 "s should be contained in [-1,1]",
151 "d == Y; ray(Y) == 0; "
152 "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "
153 "s should be contained in [-1,1]"
154 },
155 };
157 for ( unsigned int dim = 0; dim < 2; ++dim )
158 {
159 if ( are_near(ray(dim), 0) )
160 {
162 if ( initialPoint()[d] == v && finalPoint()[d] == v )
163 {
164 THROW_INFINITESOLUTIONS(0);
165 }
166 if ( (initialPoint()[d] < finalPoint()[d])
167 && (initialPoint()[d] > v || finalPoint()[d] < v) )
168 {
169 return sol;
170 }
171 if ( (initialPoint()[d] > finalPoint()[d])
172 && (finalPoint()[d] > v || initialPoint()[d] < v) )
173 {
174 return sol;
175 }
176 double ray_prj;
177 switch(d)
178 {
179 case X:
180 switch(dim)
181 {
182 case X: ray_prj = -ray(Y) * std::sin(rotation_angle());
183 break;
184 case Y: ray_prj = ray(X) * std::cos(rotation_angle());
185 break;
186 }
187 break;
188 case Y:
189 switch(dim)
190 {
191 case X: ray_prj = ray(Y) * std::cos(rotation_angle());
192 break;
193 case Y: ray_prj = ray(X) * std::sin(rotation_angle());
194 break;
195 }
196 break;
197 }
199 double s = (v - center(d)) / ray_prj;
200 if ( s < -1 || s > 1 )
201 {
202 THROW_LOGICALERROR(msg[d][dim]);
203 }
204 switch(dim)
205 {
206 case X:
207 s = std::asin(s); // return a value in [-PI/2,PI/2]
208 if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) ) )
209 {
210 if ( s < 0 ) s += 2*M_PI;
211 }
212 else
213 {
214 s = M_PI - s;
215 if (!(s < 2*M_PI) ) s -= 2*M_PI;
216 }
217 break;
218 case Y:
219 s = std::acos(s); // return a value in [0,PI]
220 if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )
221 {
222 s = 2*M_PI - s;
223 if ( !(s < 2*M_PI) ) s -= 2*M_PI;
224 }
225 break;
226 }
228 //std::cerr << "s = " << rad_to_deg(s);
229 s = map_to_01(s);
230 //std::cerr << " -> t: " << s << std::endl;
231 if ( !(s < 0 || s > 1) )
232 sol.push_back(s);
233 return sol;
234 }
235 }
237 double rotx, roty;
238 switch(d)
239 {
240 case X:
241 rotx = std::cos(rotation_angle());
242 roty = -std::sin(rotation_angle());
243 break;
244 case Y:
245 rotx = std::sin(rotation_angle());
246 roty = std::cos(rotation_angle());
247 break;
248 }
249 double rxrotx = ray(X) * rotx;
250 double c_v = center(d) - v;
252 double a = -rxrotx + c_v;
253 double b = ray(Y) * roty;
254 double c = rxrotx + c_v;
255 //std::cerr << "a = " << a << std::endl;
256 //std::cerr << "b = " << b << std::endl;
257 //std::cerr << "c = " << c << std::endl;
259 if ( are_near(a,0) )
260 {
261 sol.push_back(M_PI);
262 if ( !are_near(b,0) )
263 {
264 double s = 2 * std::atan(-c/(2*b));
265 if ( s < 0 ) s += 2*M_PI;
266 sol.push_back(s);
267 }
268 }
269 else
270 {
271 double delta = b * b - a * c;
272 //std::cerr << "delta = " << delta << std::endl;
273 if ( are_near(delta, 0) )
274 {
275 double s = 2 * std::atan(-b/a);
276 if ( s < 0 ) s += 2*M_PI;
277 sol.push_back(s);
278 }
279 else if ( delta > 0 )
280 {
281 double sq = std::sqrt(delta);
282 double s = 2 * std::atan( (-b - sq) / a );
283 if ( s < 0 ) s += 2*M_PI;
284 sol.push_back(s);
285 s = 2 * std::atan( (-b + sq) / a );
286 if ( s < 0 ) s += 2*M_PI;
287 sol.push_back(s);
288 }
289 }
291 std::vector<double> arc_sol;
292 for (unsigned int i = 0; i < sol.size(); ++i )
293 {
294 //std::cerr << "s = " << rad_to_deg(sol[i]);
295 sol[i] = map_to_01(sol[i]);
296 //std::cerr << " -> t: " << sol[i] << std::endl;
297 if ( !(sol[i] < 0 || sol[i] > 1) )
298 arc_sol.push_back(sol[i]);
299 }
300 return arc_sol;
303 // return SBasisCurve(toSBasis()).roots(v, d);
304 }
306 // D(E(t,C),t) = E(t+PI/2,O)
307 Curve* EllipticalArc::derivative() const
308 {
309 EllipticalArc* result = new EllipticalArc(*this);
310 result->m_center[X] = result->m_center[Y] = 0;
311 result->m_start_angle += M_PI/2;
312 if( !( result->m_start_angle < 2*M_PI ) )
313 {
314 result->m_start_angle -= 2*M_PI;
315 }
316 result->m_end_angle += M_PI/2;
317 if( !( result->m_end_angle < 2*M_PI ) )
318 {
319 result->m_end_angle -= 2*M_PI;
320 }
321 result->m_initial_point = result->pointAtAngle( result->start_angle() );
322 result->m_final_point = result->pointAtAngle( result->end_angle() );
323 return result;
325 }
327 std::vector<Point>
328 EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
329 {
330 unsigned int nn = n+1; // nn represents the size of the result vector.
331 std::vector<Point> result;
332 result.reserve(nn);
333 double angle = map_unit_interval_on_circular_arc(t, start_angle(),
334 end_angle(), sweep_flag());
335 EllipticalArc ea(*this);
336 ea.m_center = Point(0,0);
337 unsigned int m = std::min(nn, 4u);
338 for ( unsigned int i = 0; i < m; ++i )
339 {
340 result.push_back( ea.pointAtAngle(angle) );
341 angle += M_PI/2;
342 if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
343 }
344 m = nn / 4;
345 for ( unsigned int i = 1; i < m; ++i )
346 {
347 for ( unsigned int j = 0; j < 4; ++j )
348 result.push_back( result[j] );
349 }
350 m = nn - 4 * m;
351 for ( unsigned int i = 0; i < m; ++i )
352 {
353 result.push_back( result[i] );
354 }
355 if ( !result.empty() ) // nn != 0
356 result[0] = pointAtAngle(angle);
357 return result;
358 }
360 D2<SBasis> EllipticalArc::toSBasis() const
361 {
362 // the interval of parametrization has to be [0,1]
363 Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );
364 Linear param(start_angle(), et);
365 Coord cos_rot_angle = std::cos(rotation_angle());
366 Coord sin_rot_angle = std::sin(rotation_angle());
367 // order = 4 seems to be enough to get a perfect looking elliptical arc
368 // should it be choosen in function of the arc length anyway ?
369 // or maybe a user settable parameter: toSBasis(unsigned int order) ?
370 SBasis arc_x = ray(X) * cos(param,4);
371 SBasis arc_y = ray(Y) * sin(param,4);
372 D2<SBasis> arc;
373 arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
374 arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));
375 return arc;
376 }
379 bool EllipticalArc::containsAngle(Coord angle) const
380 {
381 if ( sweep_flag() )
382 if ( start_angle() < end_angle() )
383 return ( !( angle < start_angle() || angle > end_angle() ) );
384 else
385 return ( !( angle < start_angle() && angle > end_angle() ) );
386 else
387 if ( start_angle() > end_angle() )
388 return ( !( angle > start_angle() || angle < end_angle() ) );
389 else
390 return ( !( angle > start_angle() && angle < end_angle() ) );
391 }
394 double EllipticalArc::valueAtAngle(Coord t, Dim2 d) const
395 {
396 double sin_rot_angle = std::sin(rotation_angle());
397 double cos_rot_angle = std::cos(rotation_angle());
398 if ( d == X )
399 {
400 return ray(X) * cos_rot_angle * std::cos(t)
401 - ray(Y) * sin_rot_angle * std::sin(t)
402 + center(X);
403 }
404 else if ( d == Y )
405 {
406 return ray(X) * sin_rot_angle * std::cos(t)
407 + ray(Y) * cos_rot_angle * std::sin(t)
408 + center(Y);
409 }
410 THROW_RANGEERROR("dimension parameter out of range");
411 }
414 Curve* EllipticalArc::portion(double f, double t) const
415 {
416 if (f < 0) f = 0;
417 if (f > 1) f = 1;
418 if (t < 0) t = 0;
419 if (t > 1) t = 1;
420 if ( are_near(f, t) )
421 {
422 EllipticalArc* arc = new EllipticalArc();
423 arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);
424 arc->m_start_angle = arc->m_end_angle = m_start_angle;
425 arc->m_rot_angle = m_rot_angle;
426 arc->m_sweep = m_sweep;
427 arc->m_large_arc = m_large_arc;
428 }
429 EllipticalArc* arc = new EllipticalArc( *this );
430 arc->m_initial_point = pointAt(f);
431 arc->m_final_point = pointAt(t);
432 double sa = sweep_flag() ? sweep_angle() : -sweep_angle();
433 arc->m_start_angle = m_start_angle + sa * f;
434 if ( !(arc->m_start_angle < 2*M_PI) )
435 arc->m_start_angle -= 2*M_PI;
436 if ( arc->m_start_angle < 0 )
437 arc->m_start_angle += 2*M_PI;
438 arc->m_end_angle = m_start_angle + sa * t;
439 if ( !(arc->m_end_angle < 2*M_PI) )
440 arc->m_end_angle -= 2*M_PI;
441 if ( arc->m_end_angle < 0 )
442 arc->m_end_angle += 2*M_PI;
443 if ( f > t ) arc->m_sweep = !sweep_flag();
444 if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )
445 arc->m_large_arc = false;
446 return arc;
447 }
449 // NOTE: doesn't work with 360 deg arcs
450 void EllipticalArc::calculate_center_and_extreme_angles()
451 {
452 // as stated in the svg standard the rotation angle parameter
453 // value must be modulo 2*PI
454 m_rot_angle = std::fmod(m_rot_angle, 2*M_PI);
455 if (m_rot_angle < 0) m_rot_angle += 2*M_PI;
457 if ( are_near(initialPoint(), finalPoint()) )
458 {
459 if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
460 {
461 m_start_angle = m_end_angle = 0;
462 m_center = initialPoint();
463 return;
464 }
465 else
466 {
467 THROW_RANGEERROR("initial and final point are the same");
468 }
469 }
470 if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
471 { // but initialPoint != finalPoint
472 THROW_RANGEERROR(
473 "there is no ellipse that satisfies the given constraints: "
474 "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"
475 );
476 }
477 if ( are_near(ray(Y), 0) )
478 {
479 Point v = initialPoint() - finalPoint();
480 if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )
481 {
482 double angle = std::atan2(v[Y], v[X]);
483 if (angle < 0) angle += 2*M_PI;
484 if ( are_near( angle, rotation_angle() ) )
485 {
486 m_start_angle = 0;
487 m_end_angle = M_PI;
488 m_center = v/2 + finalPoint();
489 return;
490 }
491 angle -= M_PI;
492 if ( angle < 0 ) angle += 2*M_PI;
493 if ( are_near( angle, rotation_angle() ) )
494 {
495 m_start_angle = M_PI;
496 m_end_angle = 0;
497 m_center = v/2 + finalPoint();
498 return;
499 }
500 THROW_RANGEERROR(
501 "there is no ellipse that satisfies the given constraints: "
502 "ray(Y) == 0 "
503 "and slope(initialPoint - finalPoint) != rotation_angle "
504 "and != rotation_angle + PI"
505 );
506 }
507 if ( L2sq(v) > 4*ray(X)*ray(X) )
508 {
509 THROW_RANGEERROR(
510 "there is no ellipse that satisfies the given constraints: "
511 "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"
512 );
513 }
514 else
515 {
516 THROW_RANGEERROR(
517 "there is infinite ellipses that satisfy the given constraints: "
518 "ray(Y) == 0 and distance(initialPoint, finalPoint) < 2*ray(X)"
519 );
520 }
522 }
524 if ( are_near(ray(X), 0) )
525 {
526 Point v = initialPoint() - finalPoint();
527 if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )
528 {
529 double angle = std::atan2(v[Y], v[X]);
530 if (angle < 0) angle += 2*M_PI;
531 double rot_angle = rotation_angle() + M_PI/2;
532 if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;
533 if ( are_near( angle, rot_angle ) )
534 {
535 m_start_angle = M_PI/2;
536 m_end_angle = 3*M_PI/2;
537 m_center = v/2 + finalPoint();
538 return;
539 }
540 angle -= M_PI;
541 if ( angle < 0 ) angle += 2*M_PI;
542 if ( are_near( angle, rot_angle ) )
543 {
544 m_start_angle = 3*M_PI/2;
545 m_end_angle = M_PI/2;
546 m_center = v/2 + finalPoint();
547 return;
548 }
549 THROW_RANGEERROR(
550 "there is no ellipse that satisfies the given constraints: "
551 "ray(X) == 0 "
552 "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "
553 "and != rotation_angle + (3/2)*PI"
554 );
555 }
556 if ( L2sq(v) > 4*ray(Y)*ray(Y) )
557 {
558 THROW_RANGEERROR(
559 "there is no ellipse that satisfies the given constraints: "
560 "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"
561 );
562 }
563 else
564 {
565 THROW_RANGEERROR(
566 "there is infinite ellipses that satisfy the given constraints: "
567 "ray(X) == 0 and distance(initialPoint, finalPoint) < 2*ray(Y)"
568 );
569 }
571 }
573 double sin_rot_angle = std::sin(rotation_angle());
574 double cos_rot_angle = std::cos(rotation_angle());
576 Point sp = sweep_flag() ? initialPoint() : finalPoint();
577 Point ep = sweep_flag() ? finalPoint() : initialPoint();
579 Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,
580 -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,
581 0, 0 );
582 Matrix im = m.inverse();
583 Point sol = (ep - sp) * im;
584 double half_sum_angle = std::atan2(-sol[X], sol[Y]);
585 double half_diff_angle;
586 if ( are_near(std::fabs(half_sum_angle), M_PI/2) )
587 {
588 double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1;
589 double arg = anti_sgn_hsa * sol[X] / 2;
590 // if |arg| is a little bit > 1 acos returns nan
591 if ( are_near(arg, 1) )
592 half_diff_angle = 0;
593 else if ( are_near(arg, -1) )
594 half_diff_angle = M_PI;
595 else
596 {
597 if ( !(-1 < arg && arg < 1) )
598 THROW_RANGEERROR(
599 "there is no ellipse that satisfies the given constraints"
600 );
601 // assert( -1 < arg && arg < 1 );
602 // if it fails
603 // => there is no ellipse that satisfies the given constraints
604 half_diff_angle = std::acos( arg );
605 }
607 half_diff_angle = M_PI/2 - half_diff_angle;
608 }
609 else
610 {
611 double arg = sol[Y] / ( 2 * std::cos(half_sum_angle) );
612 // if |arg| is a little bit > 1 asin returns nan
613 if ( are_near(arg, 1) )
614 half_diff_angle = M_PI/2;
615 else if ( are_near(arg, -1) )
616 half_diff_angle = -M_PI/2;
617 else
618 {
619 if ( !(-1 < arg && arg < 1) )
620 THROW_RANGEERROR(
621 "there is no ellipse that satisfies the given constraints"
622 );
623 // assert( -1 < arg && arg < 1 );
624 // if it fails
625 // => there is no ellipse that satisfies the given constraints
626 half_diff_angle = std::asin( arg );
627 }
628 }
630 if ( ( m_large_arc && half_diff_angle > 0 )
631 || (!m_large_arc && half_diff_angle < 0 ) )
632 {
633 half_diff_angle = -half_diff_angle;
634 }
635 if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI;
636 if ( half_diff_angle < 0 ) half_diff_angle += M_PI;
638 m_start_angle = half_sum_angle - half_diff_angle;
639 m_end_angle = half_sum_angle + half_diff_angle;
640 // 0 <= m_start_angle, m_end_angle < 2PI
641 if ( m_start_angle < 0 ) m_start_angle += 2*M_PI;
642 if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI;
643 sol[0] = std::cos(m_start_angle);
644 sol[1] = std::sin(m_start_angle);
645 m_center = sp - sol * m;
646 if ( !sweep_flag() )
647 {
648 double angle = m_start_angle;
649 m_start_angle = m_end_angle;
650 m_end_angle = angle;
651 }
652 }
654 Coord EllipticalArc::map_to_02PI(Coord t) const
655 {
656 if ( sweep_flag() )
657 {
658 Coord angle = start_angle() + sweep_angle() * t;
659 if ( !(angle < 2*M_PI) )
660 angle -= 2*M_PI;
661 return angle;
662 }
663 else
664 {
665 Coord angle = start_angle() - sweep_angle() * t;
666 if ( angle < 0 ) angle += 2*M_PI;
667 return angle;
668 }
669 }
671 Coord EllipticalArc::map_to_01(Coord angle) const
672 {
673 return map_circular_arc_on_unit_interval(angle, start_angle(),
674 end_angle(), sweep_flag());
675 }
678 std::vector<double> EllipticalArc::
679 allNearestPoints( Point const& p, double from, double to ) const
680 {
681 if ( from > to ) std::swap(from, to);
682 if ( from < 0 || to > 1 )
683 {
684 THROW_RANGEERROR("[from,to] interval out of range");
685 }
686 std::vector<double> result;
687 if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )
688 {
689 result.push_back(from);
690 return result;
691 }
692 else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
693 {
694 LineSegment seg(pointAt(from), pointAt(to));
695 Point np = seg.pointAt( seg.nearestPoint(p) );
696 if ( are_near(ray(Y), 0) )
697 {
698 if ( are_near(rotation_angle(), M_PI/2)
699 || are_near(rotation_angle(), 3*M_PI/2) )
700 {
701 result = roots(np[Y], Y);
702 }
703 else
704 {
705 result = roots(np[X], X);
706 }
707 }
708 else
709 {
710 if ( are_near(rotation_angle(), M_PI/2)
711 || are_near(rotation_angle(), 3*M_PI/2) )
712 {
713 result = roots(np[X], X);
714 }
715 else
716 {
717 result = roots(np[Y], Y);
718 }
719 }
720 return result;
721 }
722 else if ( are_near(ray(X), ray(Y)) )
723 {
724 Point r = p - center();
725 if ( are_near(r, Point(0,0)) )
726 {
727 THROW_INFINITESOLUTIONS(0);
728 }
729 // TODO: implement case r != 0
730 // Point np = ray(X) * unit_vector(r);
731 // std::vector<double> solX = roots(np[X],X);
732 // std::vector<double> solY = roots(np[Y],Y);
733 // double t;
734 // if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
735 // {
736 // t = solX[0];
737 // }
738 // else
739 // {
740 // t = solX[1];
741 // }
742 // if ( !(t < from || t > to) )
743 // {
744 // result.push_back(t);
745 // }
746 // else
747 // {
748 //
749 // }
750 }
752 // solve the equation <D(E(t),t)|E(t)-p> == 0
753 // that provides min and max distance points
754 // on the ellipse E wrt the point p
755 // after the substitutions:
756 // cos(t) = (1 - s^2) / (1 + s^2)
757 // sin(t) = 2t / (1 + s^2)
758 // where s = tan(t/2)
759 // we get a 4th degree equation in s
760 /*
761 * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
762 * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
763 * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
764 * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
765 */
767 Point p_c = p - center();
768 double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
769 double cosrot = std::cos( rotation_angle() );
770 double sinrot = std::sin( rotation_angle() );
771 double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
772 Poly coeff;
773 coeff.resize(5);
774 coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
775 coeff[3] = 2 * ( rx2_ry2 + expr1 );
776 coeff[2] = 0;
777 coeff[1] = 2 * ( -rx2_ry2 + expr1 );
778 coeff[0] = -coeff[4];
780 // for ( unsigned int i = 0; i < 5; ++i )
781 // std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
783 std::vector<double> real_sol;
784 // gsl_poly_complex_solve raises an error
785 // if the leading coefficient is zero
786 if ( are_near(coeff[4], 0) )
787 {
788 real_sol.push_back(0);
789 if ( !are_near(coeff[3], 0) )
790 {
791 double sq = -coeff[1] / coeff[3];
792 if ( sq > 0 )
793 {
794 double s = std::sqrt(sq);
795 real_sol.push_back(s);
796 real_sol.push_back(-s);
797 }
798 }
799 }
800 else
801 {
802 real_sol = solve_reals(coeff);
803 }
804 // else
805 // {
806 // double sol[8];
807 // gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5);
808 // gsl_poly_complex_solve(coeff, 5, w, sol );
809 // gsl_poly_complex_workspace_free(w);
810 //
811 // for ( unsigned int i = 0; i < 4; ++i )
812 // {
813 // if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]);
814 // }
815 // }
817 for ( unsigned int i = 0; i < real_sol.size(); ++i )
818 {
819 real_sol[i] = 2 * std::atan(real_sol[i]);
820 if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;
821 }
822 // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
823 // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
824 if ( (real_sol.size() % 2) != 0 )
825 {
826 real_sol.push_back(M_PI);
827 }
829 double mindistsq1 = std::numeric_limits<double>::max();
830 double mindistsq2 = std::numeric_limits<double>::max();
831 double dsq;
832 unsigned int mi1, mi2;
833 for ( unsigned int i = 0; i < real_sol.size(); ++i )
834 {
835 dsq = distanceSq(p, pointAtAngle(real_sol[i]));
836 if ( mindistsq1 > dsq )
837 {
838 mindistsq2 = mindistsq1;
839 mi2 = mi1;
840 mindistsq1 = dsq;
841 mi1 = i;
842 }
843 else if ( mindistsq2 > dsq )
844 {
845 mindistsq2 = dsq;
846 mi2 = i;
847 }
848 }
850 double t = map_to_01( real_sol[mi1] );
851 if ( !(t < from || t > to) )
852 {
853 result.push_back(t);
854 }
856 bool second_sol = false;
857 t = map_to_01( real_sol[mi2] );
858 if ( real_sol.size() == 4 && !(t < from || t > to) )
859 {
860 if ( result.empty() || are_near(mindistsq1, mindistsq2) )
861 {
862 result.push_back(t);
863 second_sol = true;
864 }
865 }
867 // we need to test extreme points too
868 double dsq1 = distanceSq(p, pointAt(from));
869 double dsq2 = distanceSq(p, pointAt(to));
870 if ( second_sol )
871 {
872 if ( mindistsq2 > dsq1 )
873 {
874 result.clear();
875 result.push_back(from);
876 mindistsq2 = dsq1;
877 }
878 else if ( are_near(mindistsq2, dsq) )
879 {
880 result.push_back(from);
881 }
882 if ( mindistsq2 > dsq2 )
883 {
884 result.clear();
885 result.push_back(to);
886 }
887 else if ( are_near(mindistsq2, dsq2) )
888 {
889 result.push_back(to);
890 }
892 }
893 else
894 {
895 if ( result.empty() )
896 {
897 if ( are_near(dsq1, dsq2) )
898 {
899 result.push_back(from);
900 result.push_back(to);
901 }
902 else if ( dsq2 > dsq1 )
903 {
904 result.push_back(from);
905 }
906 else
907 {
908 result.push_back(to);
909 }
910 }
911 }
913 return result;
914 }
917 } // end namespace Geom
920 /*
921 Local Variables:
922 mode:c++
923 c-file-style:"stroustrup"
924 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
925 indent-tabs-mode:nil
926 fill-column:99
927 End:
928 */
929 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :