src/2geom/svg-elliptical-arc.cpp

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