src/2geom/elliptical-arc.cpp

   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  */
  29
  30
  31 #include <2geom/elliptical-arc.h>
  32 #include <2geom/bezier-curve.h>
  33 #include <2geom/poly.h>
  34 #include <2geom/sbasis-math.h>
  35
  36 #include <cfloat>
  37 #include <limits>
  38
  39
  40
  41
  42 namespace Geom
  43 {
  44
  45
  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;
  57
  58
  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]);
  68
  69
  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         }
 116
 117         return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
 118                              Point(arc_extremes[0], arc_extremes[2]) );
 119
 120 }
 121
 122
 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         }
 130
 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         }
 138
 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         };
 156
 157         for ( unsigned int dim = 0; dim < 2; ++dim )
 158         {
 159                 if ( are_near(ray(dim), 0) )
 160                 {
 161
 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                         }
 198
 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                         }
 227
 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         }
 236
 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;
 251
 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;
 258
 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         }
 290
 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;
 301
 302
 303 //      return SBasisCurve(toSBasis()).roots(v, d);
 304 }
 305
 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;
 324
 325 }
 326
 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 }
 359
 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 }
 377
 378
 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 }
 392
 393
 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 }
 412
 413
 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 }
 448
 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;
 456
 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                 }
 521
 522         }
 523
 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                 }
 570
 571         }
 572
 573     double sin_rot_angle = std::sin(rotation_angle());
 574     double cos_rot_angle = std::cos(rotation_angle());
 575
 576     Point sp = sweep_flag() ? initialPoint() : finalPoint();
 577     Point ep = sweep_flag() ? finalPoint() : initialPoint();
 578
 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         }
 606
 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     }
 629
 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;
 637
 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 }
 653
 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 }
 670
 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 }
 676
 677
 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         }
 751
 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          */
 766
 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];
 779
 780 //      for ( unsigned int i = 0; i < 5; ++i )
 781 //              std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
 782
 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 //      }
 816
 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         }
 828
 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         }
 849
 850         double t = map_to_01( real_sol[mi1] );
 851         if ( !(t < from || t > to) )
 852         {
 853                 result.push_back(t);
 854         }
 855
 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         }
 866
 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                 }
 891
 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         }
 912
 913         return result;
 914 }
 915
 916
 917 } // end namespace Geom
 918
 919
 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 :
 930
 931