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 "elliptical-arc.h"
  32 #include "bezier-curve.h"
  33 #include "poly.h"
  34
  35 #include <cfloat>
  36 #include <limits>
  37
  38
  39
  40
  41 namespace Geom
  42 {
  43
  44
  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;
  56
  57
  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]);
  67
  68
  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         }
 115
 116         return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
 117                              Point(arc_extremes[0], arc_extremes[2]) );
 118
 119 }
 120
 121
 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         }
 129
 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         }
 137
 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         };
 155
 156         for ( unsigned int dim = 0; dim < 2; ++dim )
 157         {
 158                 if ( are_near(ray(dim), 0) )
 159                 {
 160
 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                         }
 197
 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                         }
 226
 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         }
 235
 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;
 250
 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;
 257
 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         }
 289
 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;
 300
 301
 302 //      return SBasisCurve(toSBasis()).roots(v, d);
 303 }
 304
 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;
 323
 324 }
 325
 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 }
 358
 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 }
 376
 377
 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 }
 391
 392
 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 }
 411
 412
 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 }
 447
 448 // NOTE: doesn't work with 360 deg arcs
 449 void EllipticalArc::calculate_center_and_extreme_angles()
 450 {
 451     if ( are_near(initialPoint(), finalPoint()) )
 452     {
 453         if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
 454         {
 455                 m_start_angle = m_end_angle = 0;
 456                 m_center = initialPoint();
 457                 return;
 458         }
 459         else
 460         {
 461                 THROW_RANGEERROR("initial and final point are the same");
 462         }
 463     }
 464         if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
 465         { // but initialPoint != finalPoint
 466                 THROW_RANGEERROR(
 467                         "there is no ellipse that satisfies the given constraints: "
 468                         "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"
 469                 );
 470         }
 471         if ( are_near(ray(Y), 0) )
 472         {
 473                 Point v = initialPoint() - finalPoint();
 474                 if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )
 475                 {
 476                         double angle = std::atan2(v[Y], v[X]);
 477                         if (angle < 0) angle += 2*M_PI;
 478                         if ( are_near( angle, rotation_angle() ) )
 479                         {
 480                                 m_start_angle = 0;
 481                                 m_end_angle = M_PI;
 482                                 m_center = v/2 + finalPoint();
 483                                 return;
 484                         }
 485                         angle -= M_PI;
 486                         if ( angle < 0 ) angle += 2*M_PI;
 487                         if ( are_near( angle, rotation_angle() ) )
 488                         {
 489                                 m_start_angle = M_PI;
 490                                 m_end_angle = 0;
 491                                 m_center = v/2 + finalPoint();
 492                                 return;
 493                         }
 494                         THROW_RANGEERROR(
 495                                 "there is no ellipse that satisfies the given constraints: "
 496                                 "ray(Y) == 0 "
 497                                 "and slope(initialPoint - finalPoint) != rotation_angle "
 498                                 "and != rotation_angle + PI"
 499                         );
 500                 }
 501                 if ( L2sq(v) > 4*ray(X)*ray(X) )
 502                 {
 503                         THROW_RANGEERROR(
 504                                 "there is no ellipse that satisfies the given constraints: "
 505                                 "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"
 506                         );
 507                 }
 508                 else
 509                 {
 510                         THROW_RANGEERROR(
 511                                 "there is infinite ellipses that satisfy the given constraints: "
 512                                 "ray(Y) == 0  and distance(initialPoint, finalPoint) < 2*ray(X)"
 513                         );
 514                 }
 515
 516         }
 517
 518         if ( are_near(ray(X), 0) )
 519         {
 520                 Point v = initialPoint() - finalPoint();
 521                 if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )
 522                 {
 523                         double angle = std::atan2(v[Y], v[X]);
 524                         if (angle < 0) angle += 2*M_PI;
 525                         double rot_angle = rotation_angle() + M_PI/2;
 526                         if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;
 527                         if ( are_near( angle, rot_angle ) )
 528                         {
 529                                 m_start_angle = M_PI/2;
 530                                 m_end_angle = 3*M_PI/2;
 531                                 m_center = v/2 + finalPoint();
 532                                 return;
 533                         }
 534                         angle -= M_PI;
 535                         if ( angle < 0 ) angle += 2*M_PI;
 536                         if ( are_near( angle, rot_angle ) )
 537                         {
 538                                 m_start_angle = 3*M_PI/2;
 539                                 m_end_angle = M_PI/2;
 540                                 m_center = v/2 + finalPoint();
 541                                 return;
 542                         }
 543                         THROW_RANGEERROR(
 544                                 "there is no ellipse that satisfies the given constraints: "
 545                                 "ray(X) == 0 "
 546                                 "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "
 547                                 "and != rotation_angle + (3/2)*PI"
 548                         );
 549                 }
 550                 if ( L2sq(v) > 4*ray(Y)*ray(Y) )
 551                 {
 552                         THROW_RANGEERROR(
 553                                 "there is no ellipse that satisfies the given constraints: "
 554                                 "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"
 555                         );
 556                 }
 557                 else
 558                 {
 559                         THROW_RANGEERROR(
 560                                 "there is infinite ellipses that satisfy the given constraints: "
 561                                 "ray(X) == 0  and distance(initialPoint, finalPoint) < 2*ray(Y)"
 562                         );
 563                 }
 564
 565         }
 566
 567     double sin_rot_angle = std::sin(rotation_angle());
 568     double cos_rot_angle = std::cos(rotation_angle());
 569
 570     Point sp = sweep_flag() ? initialPoint() : finalPoint();
 571     Point ep = sweep_flag() ? finalPoint() : initialPoint();
 572
 573     Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,
 574              -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,
 575               0,                      0 );
 576     Matrix im = m.inverse();
 577     Point sol = (ep - sp) * im;
 578     double half_sum_angle = std::atan2(-sol[X], sol[Y]);
 579     double half_diff_angle;
 580     if ( are_near(std::fabs(half_sum_angle), M_PI/2) )
 581     {
 582         double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1;
 583         double arg = anti_sgn_hsa * sol[X] / 2;
 584         // if |arg| is a little bit > 1 acos returns nan
 585         if ( are_near(arg, 1) )
 586             half_diff_angle = 0;
 587         else if ( are_near(arg, -1) )
 588             half_diff_angle = M_PI;
 589         else
 590         {
 591                 if ( !(-1 < arg && arg < 1) )
 592                         THROW_RANGEERROR(
 593                                 "there is no ellipse that satisfies the given constraints"
 594                         );
 595             // assert( -1 < arg && arg < 1 );
 596             // if it fails
 597             // => there is no ellipse that satisfies the given constraints
 598             half_diff_angle = std::acos( arg );
 599         }
 600
 601         half_diff_angle = M_PI/2 - half_diff_angle;
 602     }
 603     else
 604     {
 605         double  arg = sol[Y] / ( 2 * std::cos(half_sum_angle) );
 606         // if |arg| is a little bit > 1 asin returns nan
 607         if ( are_near(arg, 1) )
 608             half_diff_angle = M_PI/2;
 609         else if ( are_near(arg, -1) )
 610             half_diff_angle = -M_PI/2;
 611         else
 612         {
 613                 if ( !(-1 < arg && arg < 1) )
 614                         THROW_RANGEERROR(
 615                                 "there is no ellipse that satisfies the given constraints"
 616                         );
 617             // assert( -1 < arg && arg < 1 );
 618             // if it fails
 619             // => there is no ellipse that satisfies the given constraints
 620             half_diff_angle = std::asin( arg );
 621         }
 622     }
 623
 624     if (   ( m_large_arc && half_diff_angle > 0 )
 625         || (!m_large_arc && half_diff_angle < 0 ) )
 626     {
 627         half_diff_angle = -half_diff_angle;
 628     }
 629     if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI;
 630     if ( half_diff_angle < 0 ) half_diff_angle += M_PI;
 631
 632     m_start_angle = half_sum_angle - half_diff_angle;
 633     m_end_angle =  half_sum_angle + half_diff_angle;
 634     // 0 <= m_start_angle, m_end_angle < 2PI
 635     if ( m_start_angle < 0 ) m_start_angle += 2*M_PI;
 636     if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI;
 637     sol[0] = std::cos(m_start_angle);
 638     sol[1] = std::sin(m_start_angle);
 639     m_center = sp - sol * m;
 640     if ( !sweep_flag() )
 641     {
 642         double angle = m_start_angle;
 643         m_start_angle = m_end_angle;
 644         m_end_angle = angle;
 645     }
 646 }
 647
 648 Coord EllipticalArc::map_to_02PI(Coord t) const
 649 {
 650     if ( sweep_flag() )
 651     {
 652         Coord angle = start_angle() + sweep_angle() * t;
 653         if ( !(angle < 2*M_PI) )
 654             angle -= 2*M_PI;
 655         return angle;
 656     }
 657     else
 658     {
 659         Coord angle = start_angle() - sweep_angle() * t;
 660         if ( angle < 0 ) angle += 2*M_PI;
 661         return angle;
 662     }
 663 }
 664
 665 Coord EllipticalArc::map_to_01(Coord angle) const
 666 {
 667         return map_circular_arc_on_unit_interval(angle, start_angle(),
 668                                                          end_angle(), sweep_flag());
 669 }
 670
 671
 672 std::vector<double> EllipticalArc::
 673 allNearestPoints( Point const& p, double from, double to ) const
 674 {
 675         if ( from > to ) std::swap(from, to);
 676         if ( from < 0 || to > 1 )
 677         {
 678                 THROW_RANGEERROR("[from,to] interval out of range");
 679         }
 680         std::vector<double> result;
 681         if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) )  || are_near(from, to) )
 682         {
 683                 result.push_back(from);
 684                 return result;
 685         }
 686         else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
 687         {
 688                 LineSegment seg(pointAt(from), pointAt(to));
 689                 Point np = seg.pointAt( seg.nearestPoint(p) );
 690                 if ( are_near(ray(Y), 0) )
 691                 {
 692                         if ( are_near(rotation_angle(), M_PI/2)
 693                                  || are_near(rotation_angle(), 3*M_PI/2) )
 694                         {
 695                                 result = roots(np[Y], Y);
 696                         }
 697                         else
 698                         {
 699                                 result = roots(np[X], X);
 700                         }
 701                 }
 702                 else
 703                 {
 704                         if ( are_near(rotation_angle(), M_PI/2)
 705                                  || are_near(rotation_angle(), 3*M_PI/2) )
 706                         {
 707                                 result = roots(np[X], X);
 708                         }
 709                         else
 710                         {
 711                                 result = roots(np[Y], Y);
 712                         }
 713                 }
 714                 return result;
 715         }
 716         else if ( are_near(ray(X), ray(Y)) )
 717         {
 718                 Point r = p - center();
 719                 if ( are_near(r, Point(0,0)) )
 720                 {
 721                         THROW_INFINITESOLUTIONS(0);
 722                 }
 723                 // TODO: implement case r != 0
 724 //              Point np = ray(X) * unit_vector(r);
 725 //              std::vector<double> solX = roots(np[X],X);
 726 //              std::vector<double> solY = roots(np[Y],Y);
 727 //              double t;
 728 //              if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
 729 //              {
 730 //                      t = solX[0];
 731 //              }
 732 //              else
 733 //              {
 734 //                      t = solX[1];
 735 //              }
 736 //              if ( !(t < from || t > to) )
 737 //              {
 738 //                      result.push_back(t);
 739 //              }
 740 //              else
 741 //              {
 742 //
 743 //              }
 744         }
 745
 746         // solve the equation <D(E(t),t)|E(t)-p> == 0
 747         // that provides min and max distance points
 748         // on the ellipse E wrt the point p
 749         // after the substitutions:
 750         // cos(t) = (1 - s^2) / (1 + s^2)
 751         // sin(t) = 2t / (1 + s^2)
 752         // where s = tan(t/2)
 753         // we get a 4th degree equation in s
 754         /*
 755          *      ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
 756          *      ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
 757          *      2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
 758          *      2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
 759          */
 760
 761         Point p_c = p - center();
 762         double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
 763         double cosrot = std::cos( rotation_angle() );
 764         double sinrot = std::sin( rotation_angle() );
 765         double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
 766         Poly coeff;
 767         coeff.resize(5);
 768         coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
 769         coeff[3] = 2 * ( rx2_ry2 + expr1 );
 770         coeff[2] = 0;
 771         coeff[1] = 2 * ( -rx2_ry2 + expr1 );
 772         coeff[0] = -coeff[4];
 773
 774 //      for ( unsigned int i = 0; i < 5; ++i )
 775 //              std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
 776
 777         std::vector<double> real_sol;
 778         // gsl_poly_complex_solve raises an error
 779         // if the leading coefficient is zero
 780         if ( are_near(coeff[4], 0) )
 781         {
 782                 real_sol.push_back(0);
 783                 if ( !are_near(coeff[3], 0) )
 784                 {
 785                         double sq = -coeff[1] / coeff[3];
 786                         if ( sq > 0 )
 787                         {
 788                                 double s = std::sqrt(sq);
 789                                 real_sol.push_back(s);
 790                                 real_sol.push_back(-s);
 791                         }
 792                 }
 793         }
 794         else
 795         {
 796                 real_sol = solve_reals(coeff);
 797         }
 798 //      else
 799 //      {
 800 //              double sol[8];
 801 //              gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5);
 802 //              gsl_poly_complex_solve(coeff, 5, w, sol );
 803 //              gsl_poly_complex_workspace_free(w);
 804 //
 805 //              for ( unsigned int i = 0; i < 4; ++i )
 806 //              {
 807 //                      if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]);
 808 //              }
 809 //      }
 810
 811         for ( unsigned int i = 0; i < real_sol.size(); ++i )
 812         {
 813                 real_sol[i] = 2 * std::atan(real_sol[i]);
 814                 if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;
 815         }
 816         // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
 817         // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
 818         if ( (real_sol.size() % 2) != 0 )
 819         {
 820                 real_sol.push_back(M_PI);
 821         }
 822
 823         double mindistsq1 = std::numeric_limits<double>::max();
 824         double mindistsq2 = std::numeric_limits<double>::max();
 825         double dsq;
 826         unsigned int mi1, mi2;
 827         for ( unsigned int i = 0; i < real_sol.size(); ++i )
 828         {
 829                 dsq = distanceSq(p, pointAtAngle(real_sol[i]));
 830                 if ( mindistsq1 > dsq )
 831                 {
 832                         mindistsq2 = mindistsq1;
 833                         mi2 = mi1;
 834                         mindistsq1 = dsq;
 835                         mi1 = i;
 836                 }
 837                 else if ( mindistsq2 > dsq )
 838                 {
 839                         mindistsq2 = dsq;
 840                         mi2 = i;
 841                 }
 842         }
 843
 844         double t = map_to_01( real_sol[mi1] );
 845         if ( !(t < from || t > to) )
 846         {
 847                 result.push_back(t);
 848         }
 849
 850         bool second_sol = false;
 851         t = map_to_01( real_sol[mi2] );
 852         if ( real_sol.size() == 4 && !(t < from || t > to) )
 853         {
 854         if ( result.empty() || are_near(mindistsq1, mindistsq2) )
 855         {
 856                 result.push_back(t);
 857                 second_sol = true;
 858         }
 859         }
 860
 861         // we need to test extreme points too
 862         double dsq1 = distanceSq(p, pointAt(from));
 863         double dsq2 = distanceSq(p, pointAt(to));
 864         if ( second_sol )
 865         {
 866                 if ( mindistsq2 > dsq1 )
 867                 {
 868                         result.clear();
 869                         result.push_back(from);
 870                         mindistsq2 = dsq1;
 871                 }
 872                 else if ( are_near(mindistsq2, dsq) )
 873                 {
 874                         result.push_back(from);
 875                 }
 876                 if ( mindistsq2 > dsq2 )
 877                 {
 878                         result.clear();
 879                         result.push_back(to);
 880                 }
 881                 else if ( are_near(mindistsq2, dsq2) )
 882                 {
 883                         result.push_back(to);
 884                 }
 885
 886         }
 887         else
 888         {
 889                 if ( result.empty() )
 890                 {
 891                         if ( are_near(dsq1, dsq2) )
 892                         {
 893                                 result.push_back(from);
 894                                 result.push_back(to);
 895                         }
 896                         else if ( dsq2 > dsq1 )
 897                         {
 898                                 result.push_back(from);
 899                         }
 900                         else
 901                         {
 902                                 result.push_back(to);
 903                         }
 904                 }
 905         }
 906
 907         return result;
 908 }
 909
 910
 911 } // end namespace Geom
 912
 913
 914 /*
 915   Local Variables:
 916   mode:c++
 917   c-file-style:"stroustrup"
 918   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 919   indent-tabs-mode:nil
 920   fill-column:99
 921   End:
 922 */
 923 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
 924
 925