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