src/2geom/svg-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/svg-elliptical-arc.h>
  32 #include <2geom/ellipse.h>
  33 #include <2geom/sbasis-geometric.h>
  34 #include <2geom/bezier-curve.h>
  35 #include <2geom/poly.h>
  36
  37 #include <cfloat>
  38 #include <limits>
  39
  40 #include <2geom/numeric/vector.h>
  41 #include <2geom/numeric/fitting-tool.h>
  42 #include <2geom/numeric/fitting-model.h>
  43
  44
  45
  46 namespace Geom
  47 {
  48
  49
  50 OptRect SVGEllipticalArc::boundsExact() const
  51 {
  52     if (isDegenerate() && is_svg_compliant())
  53         return chord().boundsExact();
  54
  55     std::vector<double> extremes(4);
  56     double cosrot = std::cos(rotation_angle());
  57     double sinrot = std::sin(rotation_angle());
  58     extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
  59     extremes[1] = extremes[0] + M_PI;
  60     if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;
  61     extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
  62     extremes[3] = extremes[2] + M_PI;
  63     if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;
  64
  65
  66     std::vector<double>arc_extremes(4);
  67     arc_extremes[0] = initialPoint()[X];
  68     arc_extremes[1] = finalPoint()[X];
  69     if ( arc_extremes[0] < arc_extremes[1] )
  70         std::swap(arc_extremes[0], arc_extremes[1]);
  71     arc_extremes[2] = initialPoint()[Y];
  72     arc_extremes[3] = finalPoint()[Y];
  73     if ( arc_extremes[2] < arc_extremes[3] )
  74         std::swap(arc_extremes[2], arc_extremes[3]);
  75
  76
  77     if ( start_angle() < end_angle() )
  78     {
  79         if ( sweep_flag() )
  80         {
  81             for ( unsigned int i = 0; i < extremes.size(); ++i )
  82             {
  83                 if ( start_angle() < extremes[i] && extremes[i] < end_angle() )
  84                 {
  85                     arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
  86                 }
  87             }
  88         }
  89         else
  90         {
  91             for ( unsigned int i = 0; i < extremes.size(); ++i )
  92             {
  93                 if ( start_angle() > extremes[i] || extremes[i] > end_angle() )
  94                 {
  95                     arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
  96                 }
  97             }
  98         }
  99     }
 100     else
 101     {
 102         if ( sweep_flag() )
 103         {
 104             for ( unsigned int i = 0; i < extremes.size(); ++i )
 105             {
 106                 if ( start_angle() < extremes[i] || extremes[i] < end_angle() )
 107                 {
 108                     arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
 109                 }
 110             }
 111         }
 112         else
 113         {
 114             for ( unsigned int i = 0; i < extremes.size(); ++i )
 115             {
 116                 if ( start_angle() > extremes[i] && extremes[i] > end_angle() )
 117                 {
 118                     arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
 119                 }
 120             }
 121         }
 122     }
 123
 124     return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
 125                  Point(arc_extremes[0], arc_extremes[2]) );
 126 }
 127
 128
 129 double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const
 130 {
 131     double sin_rot_angle = std::sin(rotation_angle());
 132     double cos_rot_angle = std::cos(rotation_angle());
 133     if ( d == X )
 134     {
 135         return    ray(X) * cos_rot_angle * std::cos(t)
 136                 - ray(Y) * sin_rot_angle * std::sin(t)
 137                 + center(X);
 138     }
 139     else if ( d == Y )
 140     {
 141         return    ray(X) * sin_rot_angle * std::cos(t)
 142                 + ray(Y) * cos_rot_angle * std::sin(t)
 143                 + center(Y);
 144     }
 145     THROW_RANGEERROR("dimension parameter out of range");
 146 }
 147
 148
 149 std::vector<double>
 150 SVGEllipticalArc::roots(double v, Dim2 d) const
 151 {
 152     if ( d > Y )
 153     {
 154         THROW_RANGEERROR("dimention out of range");
 155     }
 156
 157     std::vector<double> sol;
 158
 159     if (isDegenerate() && is_svg_compliant())
 160     {
 161         return chord().roots(v, d);
 162     }
 163     else
 164     {
 165         if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
 166         {
 167             if ( center(d) == v )
 168                 sol.push_back(0);
 169             return sol;
 170         }
 171
 172         const char* msg[2][2] =
 173         {
 174             { "d == X; ray(X) == 0; "
 175               "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "
 176               "s should be contained in [-1,1]",
 177               "d == X; ray(Y) == 0; "
 178               "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "
 179               "s should be contained in [-1,1]"
 180             },
 181             { "d == Y; ray(X) == 0; "
 182               "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "
 183               "s should be contained in [-1,1]",
 184               "d == Y; ray(Y) == 0; "
 185               "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "
 186               "s should be contained in [-1,1]"
 187             },
 188         };
 189
 190         for ( unsigned int dim = 0; dim < 2; ++dim )
 191         {
 192             if ( are_near(ray(dim), 0) )
 193             {
 194                 if ( initialPoint()[d] == v && finalPoint()[d] == v )
 195                 {
 196                     THROW_INFINITESOLUTIONS(0);
 197                 }
 198                 if ( (initialPoint()[d] < finalPoint()[d])
 199                      && (initialPoint()[d] > v || finalPoint()[d] < v) )
 200                 {
 201                     return sol;
 202                 }
 203                 if ( (initialPoint()[d] > finalPoint()[d])
 204                      && (finalPoint()[d] > v || initialPoint()[d] < v) )
 205                 {
 206                     return sol;
 207                 }
 208                 double ray_prj;
 209                 switch(d)
 210                 {
 211                     case X:
 212                         switch(dim)
 213                         {
 214                             case X: ray_prj = -ray(Y) * std::sin(rotation_angle());
 215                                     break;
 216                             case Y: ray_prj = ray(X) * std::cos(rotation_angle());
 217                                     break;
 218                         }
 219                         break;
 220                     case Y:
 221                         switch(dim)
 222                         {
 223                             case X: ray_prj = ray(Y) * std::cos(rotation_angle());
 224                                     break;
 225                             case Y: ray_prj = ray(X) * std::sin(rotation_angle());
 226                                     break;
 227                         }
 228                         break;
 229                 }
 230
 231                 double s = (v - center(d)) / ray_prj;
 232                 if ( s < -1 || s > 1 )
 233                 {
 234                     THROW_LOGICALERROR(msg[d][dim]);
 235                 }
 236                 switch(dim)
 237                 {
 238                     case X:
 239                         s = std::asin(s); // return a value in [-PI/2,PI/2]
 240                         if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) )  )
 241                         {
 242                             if ( s < 0 ) s += 2*M_PI;
 243                         }
 244                         else
 245                         {
 246                             s = M_PI - s;
 247                             if (!(s < 2*M_PI) ) s -= 2*M_PI;
 248                         }
 249                         break;
 250                     case Y:
 251                         s = std::acos(s); // return a value in [0,PI]
 252                         if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )
 253                         {
 254                             s = 2*M_PI - s;
 255                             if ( !(s < 2*M_PI) ) s -= 2*M_PI;
 256                         }
 257                         break;
 258                 }
 259
 260                 //std::cerr << "s = " << rad_to_deg(s);
 261                 s = map_to_01(s);
 262                 //std::cerr << " -> t: " << s << std::endl;
 263                 if ( !(s < 0 || s > 1) )
 264                     sol.push_back(s);
 265                 return sol;
 266             }
 267         }
 268
 269     }
 270
 271     double rotx, roty;
 272     switch(d)
 273     {
 274         case X:
 275             rotx = std::cos(rotation_angle());
 276             roty = -std::sin(rotation_angle());
 277             break;
 278         case Y:
 279             rotx = std::sin(rotation_angle());
 280             roty = std::cos(rotation_angle());
 281             break;
 282     }
 283     double rxrotx = ray(X) * rotx;
 284     double c_v = center(d) - v;
 285
 286     double a = -rxrotx + c_v;
 287     double b = ray(Y) * roty;
 288     double c = rxrotx + c_v;
 289     //std::cerr << "a = " << a << std::endl;
 290     //std::cerr << "b = " << b << std::endl;
 291     //std::cerr << "c = " << c << std::endl;
 292
 293     if ( are_near(a,0) )
 294     {
 295         sol.push_back(M_PI);
 296         if ( !are_near(b,0) )
 297         {
 298             double s = 2 * std::atan(-c/(2*b));
 299             if ( s < 0 ) s += 2*M_PI;
 300             sol.push_back(s);
 301         }
 302     }
 303     else
 304     {
 305         double delta = b * b - a * c;
 306         //std::cerr << "delta = " << delta << std::endl;
 307         if ( are_near(delta, 0) )
 308         {
 309             double s = 2 * std::atan(-b/a);
 310             if ( s < 0 ) s += 2*M_PI;
 311             sol.push_back(s);
 312         }
 313         else if ( delta > 0 )
 314         {
 315             double sq = std::sqrt(delta);
 316             double s = 2 * std::atan( (-b - sq) / a );
 317             if ( s < 0 ) s += 2*M_PI;
 318             sol.push_back(s);
 319             s = 2 * std::atan( (-b + sq) / a );
 320             if ( s < 0 ) s += 2*M_PI;
 321             sol.push_back(s);
 322         }
 323     }
 324
 325     std::vector<double> arc_sol;
 326     for (unsigned int i = 0; i < sol.size(); ++i )
 327     {
 328         //std::cerr << "s = " << rad_to_deg(sol[i]);
 329         sol[i] = map_to_01(sol[i]);
 330         //std::cerr << " -> t: " << sol[i] << std::endl;
 331         if ( !(sol[i] < 0 || sol[i] > 1) )
 332             arc_sol.push_back(sol[i]);
 333     }
 334     return arc_sol;
 335 }
 336
 337
 338 // D(E(t,C),t) = E(t+PI/2,O), where C is the ellipse center
 339 // the derivative doesn't rotate the ellipse but there is a translation
 340 // of the parameter t by an angle of PI/2 so the ellipse points are shifted
 341 // of such an angle in the cw direction
 342 Curve* SVGEllipticalArc::derivative() const
 343 {
 344     if (isDegenerate() && is_svg_compliant())
 345             return chord().derivative();
 346
 347     SVGEllipticalArc* result = new SVGEllipticalArc(*this);
 348     result->m_center[X] = result->m_center[Y] = 0;
 349     result->m_start_angle += M_PI/2;
 350     if( !( result->m_start_angle < 2*M_PI ) )
 351     {
 352         result->m_start_angle -= 2*M_PI;
 353     }
 354     result->m_end_angle += M_PI/2;
 355     if( !( result->m_end_angle < 2*M_PI ) )
 356     {
 357         result->m_end_angle -= 2*M_PI;
 358     }
 359     result->m_initial_point = result->pointAtAngle( result->start_angle() );
 360     result->m_final_point = result->pointAtAngle( result->end_angle() );
 361     return result;
 362 }
 363
 364
 365 std::vector<Point>
 366 SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
 367 {
 368     if (isDegenerate() && is_svg_compliant())
 369             return chord().pointAndDerivatives(t, n);
 370
 371     unsigned int nn = n+1; // nn represents the size of the result vector.
 372     std::vector<Point> result;
 373     result.reserve(nn);
 374     double angle = map_unit_interval_on_circular_arc(t, start_angle(),
 375                                                      end_angle(), sweep_flag());
 376     SVGEllipticalArc ea(*this);
 377     ea.m_center = Point(0,0);
 378     unsigned int m = std::min(nn, 4u);
 379     for ( unsigned int i = 0; i < m; ++i )
 380     {
 381         result.push_back( ea.pointAtAngle(angle) );
 382         angle += M_PI/2;
 383         if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
 384     }
 385     m = nn / 4;
 386     for ( unsigned int i = 1; i < m; ++i )
 387     {
 388         for ( unsigned int j = 0; j < 4; ++j )
 389             result.push_back( result[j] );
 390     }
 391     m = nn - 4 * m;
 392     for ( unsigned int i = 0; i < m; ++i )
 393     {
 394         result.push_back( result[i] );
 395     }
 396     if ( !result.empty() ) // nn != 0
 397         result[0] = pointAtAngle(angle);
 398     return result;
 399 }
 400
 401 bool SVGEllipticalArc::containsAngle(Coord angle) const
 402 {
 403     if ( sweep_flag() )
 404         if ( start_angle() < end_angle() )
 405             return ( !( angle < start_angle() || angle > end_angle() ) );
 406         else
 407             return ( !( angle < start_angle() && angle > end_angle() ) );
 408     else
 409         if ( start_angle() > end_angle() )
 410             return ( !( angle > start_angle() || angle < end_angle() ) );
 411         else
 412             return ( !( angle > start_angle() && angle < end_angle() ) );
 413 }
 414
 415 Curve* SVGEllipticalArc::portion(double f, double t) const
 416 {
 417     // fix input arguments
 418     if (f < 0) f = 0;
 419     if (f > 1) f = 1;
 420     if (t < 0) t = 0;
 421     if (t > 1) t = 1;
 422
 423     if ( are_near(f, t) )
 424     {
 425         SVGEllipticalArc* arc = new SVGEllipticalArc();
 426         arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);
 427         arc->m_start_angle = arc->m_end_angle = m_start_angle;
 428         arc->m_rot_angle = m_rot_angle;
 429         arc->m_sweep = m_sweep;
 430         arc->m_large_arc = m_large_arc;
 431     }
 432
 433     SVGEllipticalArc* arc = new SVGEllipticalArc( *this );
 434     arc->m_initial_point = pointAt(f);
 435     arc->m_final_point = pointAt(t);
 436     double sa = sweep_flag() ? sweep_angle() : -sweep_angle();
 437     arc->m_start_angle = m_start_angle + sa * f;
 438     if ( !(arc->m_start_angle < 2*M_PI) )
 439         arc->m_start_angle -= 2*M_PI;
 440     if ( arc->m_start_angle < 0 )
 441         arc->m_start_angle += 2*M_PI;
 442     arc->m_end_angle = m_start_angle + sa * t;
 443     if ( !(arc->m_end_angle < 2*M_PI) )
 444         arc->m_end_angle -= 2*M_PI;
 445     if ( arc->m_end_angle < 0 )
 446         arc->m_end_angle += 2*M_PI;
 447     if ( f > t ) arc->m_sweep = !sweep_flag();
 448     if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )
 449         arc->m_large_arc = false;
 450     return arc;
 451 }
 452
 453
 454 std::vector<double> SVGEllipticalArc::
 455 allNearestPoints( Point const& p, double from, double to ) const
 456 {
 457     std::vector<double> result;
 458     if (isDegenerate() && is_svg_compliant())
 459     {
 460         result.push_back( chord().nearestPoint(p, from, to) );
 461         return result;
 462     }
 463
 464     if ( from > to ) std::swap(from, to);
 465     if ( from < 0 || to > 1 )
 466     {
 467         THROW_RANGEERROR("[from,to] interval out of range");
 468     }
 469
 470     if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) )  || are_near(from, to) )
 471     {
 472         result.push_back(from);
 473         return result;
 474     }
 475     else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
 476     {
 477         LineSegment seg(pointAt(from), pointAt(to));
 478         Point np = seg.pointAt( seg.nearestPoint(p) );
 479         if ( are_near(ray(Y), 0) )
 480         {
 481             if ( are_near(rotation_angle(), M_PI/2)
 482                  || are_near(rotation_angle(), 3*M_PI/2) )
 483             {
 484                 result = roots(np[Y], Y);
 485             }
 486             else
 487             {
 488                 result = roots(np[X], X);
 489             }
 490         }
 491         else
 492         {
 493             if ( are_near(rotation_angle(), M_PI/2)
 494                  || are_near(rotation_angle(), 3*M_PI/2) )
 495             {
 496                 result = roots(np[X], X);
 497             }
 498             else
 499             {
 500                 result = roots(np[Y], Y);
 501             }
 502         }
 503         return result;
 504     }
 505     else if ( are_near(ray(X), ray(Y)) )
 506     {
 507         Point r = p - center();
 508         if ( are_near(r, Point(0,0)) )
 509         {
 510             THROW_INFINITESOLUTIONS(0);
 511         }
 512         // TODO: implement case r != 0
 513 //      Point np = ray(X) * unit_vector(r);
 514 //      std::vector<double> solX = roots(np[X],X);
 515 //      std::vector<double> solY = roots(np[Y],Y);
 516 //      double t;
 517 //      if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
 518 //      {
 519 //          t = solX[0];
 520 //      }
 521 //      else
 522 //      {
 523 //          t = solX[1];
 524 //      }
 525 //      if ( !(t < from || t > to) )
 526 //      {
 527 //          result.push_back(t);
 528 //      }
 529 //      else
 530 //      {
 531 //
 532 //      }
 533     }
 534
 535     // solve the equation <D(E(t),t)|E(t)-p> == 0
 536     // that provides min and max distance points
 537     // on the ellipse E wrt the point p
 538     // after the substitutions:
 539     // cos(t) = (1 - s^2) / (1 + s^2)
 540     // sin(t) = 2t / (1 + s^2)
 541     // where s = tan(t/2)
 542     // we get a 4th degree equation in s
 543     /*
 544      *  ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
 545      *  ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
 546      *  2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
 547      *  2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
 548      */
 549
 550     Point p_c = p - center();
 551     double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
 552     double cosrot = std::cos( rotation_angle() );
 553     double sinrot = std::sin( rotation_angle() );
 554     double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
 555     Poly coeff;
 556     coeff.resize(5);
 557     coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
 558     coeff[3] = 2 * ( rx2_ry2 + expr1 );
 559     coeff[2] = 0;
 560     coeff[1] = 2 * ( -rx2_ry2 + expr1 );
 561     coeff[0] = -coeff[4];
 562
 563 //  for ( unsigned int i = 0; i < 5; ++i )
 564 //      std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
 565
 566     std::vector<double> real_sol;
 567     // gsl_poly_complex_solve raises an error
 568     // if the leading coefficient is zero
 569     if ( are_near(coeff[4], 0) )
 570     {
 571         real_sol.push_back(0);
 572         if ( !are_near(coeff[3], 0) )
 573         {
 574             double sq = -coeff[1] / coeff[3];
 575             if ( sq > 0 )
 576             {
 577                 double s = std::sqrt(sq);
 578                 real_sol.push_back(s);
 579                 real_sol.push_back(-s);
 580             }
 581         }
 582     }
 583     else
 584     {
 585         real_sol = solve_reals(coeff);
 586     }
 587
 588     for ( unsigned int i = 0; i < real_sol.size(); ++i )
 589     {
 590         real_sol[i] = 2 * std::atan(real_sol[i]);
 591         if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;
 592     }
 593     // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
 594     // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
 595     if ( (real_sol.size() % 2) != 0 )
 596     {
 597         real_sol.push_back(M_PI);
 598     }
 599
 600     double mindistsq1 = std::numeric_limits<double>::max();
 601     double mindistsq2 = std::numeric_limits<double>::max();
 602     double dsq;
 603     unsigned int mi1, mi2;
 604     for ( unsigned int i = 0; i < real_sol.size(); ++i )
 605     {
 606         dsq = distanceSq(p, pointAtAngle(real_sol[i]));
 607         if ( mindistsq1 > dsq )
 608         {
 609             mindistsq2 = mindistsq1;
 610             mi2 = mi1;
 611             mindistsq1 = dsq;
 612             mi1 = i;
 613         }
 614         else if ( mindistsq2 > dsq )
 615         {
 616             mindistsq2 = dsq;
 617             mi2 = i;
 618         }
 619     }
 620
 621     double t = map_to_01( real_sol[mi1] );
 622     if ( !(t < from || t > to) )
 623     {
 624         result.push_back(t);
 625     }
 626
 627     bool second_sol = false;
 628     t = map_to_01( real_sol[mi2] );
 629     if ( real_sol.size() == 4 && !(t < from || t > to) )
 630     {
 631         if ( result.empty() || are_near(mindistsq1, mindistsq2) )
 632         {
 633             result.push_back(t);
 634             second_sol = true;
 635         }
 636     }
 637
 638     // we need to test extreme points too
 639     double dsq1 = distanceSq(p, pointAt(from));
 640     double dsq2 = distanceSq(p, pointAt(to));
 641     if ( second_sol )
 642     {
 643         if ( mindistsq2 > dsq1 )
 644         {
 645             result.clear();
 646             result.push_back(from);
 647             mindistsq2 = dsq1;
 648         }
 649         else if ( are_near(mindistsq2, dsq) )
 650         {
 651             result.push_back(from);
 652         }
 653         if ( mindistsq2 > dsq2 )
 654         {
 655             result.clear();
 656             result.push_back(to);
 657         }
 658         else if ( are_near(mindistsq2, dsq2) )
 659         {
 660             result.push_back(to);
 661         }
 662
 663     }
 664     else
 665     {
 666         if ( result.empty() )
 667         {
 668             if ( are_near(dsq1, dsq2) )
 669             {
 670                 result.push_back(from);
 671                 result.push_back(to);
 672             }
 673             else if ( dsq2 > dsq1 )
 674             {
 675                 result.push_back(from);
 676             }
 677             else
 678             {
 679                 result.push_back(to);
 680             }
 681         }
 682     }
 683
 684     return result;
 685 }
 686
 687
 688 /*
 689  * NOTE: this implementation follows Standard SVG 1.1 implementation guidelines
 690  * for elliptical arc curves. See Appendix F.6.
 691  */
 692 void SVGEllipticalArc::calculate_center_and_extreme_angles()
 693 {
 694     Point d = initialPoint() - finalPoint();
 695
 696     if (is_svg_compliant())
 697     {
 698         if ( initialPoint() == finalPoint() )
 699         {
 700             m_rx = m_ry = m_rot_angle = m_start_angle = m_end_angle = 0;
 701             m_center = initialPoint();
 702             m_large_arc = m_sweep = false;
 703             return;
 704         }
 705
 706         m_rx = std::fabs(m_rx);
 707         m_ry = std::fabs(m_ry);
 708
 709         if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
 710         {
 711             m_rx = L2(d) / 2;
 712             m_ry = 0;
 713             m_rot_angle = std::atan2(d[Y], d[X]);
 714             if (m_rot_angle < 0) m_rot_angle += 2*M_PI;
 715             m_start_angle = 0;
 716             m_end_angle = M_PI;
 717             m_center = middle_point(initialPoint(), finalPoint());
 718             m_large_arc = false;
 719             m_sweep = false;
 720             return;
 721         }
 722     }
 723     else
 724     {
 725         if ( are_near(initialPoint(), finalPoint()) )
 726         {
 727             if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
 728             {
 729                 m_start_angle = m_end_angle = 0;
 730                 m_center = initialPoint();
 731                 return;
 732             }
 733             else
 734             {
 735                 THROW_RANGEERROR("initial and final point are the same");
 736             }
 737         }
 738         if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
 739         { // but initialPoint != finalPoint
 740             THROW_RANGEERROR(
 741                 "there is no ellipse that satisfies the given constraints: "
 742                 "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"
 743             );
 744         }
 745         if ( are_near(ray(Y), 0) )
 746         {
 747             Point v = initialPoint() - finalPoint();
 748             if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )
 749             {
 750                 double angle = std::atan2(v[Y], v[X]);
 751                 if (angle < 0) angle += 2*M_PI;
 752                 if ( are_near( angle, rotation_angle() ) )
 753                 {
 754                     m_start_angle = 0;
 755                     m_end_angle = M_PI;
 756                     m_center = v/2 + finalPoint();
 757                     return;
 758                 }
 759                 angle -= M_PI;
 760                 if ( angle < 0 ) angle += 2*M_PI;
 761                 if ( are_near( angle, rotation_angle() ) )
 762                 {
 763                     m_start_angle = M_PI;
 764                     m_end_angle = 0;
 765                     m_center = v/2 + finalPoint();
 766                     return;
 767                 }
 768                 THROW_RANGEERROR(
 769                     "there is no ellipse that satisfies the given constraints: "
 770                     "ray(Y) == 0 "
 771                     "and slope(initialPoint - finalPoint) != rotation_angle "
 772                     "and != rotation_angle + PI"
 773                 );
 774             }
 775             if ( L2sq(v) > 4*ray(X)*ray(X) )
 776             {
 777                 THROW_RANGEERROR(
 778                     "there is no ellipse that satisfies the given constraints: "
 779                     "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"
 780                 );
 781             }
 782             else
 783             {
 784                 THROW_RANGEERROR(
 785                     "there is infinite ellipses that satisfy the given constraints: "
 786                     "ray(Y) == 0  and distance(initialPoint, finalPoint) < 2*ray(X)"
 787                 );
 788             }
 789
 790         }
 791
 792         if ( are_near(ray(X), 0) )
 793         {
 794             Point v = initialPoint() - finalPoint();
 795             if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )
 796             {
 797                 double angle = std::atan2(v[Y], v[X]);
 798                 if (angle < 0) angle += 2*M_PI;
 799                 double rot_angle = rotation_angle() + M_PI/2;
 800                 if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;
 801                 if ( are_near( angle, rot_angle ) )
 802                 {
 803                     m_start_angle = M_PI/2;
 804                     m_end_angle = 3*M_PI/2;
 805                     m_center = v/2 + finalPoint();
 806                     return;
 807                 }
 808                 angle -= M_PI;
 809                 if ( angle < 0 ) angle += 2*M_PI;
 810                 if ( are_near( angle, rot_angle ) )
 811                 {
 812                     m_start_angle = 3*M_PI/2;
 813                     m_end_angle = M_PI/2;
 814                     m_center = v/2 + finalPoint();
 815                     return;
 816                 }
 817                 THROW_RANGEERROR(
 818                     "there is no ellipse that satisfies the given constraints: "
 819                     "ray(X) == 0 "
 820                     "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "
 821                     "and != rotation_angle + (3/2)*PI"
 822                 );
 823             }
 824             if ( L2sq(v) > 4*ray(Y)*ray(Y) )
 825             {
 826                 THROW_RANGEERROR(
 827                     "there is no ellipse that satisfies the given constraints: "
 828                     "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"
 829                 );
 830             }
 831             else
 832             {
 833                 THROW_RANGEERROR(
 834                     "there is infinite ellipses that satisfy the given constraints: "
 835                     "ray(X) == 0  and distance(initialPoint, finalPoint) < 2*ray(Y)"
 836                 );
 837             }
 838
 839         }
 840
 841     }
 842
 843     double sin_rot_angle = std::sin(rotation_angle());
 844     double cos_rot_angle = std::cos(rotation_angle());
 845
 846
 847     Matrix m( cos_rot_angle, -sin_rot_angle,
 848               sin_rot_angle, cos_rot_angle,
 849               0,             0              );
 850
 851     Point p = (d / 2) * m;
 852     double rx2 = m_rx * m_rx;
 853     double ry2 = m_ry * m_ry;
 854     double rxpy = m_rx * p[Y];
 855     double rypx = m_ry * p[X];
 856     double rx2py2 = rxpy * rxpy;
 857     double ry2px2 = rypx * rypx;
 858     double num = rx2 * ry2;
 859     double den = rx2py2 + ry2px2;
 860     assert(den != 0);
 861     double rad = num / den;
 862     Point c(0,0);
 863     if (rad > 1)
 864     {
 865         rad -= 1;
 866         rad = std::sqrt(rad);
 867
 868         if (m_large_arc == m_sweep) rad = -rad;
 869         c = rad * Point(rxpy / m_ry, -rypx / m_rx);
 870
 871         m[1] = -m[1];
 872         m[2] = -m[2];
 873
 874         m_center = c * m + middle_point(initialPoint(), finalPoint());
 875     }
 876     else if (rad == 1 || is_svg_compliant())
 877     {
 878         double lamda = std::sqrt(1 / rad);
 879         m_rx *= lamda;
 880         m_ry *= lamda;
 881         m_center = middle_point(initialPoint(), finalPoint());
 882     }
 883     else
 884     {
 885         THROW_RANGEERROR(
 886             "there is no ellipse that satisfies the given constraints"
 887         );
 888     }
 889
 890     Point sp((p[X] - c[X]) / m_rx, (p[Y] - c[Y]) / m_ry);
 891     Point ep((-p[X] - c[X]) / m_rx, (-p[Y] - c[Y]) / m_ry);
 892     Point v(1, 0);
 893     m_start_angle = angle_between(v, sp);
 894     double sweep_angle = angle_between(sp, ep);
 895     if (!m_sweep && sweep_angle > 0) sweep_angle -= 2*M_PI;
 896     if (m_sweep && sweep_angle < 0) sweep_angle += 2*M_PI;
 897
 898     if (m_start_angle < 0) m_start_angle += 2*M_PI;
 899     m_end_angle = m_start_angle + sweep_angle;
 900     if (m_end_angle < 0) m_end_angle += 2*M_PI;
 901     if (m_end_angle >= 2*M_PI) m_end_angle -= 2*M_PI;
 902 }
 903
 904
 905 D2<SBasis> SVGEllipticalArc::toSBasis() const
 906 {
 907
 908     if (isDegenerate() && is_svg_compliant())
 909         return chord().toSBasis();
 910
 911     D2<SBasis> arc;
 912     // the interval of parametrization has to be [0,1]
 913     Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );
 914     Linear param(start_angle(), et);
 915     Coord cos_rot_angle = std::cos(rotation_angle());
 916     Coord sin_rot_angle = std::sin(rotation_angle());
 917     // order = 4 seems to be enough to get a perfect looking elliptical arc
 918     SBasis arc_x = ray(X) * cos(param,4);
 919     SBasis arc_y = ray(Y) * sin(param,4);
 920     arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
 921     arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));
 922
 923     // ensure that endpoints remain exact
 924     for ( int d = 0 ; d < 2 ; d++ ) {
 925         arc[d][0][0] = initialPoint()[d];
 926         arc[d][0][1] = finalPoint()[d];
 927     }
 928
 929     return arc;
 930 }
 931
 932
 933 Curve* SVGEllipticalArc::transformed(Matrix const& m) const
 934 {
 935     Ellipse e(center(X), center(Y), ray(X), ray(Y), rotation_angle());
 936     Ellipse et = e.transformed(m);
 937     Point inner_point = pointAt(0.5);
 938     SVGEllipticalArc ea = et.arc( initialPoint() * m,
 939                                   inner_point * m,
 940                                   finalPoint() * m,
 941                                   is_svg_compliant() );
 942     return ea.duplicate();
 943 }
 944
 945 /*
 946  * helper routine to convert the parameter t value
 947  * btw [0,1] and [0,2PI] domain and back
 948  *
 949  */
 950 Coord SVGEllipticalArc::map_to_02PI(Coord t) const
 951 {
 952     Coord angle = start_angle();
 953     if ( sweep_flag() )
 954     {
 955         angle += sweep_angle() * t;
 956     }
 957     else
 958     {
 959         angle -= sweep_angle() * t;
 960     }
 961     angle = std::fmod(angle, 2*M_PI);
 962     if ( angle < 0 ) angle += 2*M_PI;
 963     return angle;
 964 }
 965
 966 Coord SVGEllipticalArc::map_to_01(Coord angle) const
 967 {
 968     return map_circular_arc_on_unit_interval(angle, start_angle(),
 969                                              end_angle(), sweep_flag());
 970 }
 971
 972
 973 namespace detail
 974 {
 975 /*
 976  * ellipse_equation
 977  *
 978  * this is an helper struct, it provides two routines:
 979  * the first one evaluates the implicit form of an ellipse on a given point
 980  * the second one computes the normal versor at a given point of an ellipse
 981  * in implicit form
 982  */
 983 struct ellipse_equation
 984 {
 985     ellipse_equation(double a, double b, double c, double d, double e, double f)
 986         : A(a), B(b), C(c), D(d), E(e), F(f)
 987     {
 988     }
 989
 990     double operator()(double x, double y) const
 991     {
 992         // A * x * x + B * x * y + C * y * y + D * x + E * y + F;
 993         return (A * x + B * y + D) * x + (C * y + E) * y + F;
 994     }
 995
 996     double operator()(Point const& p) const
 997     {
 998         return (*this)(p[X], p[Y]);
 999     }
1000
1001     Point normal(double x, double y) const
1002     {
1003         Point n( 2 * A * x + B * y + D, 2 * C * y + B * x + E );
1004         return unit_vector(n);
1005     }
1006
1007     Point normal(Point const& p) const
1008     {
1009         return normal(p[X], p[Y]);
1010     }
1011
1012     double A, B, C, D, E, F;
1013 };
1014
1015 }
1016
1017 make_elliptical_arc::
1018 make_elliptical_arc( SVGEllipticalArc& _ea,
1019                      curve_type const& _curve,
1020                      unsigned int _total_samples,
1021                      double _tolerance )
1022     : ea(_ea), curve(_curve),
1023       dcurve( unitVector(derivative(curve)) ),
1024       model(), fitter(model, _total_samples),
1025       tolerance(_tolerance), tol_at_extr(tolerance/2),
1026       tol_at_center(0.1), angle_tol(0.1),
1027       initial_point(curve.at0()), final_point(curve.at1()),
1028       N(_total_samples), last(N-1), partitions(N-1), p(N),
1029       svg_compliant(true)
1030 {
1031 }
1032
1033 /*
1034  * check that the coefficients computed by the fit method satisfy
1035  * the tolerance parameters at the k-th sample point
1036  */
1037 bool
1038 make_elliptical_arc::
1039 bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,
1040                 double e1x, double e1y, double e2 )
1041 {
1042     dist_err = std::fabs( ee(p[k]) );
1043     dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 );
1044     // check that the angle btw the tangent versor to the input curve
1045     // and the normal versor of the elliptical arc, both evaluate
1046     // at the k-th sample point, are really othogonal
1047     angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) );
1048     //angle_err *= angle_err;
1049     return ( dist_err  > dist_bound || angle_err > angle_tol );
1050 }
1051
1052 /*
1053  * check that the coefficients computed by the fit method satisfy
1054  * the tolerance parameters at each sample point
1055  */
1056 bool
1057 make_elliptical_arc::
1058 check_bound(double A, double B, double C, double D, double E, double F)
1059 {
1060     detail::ellipse_equation ee(A, B, C, D, E, F);
1061
1062     // check error magnitude at the end-points
1063     double e1x = (2*A + B) * tol_at_extr;
1064     double e1y = (B + 2*C) * tol_at_extr;
1065     double e2 = ((D + E)  + (A + B + C) * tol_at_extr) * tol_at_extr;
1066     if (bound_exceeded(0, ee, e1x, e1y, e2))
1067     {
1068         print_bound_error(0);
1069         return false;
1070     }
1071     if (bound_exceeded(0, ee, e1x, e1y, e2))
1072     {
1073         print_bound_error(last);
1074         return false;
1075     }
1076
1077     // e1x = derivative((ee(x,y), x) | x->tolerance, y->tolerance
1078     e1x = (2*A + B) * tolerance;
1079     // e1y = derivative((ee(x,y), y) | x->tolerance, y->tolerance
1080     e1y = (B + 2*C) * tolerance;
1081     // e2 = ee(tolerance, tolerance) - F;
1082     e2 = ((D + E)  + (A + B + C) * tolerance) * tolerance;
1083 //  std::cerr << "e1x = " << e1x << std::endl;
1084 //  std::cerr << "e1y = " << e1y << std::endl;
1085 //  std::cerr << "e2 = " << e2 << std::endl;
1086
1087     // check error magnitude at sample points
1088     for ( unsigned int k = 1; k < last; ++k )
1089     {
1090         if ( bound_exceeded(k, ee, e1x, e1y, e2) )
1091         {
1092             print_bound_error(k);
1093             return false;
1094         }
1095     }
1096
1097     return true;
1098 }
1099
1100 /*
1101  * fit
1102  *
1103  * supply the samples to the fitter and compute
1104  * the ellipse implicit equation coefficients
1105  */
1106 void make_elliptical_arc::fit()
1107 {
1108     for (unsigned int k = 0; k < N; ++k)
1109     {
1110         p[k] = curve( k / partitions );
1111         fitter.append(p[k]);
1112     }
1113     fitter.update();
1114
1115     NL::Vector z(N, 0.0);
1116     fitter.result(z);
1117 }
1118
1119 bool make_elliptical_arc::make_elliptiarc()
1120 {
1121     const NL::Vector & coeff = fitter.result();
1122     Ellipse e;
1123     try
1124     {
1125         e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
1126     }
1127     catch(LogicalError exc)
1128     {
1129         return false;
1130     }
1131
1132     Point inner_point = curve(0.5);
1133
1134     if (svg_compliant_flag())
1135     {
1136         ea = e.arc(initial_point, inner_point, final_point);
1137     }
1138     else
1139     {
1140         try
1141         {
1142             ea = e.arc(initial_point, inner_point, final_point, false);
1143         }
1144         catch(RangeError exc)
1145         {
1146             return false;
1147         }
1148     }
1149
1150
1151     if ( !are_near( e.center(),
1152                     ea.center(),
1153                     tol_at_center * std::min(e.ray(X),e.ray(Y))
1154                   )
1155        )
1156     {
1157         return false;
1158     }
1159     return true;
1160 }
1161
1162
1163
1164 } // end namespace Geom
1165
1166
1167
1168
1169 /*
1170   Local Variables:
1171   mode:c++
1172   c-file-style:"stroustrup"
1173   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1174   indent-tabs-mode:nil
1175   fill-column:99
1176   End:
1177 */
1178 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
1179