![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
index 6583c948f2b16a10b4e53f4eea99717dcffe0240..b3420fba6c9426ad8f6adff23653b1ec7a01a38f 100644 (file)
-/*
- * SVG Elliptical Arc Class
- *
- * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- */
-
-
-#include "path.h"
-
-
-namespace Geom
-{
-
-D2<SBasis> SVGEllipticalArc::toSBasis() const
-{
- // the interval of parametrization has to be [0,1]
- Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );
- Linear param(start_angle(), et);
- Coord cos_rot_angle = std::cos(rotation_angle());
- Coord sin_rot_angle = std::sin(rotation_angle());
- // order = 4 seems to be enough to get a perfect looking elliptical arc
- // should it be choosen in function of the arc length anyway ?
- // or maybe a user settable parameter: toSBasis(unsigned int order) ?
- SBasis arc_x = ray(X) * cos(param,4);
- SBasis arc_y = ray(Y) * sin(param,4);
- D2<SBasis> arc;
- arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
- arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));
- return arc;
-}
-
-double SVGEllipticalArc::valueAt(Coord t, Dim2 d) const
-{
- Coord tt = from_01_to_02PI(t);
- double sin_rot_angle = std::sin(rotation_angle());
- double cos_rot_angle = std::cos(rotation_angle());
- if ( d == X )
- {
- return ray(X) * cos_rot_angle * std::cos(tt)
- - ray(Y) * sin_rot_angle * std::sin(tt)
- + center(X);
- }
- else
- {
- return ray(X) * sin_rot_angle * std::cos(tt)
- + ray(Y) * cos_rot_angle * std::sin(tt)
- + center(Y);
- }
-}
-
-
-Curve* SVGEllipticalArc::portion(double f, double t) const
-{
- if (f < 0) f = 0;
- if (f > 1) f = 1;
- if (t < 0) t = 0;
- if (t > 1) t = 1;
- SVGEllipticalArc* arc = new SVGEllipticalArc( *this );
- arc->m_initial_point = pointAt(f);
- arc->m_final_point = pointAt(t);
- double sa = sweep_flag() ? sweep_angle() : -sweep_angle();
- arc->m_start_angle = m_start_angle + sa * f;
- if ( !(arc->m_start_angle < 2*M_PI) )
- arc->m_start_angle -= 2*M_PI;
- if ( !(arc->m_start_angle > 0) )
- arc->m_start_angle += 2*M_PI;
- arc->m_end_angle = m_start_angle + sa * t;
- if ( !(arc->m_end_angle < 2*M_PI) )
- arc->m_end_angle -= 2*M_PI;
- if ( !(arc->m_end_angle > 0) )
- arc->m_end_angle += 2*M_PI;
- if ( f > t ) arc->m_sweep = !sweep_flag();
- if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )
- arc->m_large_arc = false;
- return arc;
-}
-
-// NOTE: doesn't work with 360 deg arcs
-void SVGEllipticalArc::calculate_center_and_extreme_angles() throw(RangeError)
-{
- double sin_rot_angle = std::sin(rotation_angle());
- double cos_rot_angle = std::cos(rotation_angle());
-
- Point sp = sweep_flag() ? initialPoint() : finalPoint();
- Point ep = sweep_flag() ? finalPoint() : initialPoint();
-
- Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,
- -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,
- 0, 0 );
- Matrix im = m.inverse();
- Point sol = (ep - sp) * im;
- double half_sum_angle = std::atan2(-sol[X], sol[Y]);
- double half_diff_angle;
- if ( are_near(std::fabs(half_sum_angle), M_PI/2) )
- {
- double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1;
- double arg = anti_sgn_hsa * sol[X] / 2;
- // if |arg| is a little bit > 1 acos returns nan
- if ( are_near(arg, 1) )
- half_diff_angle = 0;
- else if ( are_near(arg, -1) )
- half_diff_angle = M_PI;
- else
- {
- if ( !(-1 < arg && arg < 1) )
- throwRangeError("there is no ellipse that satisfies the given "
- "constraints");
- // assert( -1 < arg && arg < 1 );
- // if it fails
- // => there is no ellipse that satisfies the given constraints
- half_diff_angle = std::acos( arg );
- }
-
- half_diff_angle = M_PI/2 - half_diff_angle;
- }
- else
- {
- double arg = sol[Y] / ( 2 * std::cos(half_sum_angle) );
- // if |arg| is a little bit > 1 asin returns nan
- if ( are_near(arg, 1) )
- half_diff_angle = M_PI/2;
- else if ( are_near(arg, -1) )
- half_diff_angle = -M_PI/2;
- else
- {
- if ( !(-1 < arg && arg < 1) )
- throwRangeError("there is no ellipse that satisfies the given "
- "constraints");
- // assert( -1 < arg && arg < 1 );
- // if it fails
- // => there is no ellipse that satisfies the given constraints
- half_diff_angle = std::asin( arg );
- }
- }
-
- if ( ( m_large_arc && half_diff_angle > 0 )
- || (!m_large_arc && half_diff_angle < 0 ) )
- {
- half_diff_angle = -half_diff_angle;
- }
- if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI;
- if ( half_diff_angle < 0 ) half_diff_angle += M_PI;
-
- m_start_angle = half_sum_angle - half_diff_angle;
- m_end_angle = half_sum_angle + half_diff_angle;
- // 0 <= m_start_angle, m_end_angle < 2PI
- if ( m_start_angle < 0 ) m_start_angle += 2*M_PI;
- if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI;
- sol[0] = std::cos(m_start_angle);
- sol[1] = std::sin(m_start_angle);
- m_center = sp - sol * m;
- if ( !sweep_flag() )
- {
- double angle = m_start_angle;
- m_start_angle = m_end_angle;
- m_end_angle = angle;
- }
-}
-
-Coord SVGEllipticalArc::from_01_to_02PI(Coord t) const
-{
- if ( sweep_flag() )
- {
- Coord angle = start_angle() + sweep_angle() * t;
- if ( !(angle < 2*M_PI) )
- angle -= 2*M_PI;
- return angle;
- }
- else
- {
- Coord angle = start_angle() - sweep_angle() * t;
- if ( angle < 0 ) angle += 2*M_PI;
- return angle;
- }
-}
-
-
-} // end namespace Geom
-
-
-/*
- Local Variables:
- mode:c++
- c-file-style:"stroustrup"
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
- indent-tabs-mode:nil
- fill-column:99
- End:
-*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
-
-
+/*\r
+ * SVG Elliptical Arc Class\r
+ *\r
+ * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#include "path.h"\r
+#include "angle.h"\r
+\r
+#include <gsl/gsl_poly.h>\r
+\r
+#include <cfloat>\r
+\r
+\r
+\r
+\r
+namespace Geom\r
+{\r
+\r
+\r
+Rect SVGEllipticalArc::boundsExact() const\r
+{\r
+ std::vector<double> extremes(4);\r
+ double cosrot = std::cos(rotation_angle());\r
+ double sinrot = std::sin(rotation_angle());\r
+ extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );\r
+ extremes[1] = extremes[0] + M_PI;\r
+ if ( extremes[0] < 0 ) extremes[0] += 2*M_PI; \r
+ extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );\r
+ extremes[3] = extremes[2] + M_PI;\r
+ if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;\r
+ \r
+ \r
+ std::vector<double>arc_extremes(4);\r
+ arc_extremes[0] = initialPoint()[X];\r
+ arc_extremes[1] = finalPoint()[X];\r
+ if ( arc_extremes[0] < arc_extremes[1] ) \r
+ std::swap(arc_extremes[0], arc_extremes[1]);\r
+ arc_extremes[2] = initialPoint()[Y];\r
+ arc_extremes[3] = finalPoint()[Y];\r
+ if ( arc_extremes[2] < arc_extremes[3] ) \r
+ std::swap(arc_extremes[2], arc_extremes[3]);\r
+ \r
+ \r
+ if ( start_angle() < end_angle() )\r
+ {\r
+ if ( sweep_flag() )\r
+ {\r
+ for ( unsigned int i = 0; i < extremes.size(); ++i )\r
+ {\r
+ if ( start_angle() < extremes[i] && extremes[i] < end_angle() )\r
+ {\r
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ for ( unsigned int i = 0; i < extremes.size(); ++i )\r
+ {\r
+ if ( start_angle() > extremes[i] || extremes[i] > end_angle() )\r
+ {\r
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
+ }\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if ( sweep_flag() )\r
+ {\r
+ for ( unsigned int i = 0; i < extremes.size(); ++i )\r
+ {\r
+ if ( start_angle() < extremes[i] || extremes[i] < end_angle() )\r
+ {\r
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ for ( unsigned int i = 0; i < extremes.size(); ++i )\r
+ {\r
+ if ( start_angle() > extremes[i] && extremes[i] > end_angle() )\r
+ {\r
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
+ }\r
+ } \r
+ }\r
+ }\r
+ \r
+ return Rect( Point(arc_extremes[1], arc_extremes[3]) , \r
+ Point(arc_extremes[0], arc_extremes[2]) );\r
+\r
+}\r
+\r
+\r
+std::vector<double> \r
+SVGEllipticalArc::roots(double v, Dim2 d) const\r
+{\r
+ if ( d > Y )\r
+ {\r
+ THROW_RANGEERROR("dimention out of range");\r
+ }\r
+ \r
+ std::vector<double> sol;\r
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
+ {\r
+ if ( center(d) == v )\r
+ sol.push_back(0);\r
+ return sol;\r
+ }\r
+ \r
+ const char* msg[2][2] = \r
+ {\r
+ { "d == X; ray(X) == 0; "\r
+ "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "\r
+ "s should be contained in [-1,1]",\r
+ "d == X; ray(Y) == 0; "\r
+ "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "\r
+ "s should be contained in [-1,1]"\r
+ },\r
+ { "d == Y; ray(X) == 0; "\r
+ "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "\r
+ "s should be contained in [-1,1]",\r
+ "d == Y; ray(Y) == 0; "\r
+ "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "\r
+ "s should be contained in [-1,1]"\r
+ },\r
+ }; \r
+ \r
+ for ( unsigned int dim = 0; dim < 2; ++dim )\r
+ {\r
+ if ( are_near(ray(dim), 0) )\r
+ {\r
+ \r
+ if ( initialPoint()[d] == v && finalPoint()[d] == v )\r
+ {\r
+ THROW_EXCEPTION("infinite solutions");\r
+ }\r
+ if ( (initialPoint()[d] < finalPoint()[d])\r
+ && (initialPoint()[d] > v || finalPoint()[d] < v) )\r
+ {\r
+ return sol;\r
+ }\r
+ if ( (initialPoint()[d] > finalPoint()[d])\r
+ && (finalPoint()[d] > v || initialPoint()[d] < v) )\r
+ {\r
+ return sol;\r
+ }\r
+ double ray_prj;\r
+ switch(d)\r
+ {\r
+ case X: \r
+ switch(dim)\r
+ { \r
+ case X: ray_prj = -ray(Y) * std::sin(rotation_angle());\r
+ break;\r
+ case Y: ray_prj = ray(X) * std::cos(rotation_angle());\r
+ break;\r
+ }\r
+ break;\r
+ case Y:\r
+ switch(dim)\r
+ { \r
+ case X: ray_prj = ray(Y) * std::cos(rotation_angle());\r
+ break;\r
+ case Y: ray_prj = ray(X) * std::sin(rotation_angle());\r
+ break;\r
+ }\r
+ break;\r
+ }\r
+ \r
+ double s = (v - center(d)) / ray_prj;\r
+ if ( s < -1 || s > 1 )\r
+ {\r
+ THROW_LOGICALERROR(msg[d][dim]);\r
+ }\r
+ switch(dim)\r
+ { \r
+ case X: \r
+ s = std::asin(s); // return a value in [-PI/2,PI/2]\r
+ if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) ) )\r
+ {\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ }\r
+ else\r
+ {\r
+ s = M_PI - s;\r
+ if (!(s < 2*M_PI) ) s -= 2*M_PI;\r
+ }\r
+ break;\r
+ case Y: \r
+ s = std::acos(s); // return a value in [0,PI]\r
+ if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )\r
+ {\r
+ s = 2*M_PI - s;\r
+ if ( !(s < 2*M_PI) ) s -= 2*M_PI;\r
+ }\r
+ break;\r
+ }\r
+ \r
+ //std::cerr << "s = " << rad_to_deg(s);\r
+ s = map_to_01(s);\r
+ //std::cerr << " -> t: " << s << std::endl;\r
+ if ( !(s < 0 || s > 1) )\r
+ sol.push_back(s);\r
+ return sol;\r
+ }\r
+ }\r
+ \r
+ double rotx, roty;\r
+ switch(d)\r
+ {\r
+ case X: \r
+ rotx = std::cos(rotation_angle());\r
+ roty = -std::sin(rotation_angle());\r
+ break;\r
+ case Y:\r
+ rotx = std::sin(rotation_angle());\r
+ roty = std::cos(rotation_angle());\r
+ break;\r
+ }\r
+ double rxrotx = ray(X) * rotx;\r
+ double c_v = center(d) - v;\r
+\r
+ double a = -rxrotx + c_v;\r
+ double b = ray(Y) * roty;\r
+ double c = rxrotx + c_v;\r
+ //std::cerr << "a = " << a << std::endl;\r
+ //std::cerr << "b = " << b << std::endl;\r
+ //std::cerr << "c = " << c << std::endl;\r
+ \r
+ if ( are_near(a,0) )\r
+ {\r
+ sol.push_back(M_PI);\r
+ if ( !are_near(b,0) )\r
+ {\r
+ double s = 2 * std::atan(-c/(2*b));\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ sol.push_back(s);\r
+ }\r
+ }\r
+ else\r
+ {\r
+ double delta = b * b - a * c;\r
+ //std::cerr << "delta = " << delta << std::endl;\r
+ if ( are_near(delta, 0) )\r
+ {\r
+ double s = 2 * std::atan(-b/a);\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ sol.push_back(s);\r
+ }\r
+ else if ( delta > 0 )\r
+ { \r
+ double sq = std::sqrt(delta);\r
+ double s = 2 * std::atan( (-b - sq) / a );\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ sol.push_back(s);\r
+ s = 2 * std::atan( (-b + sq) / a );\r
+ if ( s < 0 ) s += 2*M_PI;\r
+ sol.push_back(s);\r
+ }\r
+ }\r
+ \r
+ std::vector<double> arc_sol;\r
+ for (unsigned int i = 0; i < sol.size(); ++i )\r
+ {\r
+ //std::cerr << "s = " << rad_to_deg(sol[i]);\r
+ sol[i] = map_to_01(sol[i]);\r
+ //std::cerr << " -> t: " << sol[i] << std::endl;\r
+ if ( !(sol[i] < 0 || sol[i] > 1) )\r
+ arc_sol.push_back(sol[i]);\r
+ }\r
+ return arc_sol;\r
+ \r
+ \r
+// return SBasisCurve(toSBasis()).roots(v, d);\r
+}\r
+\r
+// D(E(t,C),t) = E(t+PI/2,O)\r
+Curve* SVGEllipticalArc::derivative() const\r
+{\r
+ SVGEllipticalArc* result = new SVGEllipticalArc(*this);\r
+ result->m_center[X] = result->m_center[Y] = 0;\r
+ result->m_start_angle += M_PI/2;\r
+ if( !( result->m_start_angle < 2*M_PI ) )\r
+ {\r
+ result->m_start_angle -= 2*M_PI;\r
+ }\r
+ result->m_end_angle += M_PI/2;\r
+ if( !( result->m_end_angle < 2*M_PI ) )\r
+ {\r
+ result->m_end_angle -= 2*M_PI;\r
+ }\r
+ result->m_initial_point = result->pointAtAngle( result->start_angle() );\r
+ result->m_final_point = result->pointAtAngle( result->end_angle() );\r
+ return result;\r
+ \r
+}\r
+\r
+std::vector<Point> \r
+SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const\r
+{\r
+ std::vector<Point> result;\r
+ result.reserve(n);\r
+ double angle = map_unit_interval_on_circular_arc(t, start_angle(), \r
+ end_angle(), sweep_flag());\r
+ SVGEllipticalArc ea(*this);\r
+ ea.m_center = Point(0,0);\r
+ unsigned int m = std::min(n, 4u);\r
+ for ( unsigned int i = 0; i < m; ++i )\r
+ {\r
+ result.push_back( ea.pointAtAngle(angle) );\r
+ angle += M_PI/2;\r
+ if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;\r
+ }\r
+ m = n / 4;\r
+ for ( unsigned int i = 1; i < m; ++i )\r
+ {\r
+ for ( unsigned int j = 0; j < 4; ++j )\r
+ result.push_back( result[j] );\r
+ }\r
+ m = n - 4 * m;\r
+ for ( unsigned int i = 0; i < m; ++i )\r
+ {\r
+ result.push_back( result[i] );\r
+ }\r
+ if ( !result.empty() ) // n != 0\r
+ result[0] = pointAtAngle(angle);\r
+ return result;\r
+}\r
+\r
+D2<SBasis> SVGEllipticalArc::toSBasis() const\r
+{\r
+ // the interval of parametrization has to be [0,1]\r
+ Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );\r
+ Linear param(start_angle(), et);\r
+ Coord cos_rot_angle = std::cos(rotation_angle());\r
+ Coord sin_rot_angle = std::sin(rotation_angle());\r
+ // order = 4 seems to be enough to get a perfect looking elliptical arc\r
+ // should it be choosen in function of the arc length anyway ?\r
+ // or maybe a user settable parameter: toSBasis(unsigned int order) ?\r
+ SBasis arc_x = ray(X) * cos(param,4);\r
+ SBasis arc_y = ray(Y) * sin(param,4);\r
+ D2<SBasis> arc;\r
+ arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));\r
+ arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));\r
+ return arc;\r
+}\r
+\r
+\r
+bool SVGEllipticalArc::containsAngle(Coord angle) const\r
+{\r
+ if ( sweep_flag() )\r
+ if ( start_angle() < end_angle() )\r
+ if ( !( angle < start_angle() || angle > end_angle() ) )\r
+ return true;\r
+ else\r
+ return false;\r
+ else\r
+ if ( !( angle < start_angle() && angle > end_angle() ) )\r
+ return true;\r
+ else\r
+ return false;\r
+ else\r
+ if ( start_angle() > end_angle() )\r
+ if ( !( angle > start_angle() || angle < end_angle() ) )\r
+ return true;\r
+ else\r
+ return false;\r
+ else\r
+ if ( !( angle > start_angle() && angle < end_angle() ) )\r
+ return true;\r
+ else\r
+ return false;\r
+}\r
+\r
+\r
+double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const\r
+{\r
+ double sin_rot_angle = std::sin(rotation_angle());\r
+ double cos_rot_angle = std::cos(rotation_angle());\r
+ if ( d == X )\r
+ {\r
+ return ray(X) * cos_rot_angle * std::cos(t) \r
+ - ray(Y) * sin_rot_angle * std::sin(t) \r
+ + center(X);\r
+ }\r
+ else if ( d == Y )\r
+ {\r
+ return ray(X) * sin_rot_angle * std::cos(t) \r
+ + ray(Y) * cos_rot_angle * std::sin(t) \r
+ + center(Y);\r
+ }\r
+ THROW_RANGEERROR("dimension parameter out of range");\r
+}\r
+\r
+\r
+Curve* SVGEllipticalArc::portion(double f, double t) const \r
+{\r
+ if (f < 0) f = 0;\r
+ if (f > 1) f = 1;\r
+ if (t < 0) t = 0;\r
+ if (t > 1) t = 1;\r
+ if ( are_near(f, t) )\r
+ {\r
+ SVGEllipticalArc* arc = new SVGEllipticalArc();\r
+ arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);\r
+ arc->m_start_angle = arc->m_end_angle = m_start_angle;\r
+ arc->m_rot_angle = m_rot_angle;\r
+ arc->m_sweep = m_sweep;\r
+ arc->m_large_arc = m_large_arc;\r
+ }\r
+ SVGEllipticalArc* arc = new SVGEllipticalArc( *this );\r
+ arc->m_initial_point = pointAt(f);\r
+ arc->m_final_point = pointAt(t);\r
+ double sa = sweep_flag() ? sweep_angle() : -sweep_angle();\r
+ arc->m_start_angle = m_start_angle + sa * f;\r
+ if ( !(arc->m_start_angle < 2*M_PI) )\r
+ arc->m_start_angle -= 2*M_PI;\r
+ if ( arc->m_start_angle < 0 )\r
+ arc->m_start_angle += 2*M_PI;\r
+ arc->m_end_angle = m_start_angle + sa * t;\r
+ if ( !(arc->m_end_angle < 2*M_PI) )\r
+ arc->m_end_angle -= 2*M_PI;\r
+ if ( arc->m_end_angle < 0 )\r
+ arc->m_end_angle += 2*M_PI;\r
+ if ( f > t ) arc->m_sweep = !sweep_flag();\r
+ if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )\r
+ arc->m_large_arc = false;\r
+ return arc;\r
+}\r
+\r
+// NOTE: doesn't work with 360 deg arcs\r
+void SVGEllipticalArc::calculate_center_and_extreme_angles()\r
+{\r
+ if ( are_near(initialPoint(), finalPoint()) )\r
+ {\r
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
+ {\r
+ m_start_angle = m_end_angle = 0;\r
+ m_center = initialPoint();\r
+ return;\r
+ }\r
+ else\r
+ {\r
+ THROW_RANGEERROR("initial and final point are the same");\r
+ }\r
+ }\r
+ if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
+ { // but initialPoint != finalPoint\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"\r
+ );\r
+ }\r
+ if ( are_near(ray(Y), 0) )\r
+ {\r
+ Point v = initialPoint() - finalPoint();\r
+ if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )\r
+ {\r
+ double angle = std::atan2(v[Y], v[X]);\r
+ if (angle < 0) angle += 2*M_PI;\r
+ if ( are_near( angle, rotation_angle() ) )\r
+ {\r
+ m_start_angle = 0;\r
+ m_end_angle = M_PI;\r
+ m_center = v/2 + finalPoint();\r
+ return;\r
+ }\r
+ angle -= M_PI;\r
+ if ( angle < 0 ) angle += 2*M_PI;\r
+ if ( are_near( angle, rotation_angle() ) )\r
+ {\r
+ m_start_angle = M_PI;\r
+ m_end_angle = 0;\r
+ m_center = v/2 + finalPoint();\r
+ return;\r
+ }\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(Y) == 0 "\r
+ "and slope(initialPoint - finalPoint) != rotation_angle "\r
+ "and != rotation_angle + PI"\r
+ );\r
+ }\r
+ if ( L2sq(v) > 4*ray(X)*ray(X) )\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"\r
+ );\r
+ }\r
+ else\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is infinite ellipses that satisfy the given constraints: "\r
+ "ray(Y) == 0 and distance(initialPoint, finalPoint) < 2*ray(X)"\r
+ );\r
+ }\r
+ \r
+ }\r
+ \r
+ if ( are_near(ray(X), 0) )\r
+ {\r
+ Point v = initialPoint() - finalPoint();\r
+ if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )\r
+ {\r
+ double angle = std::atan2(v[Y], v[X]);\r
+ if (angle < 0) angle += 2*M_PI;\r
+ double rot_angle = rotation_angle() + M_PI/2;\r
+ if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;\r
+ if ( are_near( angle, rot_angle ) )\r
+ {\r
+ m_start_angle = M_PI/2;\r
+ m_end_angle = 3*M_PI/2;\r
+ m_center = v/2 + finalPoint();\r
+ return;\r
+ }\r
+ angle -= M_PI;\r
+ if ( angle < 0 ) angle += 2*M_PI;\r
+ if ( are_near( angle, rot_angle ) )\r
+ {\r
+ m_start_angle = 3*M_PI/2;\r
+ m_end_angle = M_PI/2;\r
+ m_center = v/2 + finalPoint();\r
+ return;\r
+ }\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(X) == 0 "\r
+ "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "\r
+ "and != rotation_angle + (3/2)*PI"\r
+ );\r
+ }\r
+ if ( L2sq(v) > 4*ray(Y)*ray(Y) )\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints: "\r
+ "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"\r
+ );\r
+ }\r
+ else\r
+ {\r
+ THROW_RANGEERROR(\r
+ "there is infinite ellipses that satisfy the given constraints: "\r
+ "ray(X) == 0 and distance(initialPoint, finalPoint) < 2*ray(Y)"\r
+ );\r
+ }\r
+ \r
+ }\r
+ \r
+ double sin_rot_angle = std::sin(rotation_angle());\r
+ double cos_rot_angle = std::cos(rotation_angle());\r
+\r
+ Point sp = sweep_flag() ? initialPoint() : finalPoint();\r
+ Point ep = sweep_flag() ? finalPoint() : initialPoint();\r
+\r
+ Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,\r
+ -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,\r
+ 0, 0 );\r
+ Matrix im = m.inverse();\r
+ Point sol = (ep - sp) * im;\r
+ double half_sum_angle = std::atan2(-sol[X], sol[Y]);\r
+ double half_diff_angle;\r
+ if ( are_near(std::fabs(half_sum_angle), M_PI/2) )\r
+ {\r
+ double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1;\r
+ double arg = anti_sgn_hsa * sol[X] / 2;\r
+ // if |arg| is a little bit > 1 acos returns nan\r
+ if ( are_near(arg, 1) )\r
+ half_diff_angle = 0;\r
+ else if ( are_near(arg, -1) )\r
+ half_diff_angle = M_PI;\r
+ else\r
+ {\r
+ if ( !(-1 < arg && arg < 1) )\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints"\r
+ );\r
+ // assert( -1 < arg && arg < 1 );\r
+ // if it fails \r
+ // => there is no ellipse that satisfies the given constraints\r
+ half_diff_angle = std::acos( arg );\r
+ }\r
+\r
+ half_diff_angle = M_PI/2 - half_diff_angle;\r
+ }\r
+ else\r
+ {\r
+ double arg = sol[Y] / ( 2 * std::cos(half_sum_angle) );\r
+ // if |arg| is a little bit > 1 asin returns nan\r
+ if ( are_near(arg, 1) ) \r
+ half_diff_angle = M_PI/2;\r
+ else if ( are_near(arg, -1) )\r
+ half_diff_angle = -M_PI/2;\r
+ else\r
+ {\r
+ if ( !(-1 < arg && arg < 1) )\r
+ THROW_RANGEERROR(\r
+ "there is no ellipse that satisfies the given constraints"\r
+ );\r
+ // assert( -1 < arg && arg < 1 );\r
+ // if it fails \r
+ // => there is no ellipse that satisfies the given constraints\r
+ half_diff_angle = std::asin( arg );\r
+ }\r
+ }\r
+\r
+ if ( ( m_large_arc && half_diff_angle > 0 ) \r
+ || (!m_large_arc && half_diff_angle < 0 ) )\r
+ {\r
+ half_diff_angle = -half_diff_angle;\r
+ }\r
+ if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI;\r
+ if ( half_diff_angle < 0 ) half_diff_angle += M_PI;\r
+ \r
+ m_start_angle = half_sum_angle - half_diff_angle;\r
+ m_end_angle = half_sum_angle + half_diff_angle;\r
+ // 0 <= m_start_angle, m_end_angle < 2PI\r
+ if ( m_start_angle < 0 ) m_start_angle += 2*M_PI;\r
+ if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI;\r
+ sol[0] = std::cos(m_start_angle);\r
+ sol[1] = std::sin(m_start_angle);\r
+ m_center = sp - sol * m;\r
+ if ( !sweep_flag() )\r
+ {\r
+ double angle = m_start_angle;\r
+ m_start_angle = m_end_angle;\r
+ m_end_angle = angle;\r
+ }\r
+}\r
+\r
+Coord SVGEllipticalArc::map_to_02PI(Coord t) const\r
+{\r
+ if ( sweep_flag() )\r
+ {\r
+ Coord angle = start_angle() + sweep_angle() * t;\r
+ if ( !(angle < 2*M_PI) )\r
+ angle -= 2*M_PI;\r
+ return angle;\r
+ }\r
+ else\r
+ {\r
+ Coord angle = start_angle() - sweep_angle() * t;\r
+ if ( angle < 0 ) angle += 2*M_PI;\r
+ return angle;\r
+ }\r
+}\r
+\r
+Coord SVGEllipticalArc::map_to_01(Coord angle) const \r
+{\r
+ return map_circular_arc_on_unit_interval(angle, start_angle(), \r
+ end_angle(), sweep_flag());\r
+}\r
+\r
+\r
+std::vector<double> SVGEllipticalArc::\r
+allNearestPoints( Point const& p, double from, double to ) const\r
+{\r
+ if ( from > to ) std::swap(from, to);\r
+ if ( from < 0 || to > 1 )\r
+ {\r
+ THROW_RANGEERROR("[from,to] interval out of range");\r
+ }\r
+ std::vector<double> result;\r
+ if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )\r
+ {\r
+ result.push_back(from);\r
+ return result;\r
+ }\r
+ else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )\r
+ {\r
+ LineSegment seg(pointAt(from), pointAt(to));\r
+ Point np = seg.pointAt( seg.nearestPoint(p) );\r
+ if ( are_near(ray(Y), 0) )\r
+ {\r
+ if ( are_near(rotation_angle(), M_PI/2) \r
+ || are_near(rotation_angle(), 3*M_PI/2) )\r
+ {\r
+ result = roots(np[Y], Y);\r
+ }\r
+ else\r
+ {\r
+ result = roots(np[X], X);\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if ( are_near(rotation_angle(), M_PI/2) \r
+ || are_near(rotation_angle(), 3*M_PI/2) )\r
+ {\r
+ result = roots(np[X], X);\r
+ }\r
+ else\r
+ {\r
+ result = roots(np[Y], Y);\r
+ }\r
+ }\r
+ return result;\r
+ }\r
+ else if ( are_near(ray(X), ray(Y)) )\r
+ {\r
+ Point r = p - center();\r
+ if ( are_near(r, Point(0,0)) )\r
+ {\r
+ THROW_EXCEPTION("infinite nearest points");\r
+ }\r
+ // TODO: implement case r != 0\r
+// Point np = ray(X) * unit_vector(r);\r
+// std::vector<double> solX = roots(np[X],X);\r
+// std::vector<double> solY = roots(np[Y],Y);\r
+// double t;\r
+// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))\r
+// {\r
+// t = solX[0];\r
+// }\r
+// else\r
+// {\r
+// t = solX[1];\r
+// }\r
+// if ( !(t < from || t > to) )\r
+// {\r
+// result.push_back(t);\r
+// }\r
+// else\r
+// {\r
+// \r
+// }\r
+ }\r
+ \r
+ // solve the equation <D(E(t),t)|E(t)-p> == 0\r
+ // that provides min and max distance points \r
+ // on the ellipse E wrt the point p\r
+ // after the substitutions: \r
+ // cos(t) = (1 - s^2) / (1 + s^2)\r
+ // sin(t) = 2t / (1 + s^2)\r
+ // where s = tan(t/2)\r
+ // we get a 4th degree equation in s\r
+ /*\r
+ * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) + \r
+ * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) + \r
+ * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) + \r
+ * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])\r
+ */\r
+\r
+ Point p_c = p - center();\r
+ double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));\r
+ double cosrot = std::cos( rotation_angle() );\r
+ double sinrot = std::sin( rotation_angle() );\r
+ double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);\r
+ double coeff[5];\r
+ coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );\r
+ coeff[3] = 2 * ( rx2_ry2 + expr1 );\r
+ coeff[2] = 0;\r
+ coeff[1] = 2 * ( -rx2_ry2 + expr1 );\r
+ coeff[0] = -coeff[4];\r
+ \r
+// for ( unsigned int i = 0; i < 5; ++i )\r
+// std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;\r
+ \r
+ std::vector<double> real_sol;\r
+ // gsl_poly_complex_solve raises an error \r
+ // if the leading coefficient is zero\r
+ if ( are_near(coeff[4], 0) ) \r
+ {\r
+ real_sol.push_back(0);\r
+ if ( !are_near(coeff[3], 0) )\r
+ {\r
+ double sq = -coeff[1] / coeff[3];\r
+ if ( sq > 0 )\r
+ {\r
+ double s = std::sqrt(sq);\r
+ real_sol.push_back(s);\r
+ real_sol.push_back(-s);\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ double sol[8];\r
+ gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5);\r
+ gsl_poly_complex_solve(coeff, 5, w, sol );\r
+ gsl_poly_complex_workspace_free(w);\r
+ \r
+ for ( unsigned int i = 0; i < 4; ++i )\r
+ {\r
+ if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]);\r
+ }\r
+ }\r
+ \r
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )\r
+ {\r
+ real_sol[i] = 2 * std::atan(real_sol[i]);\r
+ if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;\r
+ }\r
+ // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0\r
+ // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity\r
+ if ( (real_sol.size() % 2) != 0 )\r
+ {\r
+ real_sol.push_back(M_PI);\r
+ }\r
+ \r
+ double mindistsq1 = std::numeric_limits<double>::max();\r
+ double mindistsq2 = std::numeric_limits<double>::max();\r
+ double dsq;\r
+ unsigned int mi1, mi2;\r
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )\r
+ {\r
+ dsq = distanceSq(p, pointAtAngle(real_sol[i]));\r
+ if ( mindistsq1 > dsq )\r
+ {\r
+ mindistsq2 = mindistsq1;\r
+ mi2 = mi1;\r
+ mindistsq1 = dsq;\r
+ mi1 = i;\r
+ }\r
+ else if ( mindistsq2 > dsq )\r
+ {\r
+ mindistsq2 = dsq;\r
+ mi2 = i;\r
+ }\r
+ }\r
+ \r
+ double t = map_to_01( real_sol[mi1] );\r
+ if ( !(t < from || t > to) )\r
+ {\r
+ result.push_back(t);\r
+ }\r
+ \r
+ bool second_sol = false; \r
+ t = map_to_01( real_sol[mi2] );\r
+ if ( real_sol.size() == 4 && !(t < from || t > to) )\r
+ {\r
+ if ( result.empty() || are_near(mindistsq1, mindistsq2) )\r
+ {\r
+ result.push_back(t);\r
+ second_sol = true;\r
+ }\r
+ }\r
+ \r
+ // we need to test extreme points too\r
+ double dsq1 = distanceSq(p, pointAt(from));\r
+ double dsq2 = distanceSq(p, pointAt(to));\r
+ if ( second_sol )\r
+ {\r
+ if ( mindistsq2 > dsq1 )\r
+ {\r
+ result.clear();\r
+ result.push_back(from);\r
+ mindistsq2 = dsq1;\r
+ }\r
+ else if ( are_near(mindistsq2, dsq) )\r
+ {\r
+ result.push_back(from);\r
+ }\r
+ if ( mindistsq2 > dsq2 )\r
+ {\r
+ result.clear();\r
+ result.push_back(to);\r
+ }\r
+ else if ( are_near(mindistsq2, dsq2) )\r
+ {\r
+ result.push_back(to);\r
+ }\r
+ \r
+ }\r
+ else\r
+ {\r
+ if ( result.empty() )\r
+ {\r
+ if ( are_near(dsq1, dsq2) )\r
+ {\r
+ result.push_back(from);\r
+ result.push_back(to);\r
+ }\r
+ else if ( dsq2 > dsq1 )\r
+ {\r
+ result.push_back(from);\r
+ }\r
+ else\r
+ {\r
+ result.push_back(to);\r
+ }\r
+ }\r
+ }\r
+ \r
+ return result;\r
+}\r
+\r
+\r
+} // end namespace Geom\r
+\r
+\r
+/*\r
+ Local Variables:\r
+ mode:c++\r
+ c-file-style:"stroustrup"\r
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
+ indent-tabs-mode:nil\r
+ fill-column:99\r
+ End:\r
+*/\r
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+\r
+\r