diff options
Diffstat (limited to 'src/2geom/elliptical-arc.cpp')
| -rw-r--r-- | src/2geom/elliptical-arc.cpp | 1045 |
1 files changed, 1045 insertions, 0 deletions
diff --git a/src/2geom/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp new file mode 100644 index 0000000..63e534c --- /dev/null +++ b/src/2geom/elliptical-arc.cpp @@ -0,0 +1,1045 @@ +/* + * SVG Elliptical Arc Class + * + * Authors: + * Marco Cecchetti <mrcekets at gmail.com> + * Krzysztof Kosiński <tweenk.pl@gmail.com> + * Copyright 2008-2009 Authors + * + * 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 <cfloat> +#include <limits> +#include <memory> + +#include <2geom/bezier-curve.h> +#include <2geom/ellipse.h> +#include <2geom/elliptical-arc.h> +#include <2geom/path-sink.h> +#include <2geom/sbasis-geometric.h> +#include <2geom/transforms.h> +#include <2geom/utils.h> + +#include <2geom/numeric/vector.h> +#include <2geom/numeric/fitting-tool.h> +#include <2geom/numeric/fitting-model.h> + +namespace Geom +{ + +/** + * @class EllipticalArc + * @brief Elliptical arc curve + * + * Elliptical arc is a curve taking the shape of a section of an ellipse. + * + * The arc function has two forms: the regular one, mapping the unit interval to points + * on 2D plane (the linear domain), and a second form that maps some interval + * \f$A \subseteq [0,2\pi)\f$ to the same points (the angular domain). The interval \f$A\f$ + * determines which part of the ellipse forms the arc. The arc is said to contain an angle + * if its angular domain includes that angle (and therefore it is defined for that angle). + * + * The angular domain considers each ellipse to be + * a rotated, scaled and translated unit circle: 0 corresponds to \f$(1,0)\f$ on the unit circle, + * \f$\pi/2\f$ corresponds to \f$(0,1)\f$, \f$\pi\f$ to \f$(-1,0)\f$ and \f$3\pi/2\f$ + * to \f$(0,-1)\f$. After the angle is mapped to a point from a unit circle, the point is + * transformed using a matrix of this form + * \f[ M = \left[ \begin{array}{ccc} + r_X \cos(\theta) & -r_Y \sin(\theta) & 0 \\ + r_X \sin(\theta) & r_Y \cos(\theta) & 0 \\ + c_X & c_Y & 1 \end{array} \right] \f] + * where \f$r_X, r_Y\f$ are the X and Y rays of the ellipse, \f$\theta\f$ is its angle of rotation, + * and \f$c_X, c_Y\f$ the coordinates of the ellipse's center - thus mapping the angle + * to some point on the ellipse. Note that for example the point at angluar coordinate 0, + * the center and the point at angular coordinate \f$\pi/4\f$ do not necessarily + * create an angle of \f$\pi/4\f$ radians; it is only the case if both axes of the ellipse + * are of the same length (i.e. it is a circle). + * + * @image html ellipse-angular-coordinates.png "An illustration of the angular domain" + * + * Each arc is defined by five variables: The initial and final point, the ellipse's rays, + * and the ellipse's rotation. Each set of those parameters corresponds to four different arcs, + * with two of them larger than half an ellipse and two of them turning clockwise while traveling + * from initial to final point. The two flags disambiguate between them: "large arc flag" selects + * the bigger arc, while the "sweep flag" selects the arc going in the direction of positive + * angles. Angles always increase when going from the +X axis in the direction of the +Y axis, + * so if Y grows downwards, this means clockwise. + * + * @image html elliptical-arc-flags.png "Meaning of arc flags (Y grows downwards)" + * + * @ingroup Curves + */ + + +/** @brief Compute bounds of an elliptical arc. + * The bounds computation works as follows. The extreme X and Y points + * are either the endpoints or local minima / maxima of the ellipse. + * We already have endpoints, and we compute the local extremes. + * The local extremes correspond to two angles separated by \f$\pi\f$. + * Once we compute these angles, we check whether they belong to the arc, + * and if they do, we evaluate the ellipse at these angles. + * The bounding box of the arc is equal to the bounding box of the endpoints + * and the local extrema that belong to the arc. + */ +Rect EllipticalArc::boundsExact() const +{ + if (isChord()) { + return { _initial_point, _final_point }; + } + + if (_angles.isFull()) { + return _ellipse.boundsExact(); + } + + auto const trans = unitCircleTransform(); + + auto proj_bounds = [&] (Dim2 d) { + // The dth coordinate function pulls back to trans[d] * x + trans[d + 2] * y + trans[d + 4] + // in the coordinate system where the ellipse is a unit circle. We compute its range of + // values on the unit circle arc. + auto result = Interval(_initial_point[d], _final_point[d]); + + auto const v = Point(trans[d], trans[d + 2]); + auto const r = v.length(); + auto const mid = trans[d + 4]; + auto const angle = Angle(v); + + if (_angles.contains(angle)) { + result.expandTo(mid + r); + } + if (_angles.contains(angle + M_PI)) { + result.expandTo(mid - r); + } + + return result; + }; + + return { proj_bounds(X), proj_bounds(Y) }; +} + +void EllipticalArc::expandToTransformed(Rect &bbox, Affine const &transform) const +{ + bbox.expandTo(_final_point * transform); + + if (isChord() || bbox.contains(_ellipse.boundsFast())) { + return; + } + + auto const trans = unitCircleTransform() * transform; + + for (auto d : { X, Y }) { + // See boundsExact() for explanation. + auto const v = Point(trans[d], trans[d + 2]); + auto const r = v.length(); + auto const mid = trans[d + 4]; + + if (_angles.isFull()) { + bbox[d].unionWith(Interval(mid - r, mid + r)); + } else { + auto const angle = Angle(v); + if (_angles.contains(angle)) { + bbox[d].expandTo(mid + r); + } + if (_angles.contains(angle + M_PI)) { + bbox[d].expandTo(mid - r); + } + } + } +} + +Point EllipticalArc::pointAtAngle(Coord t) const +{ + Point ret = _ellipse.pointAt(t); + return ret; +} + +Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const +{ + return _ellipse.valueAt(t, d); +} + +std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const +{ + std::vector<Coord> sol; + + if (isChord()) { + sol = chord().roots(v, d); + return sol; + } + + Interval unit_interval(0, 1); + + double rotx, roty; + if (d == X) { + sincos(rotationAngle(), roty, rotx); + roty = -roty; + } else { + sincos(rotationAngle(), rotx, roty); + } + + double rxrotx = ray(X) * rotx; + double c_v = center(d) - v; + + double a = -rxrotx + c_v; + double b = ray(Y) * roty; + double c = rxrotx + c_v; + //std::cerr << "a = " << a << std::endl; + //std::cerr << "b = " << b << std::endl; + //std::cerr << "c = " << c << std::endl; + + if (a == 0) + { + sol.push_back(M_PI); + if (b != 0) + { + double s = 2 * std::atan(-c/(2*b)); + if ( s < 0 ) s += 2*M_PI; + sol.push_back(s); + } + } + else + { + double delta = b * b - a * c; + //std::cerr << "delta = " << delta << std::endl; + if (delta == 0) { + double s = 2 * std::atan(-b/a); + if ( s < 0 ) s += 2*M_PI; + sol.push_back(s); + } + else if ( delta > 0 ) + { + double sq = std::sqrt(delta); + double s = 2 * std::atan( (-b - sq) / a ); + if ( s < 0 ) s += 2*M_PI; + sol.push_back(s); + s = 2 * std::atan( (-b + sq) / a ); + if ( s < 0 ) s += 2*M_PI; + sol.push_back(s); + } + } + + std::vector<double> arc_sol; + for (double & i : sol) { + //std::cerr << "s = " << deg_from_rad(sol[i]); + i = timeAtAngle(i); + //std::cerr << " -> t: " << sol[i] << std::endl; + if (unit_interval.contains(i)) { + arc_sol.push_back(i); + } + } + return arc_sol; +} + + +// D(E(t,C),t) = E(t+PI/2,O), where C is the ellipse center +// the derivative doesn't rotate the ellipse but there is a translation +// of the parameter t by an angle of PI/2 so the ellipse points are shifted +// of such an angle in the cw direction +Curve *EllipticalArc::derivative() const +{ + if (isChord()) { + return chord().derivative(); + } + + EllipticalArc *result = static_cast<EllipticalArc*>(duplicate()); + result->_ellipse.setCenter(0, 0); + result->_angles.setInitial(result->_angles.initialAngle() + M_PI/2); + result->_angles.setFinal(result->_angles.finalAngle() + M_PI/2); + result->_initial_point = result->pointAtAngle( result->initialAngle() ); + result->_final_point = result->pointAtAngle( result->finalAngle() ); + return result; +} + + +std::vector<Point> +EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const +{ + if (isChord()) { + return chord().pointAndDerivatives(t, n); + } + + unsigned int nn = n+1; // nn represents the size of the result vector. + std::vector<Point> result; + result.reserve(nn); + double angle = angleAt(t); + std::unique_ptr<EllipticalArc> ea( static_cast<EllipticalArc*>(duplicate()) ); + ea->_ellipse.setCenter(0, 0); + unsigned int m = std::min(nn, 4u); + for ( unsigned int i = 0; i < m; ++i ) + { + result.push_back( ea->pointAtAngle(angle) ); + angle += (sweep() ? M_PI/2 : -M_PI/2); + if ( !(angle < 2*M_PI) ) angle -= 2*M_PI; + } + m = nn / 4; + for ( unsigned int i = 1; i < m; ++i ) + { + for ( unsigned int j = 0; j < 4; ++j ) + result.push_back( result[j] ); + } + m = nn - 4 * m; + for ( unsigned int i = 0; i < m; ++i ) + { + result.push_back( result[i] ); + } + if ( !result.empty() ) // nn != 0 + result[0] = pointAtAngle(angle); + return result; +} + +Point EllipticalArc::pointAt(Coord t) const +{ + if (t == 0.0) { + return initialPoint(); + } + if (t == 1.0) { + return finalPoint(); + } + if (isChord()) { + return chord().pointAt(t); + } + return _ellipse.pointAt(angleAt(t)); +} + +Coord EllipticalArc::valueAt(Coord t, Dim2 d) const +{ + if (isChord()) return chord().valueAt(t, d); + return valueAtAngle(angleAt(t), d); +} + +Curve* EllipticalArc::portion(double f, double t) const +{ + // fix input arguments + f = std::clamp(f, 0.0, 1.0); + t = std::clamp(t, 0.0, 1.0); + + if (f == t) { + EllipticalArc *arc = new EllipticalArc(); + arc->_initial_point = arc->_final_point = pointAt(f); + return arc; + } + if (f == 0.0 && t == 1.0) { + return duplicate(); + } + if (f == 1.0 && t == 0.0) { + return reverse(); + } + + EllipticalArc *arc = static_cast<EllipticalArc*>(duplicate()); + arc->_initial_point = pointAt(f); + arc->_final_point = pointAt(t); + arc->_angles.setAngles(angleAt(f), angleAt(t)); + if (f > t) arc->_angles.setSweep(!sweep()); + if ( _large_arc && fabs(angularExtent() * (t-f)) <= M_PI) { + arc->_large_arc = false; + } + return arc; +} + +// the arc is the same but traversed in the opposite direction +Curve *EllipticalArc::reverse() const +{ + using std::swap; + EllipticalArc *rarc = static_cast<EllipticalArc*>(duplicate()); + rarc->_angles.reverse(); + swap(rarc->_initial_point, rarc->_final_point); + return rarc; +} + +#ifdef HAVE_GSL // GSL is required for function "solve_reals" +std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from, double to ) const +{ + std::vector<double> result; + + if ( from > to ) std::swap(from, to); + if ( from < 0 || to > 1 ) + { + THROW_RANGEERROR("[from,to] interval out of range"); + } + + if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) ) + { + result.push_back(from); + return result; + } + else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) ) + { + LineSegment seg(pointAt(from), pointAt(to)); + Point np = seg.pointAt( seg.nearestTime(p) ); + if ( are_near(ray(Y), 0) ) + { + if ( are_near(rotationAngle(), M_PI/2) + || are_near(rotationAngle(), 3*M_PI/2) ) + { + result = roots(np[Y], Y); + } + else + { + result = roots(np[X], X); + } + } + else + { + if ( are_near(rotationAngle(), M_PI/2) + || are_near(rotationAngle(), 3*M_PI/2) ) + { + result = roots(np[X], X); + } + else + { + result = roots(np[Y], Y); + } + } + return result; + } + else if ( are_near(ray(X), ray(Y)) ) + { + Point r = p - center(); + if ( are_near(r, Point(0,0)) ) + { + THROW_INFINITESOLUTIONS(0); + } + // TODO: implement case r != 0 +// Point np = ray(X) * unit_vector(r); +// std::vector<double> solX = roots(np[X],X); +// std::vector<double> solY = roots(np[Y],Y); +// double t; +// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1])) +// { +// t = solX[0]; +// } +// else +// { +// t = solX[1]; +// } +// if ( !(t < from || t > to) ) +// { +// result.push_back(t); +// } +// else +// { +// +// } + } + + // solve the equation <D(E(t),t)|E(t)-p> == 0 + // that provides min and max distance points + // on the ellipse E wrt the point p + // after the substitutions: + // cos(t) = (1 - s^2) / (1 + s^2) + // sin(t) = 2t / (1 + s^2) + // where s = tan(t/2) + // we get a 4th degree equation in s + /* + * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) + + * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) + + * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) + + * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) + */ + + Point p_c = p - center(); + double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y)); + double sinrot, cosrot; + sincos(rotationAngle(), sinrot, cosrot); + double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot); + Poly coeff; + coeff.resize(5); + coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot ); + coeff[3] = 2 * ( rx2_ry2 + expr1 ); + coeff[2] = 0; + coeff[1] = 2 * ( -rx2_ry2 + expr1 ); + coeff[0] = -coeff[4]; + +// for ( unsigned int i = 0; i < 5; ++i ) +// std::cerr << "c[" << i << "] = " << coeff[i] << std::endl; + + std::vector<double> real_sol; + // gsl_poly_complex_solve raises an error + // if the leading coefficient is zero + if ( are_near(coeff[4], 0) ) + { + real_sol.push_back(0); + if ( !are_near(coeff[3], 0) ) + { + double sq = -coeff[1] / coeff[3]; + if ( sq > 0 ) + { + double s = std::sqrt(sq); + real_sol.push_back(s); + real_sol.push_back(-s); + } + } + } + else + { + real_sol = solve_reals(coeff); + } + + for (double & i : real_sol) + { + i = 2 * std::atan(i); + if ( i < 0 ) i += 2*M_PI; + } + // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0 + // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity + if ( (real_sol.size() % 2) != 0 ) + { + real_sol.push_back(M_PI); + } + + double mindistsq1 = std::numeric_limits<double>::max(); + double mindistsq2 = std::numeric_limits<double>::max(); + double dsq = 0; + unsigned int mi1 = 0, mi2 = 0; + for ( unsigned int i = 0; i < real_sol.size(); ++i ) + { + dsq = distanceSq(p, pointAtAngle(real_sol[i])); + if ( mindistsq1 > dsq ) + { + mindistsq2 = mindistsq1; + mi2 = mi1; + mindistsq1 = dsq; + mi1 = i; + } + else if ( mindistsq2 > dsq ) + { + mindistsq2 = dsq; + mi2 = i; + } + } + + double t = timeAtAngle(real_sol[mi1]); + if ( !(t < from || t > to) ) + { + result.push_back(t); + } + + bool second_sol = false; + t = timeAtAngle(real_sol[mi2]); + if ( real_sol.size() == 4 && !(t < from || t > to) ) + { + if ( result.empty() || are_near(mindistsq1, mindistsq2) ) + { + result.push_back(t); + second_sol = true; + } + } + + // we need to test extreme points too + double dsq1 = distanceSq(p, pointAt(from)); + double dsq2 = distanceSq(p, pointAt(to)); + if ( second_sol ) + { + if ( mindistsq2 > dsq1 ) + { + result.clear(); + result.push_back(from); + mindistsq2 = dsq1; + } + else if ( are_near(mindistsq2, dsq) ) + { + result.push_back(from); + } + if ( mindistsq2 > dsq2 ) + { + result.clear(); + result.push_back(to); + } + else if ( are_near(mindistsq2, dsq2) ) + { + result.push_back(to); + } + + } + else + { + if ( result.empty() ) + { + if ( are_near(dsq1, dsq2) ) + { + result.push_back(from); + result.push_back(to); + } + else if ( dsq2 > dsq1 ) + { + result.push_back(from); + } + else + { + result.push_back(to); + } + } + } + + return result; +} +#endif + +/** @brief Convert the passed intersections to curve time parametrization + * and filter out any invalid intersections. + */ +std::vector<ShapeIntersection> EllipticalArc::_filterIntersections(std::vector<ShapeIntersection> &&xs, + bool is_first) const +{ + std::vector<ShapeIntersection> result; + result.reserve(xs.size()); + for (auto &xing : xs) { + if (_validateIntersection(xing, is_first)) { + result.emplace_back(std::move(xing)); + } + } + return result; +} + +/** @brief Convert the passed intersection to curve time and check whether the intersection + * is numerically sane. + * + * @param xing The intersection to convert to curve time and to validate. + * @param is_first If true, this arc is the first of the intersected curves; if false, it's second. + * @return Whether the intersection is valid. + * + * Note that the intersection is guaranteed to be converted only if the return value is true. + */ +bool EllipticalArc::_validateIntersection(ShapeIntersection &xing, bool is_first) const +{ + static auto const UNIT_INTERVAL = Interval(0, 1); + constexpr auto EPS = 1e-4; + + Coord &t = is_first ? xing.first : xing.second; + if (!are_near_rel(_ellipse.pointAt(t), xing.point(), EPS)) { + return false; + } + + t = timeAtAngle(t); + if (!UNIT_INTERVAL.contains(t)) { + return false; + } + if (!are_near_rel(pointAt(t), xing.point(), EPS)) { + return false; + } + return true; +} + +std::vector<CurveIntersection> EllipticalArc::intersect(Curve const &other, Coord eps) const +{ + if (isLineSegment()) { + LineSegment ls(_initial_point, _final_point); + return ls.intersect(other, eps); + } + + if (other.isLineSegment()) { + LineSegment ls(other.initialPoint(), other.finalPoint()); + return _filterIntersections(_ellipse.intersect(ls), true); + } + + if (auto bez = dynamic_cast<BezierCurve const *>(&other)) { + return _filterIntersections(_ellipse.intersect(bez->fragment()), true); + } + + if (auto arc = dynamic_cast<EllipticalArc const *>(&other)) { + std::vector<CurveIntersection> crossings; + try { + crossings = _ellipse.intersect(arc->_ellipse); + } catch (InfinitelyManySolutions &) { + // This could happen if the two arcs come from the same ellipse. + return _intersectSameEllipse(arc); + } + return arc->_filterIntersections(_filterIntersections(std::move(crossings), true), false); + } + + // in case someone wants to make a custom curve type + auto result = other.intersect(*this, eps); + transpose_in_place(result); + return result; +} + +/** @brief Check if two arcs on the same ellipse intersect/overlap. + * + * @param other Another elliptical arc on the same ellipse as this one. + * @return If the arcs overlap, the returned vector contains synthesized intersections + * at the start and end of the overlap. + * If the arcs do not overlap, an empty vector is returned. + */ +std::vector<ShapeIntersection> EllipticalArc::_intersectSameEllipse(EllipticalArc const *other) const +{ + assert(_ellipse == other->_ellipse); + auto const &other_angles = other->angularInterval(); + std::vector<ShapeIntersection> result; + + /// A closure to create an "intersection" at the prescribed angle. + auto const synthesize_intersection = [&](Angle angle) { + auto const time = timeAtAngle(angle); + if (result.end() == std::find_if(result.begin(), result.end(), + [=](ShapeIntersection const &xing) -> bool { + return xing.first == time; + })) + { + result.emplace_back(time, other->timeAtAngle(angle), _ellipse.pointAt(angle)); + } + }; + + for (auto a : {_angles.initialAngle(), _angles.finalAngle()}) { + if (other_angles.contains(a)) { + synthesize_intersection(a); + } + } + for (auto a : {other_angles.initialAngle(), other_angles.finalAngle()}) { + if (_angles.contains(a)) { + synthesize_intersection(a); + } + } + return result; +} + +void EllipticalArc::_updateCenterAndAngles() +{ + // See: http://www.w3.org/TR/SVG/implnote.html#ArcImplementationNotes + Point d = initialPoint() - finalPoint(); + Point mid = middle_point(initialPoint(), finalPoint()); + + auto degenerate_ellipse = [&] { + _ellipse = Ellipse(); + _ellipse.setCenter(initialPoint()); + _angles = AngleInterval(); + _large_arc = false; + }; + + // if ip = fp, the arc contains no other points + if (initialPoint() == finalPoint()) { + degenerate_ellipse(); + return; + } + + // rays should be positive + _ellipse.setRays(std::fabs(ray(X)), std::fabs(ray(Y))); + + if (isChord()) { + _ellipse.setRays(L2(d) / 2, 0); + _ellipse.setRotationAngle(atan2(d)); + _ellipse.setCenter(mid); + _angles.setAngles(0, M_PI); + _angles.setSweep(false); + _large_arc = false; + return; + } + + Rotate rot(rotationAngle()); // the matrix in F.6.5.3 + Rotate invrot = rot.inverse(); // the matrix in F.6.5.1 + + Point r = rays(); + Point p = d / 2 * invrot; // x', y' in F.6.5.1 + Point c(0,0); // cx', cy' in F.6.5.2 + + // Correct out-of-range radii + Coord lambda = hypot(p[X]/r[X], p[Y]/r[Y]); + if (lambda > 1) { + r *= lambda; + _ellipse.setRays(r); + _ellipse.setCenter(mid); + } else { + // evaluate F.6.5.2 + Coord rxry = r[X]*r[X] * r[Y]*r[Y]; + Coord pxry = p[X]*p[X] * r[Y]*r[Y]; + Coord rxpy = r[X]*r[X] * p[Y]*p[Y]; + Coord const denominator = rxpy + pxry; + if (denominator == 0.0) { + degenerate_ellipse(); + return; + } + Coord rad = (rxry - pxry - rxpy) / denominator; + // normally rad should never be negative, but numerical inaccuracy may cause this + if (rad > 0) { + rad = std::sqrt(rad); + if (sweep() == _large_arc) { + rad = -rad; + } + c = rad * Point(r[X]*p[Y]/r[Y], -r[Y]*p[X]/r[X]); + _ellipse.setCenter(c * rot + mid); + } else { + _ellipse.setCenter(mid); + } + } + + // Compute start and end angles. + // If the ellipse was enlarged, c will be zero - this is correct. + Point sp((p[X] - c[X]) / r[X], (p[Y] - c[Y]) / r[Y]); + Point ep((-p[X] - c[X]) / r[X], (-p[Y] - c[Y]) / r[Y]); + Point v(1, 0); + + _angles.setInitial(angle_between(v, sp)); + _angles.setFinal(angle_between(v, ep)); + + /*double sweep_angle = angle_between(sp, ep); + if (!sweep() && sweep_angle > 0) sweep_angle -= 2*M_PI; + if (sweep() && sweep_angle < 0) sweep_angle += 2*M_PI;*/ +} + +D2<SBasis> EllipticalArc::toSBasis() const +{ + if (isChord()) { + return chord().toSBasis(); + } + + D2<SBasis> arc; + // the interval of parametrization has to be [0,1] + Coord et = initialAngle().radians() + sweepAngle(); + Linear param(initialAngle().radians(), et); + Coord cosrot, sinrot; + sincos(rotationAngle(), sinrot, cosrot); + + // order = 4 seems to be enough to get a perfect looking elliptical arc + SBasis arc_x = ray(X) * cos(param,4); + SBasis arc_y = ray(Y) * sin(param,4); + arc[0] = arc_x * cosrot - arc_y * sinrot + Linear(center(X), center(X)); + arc[1] = arc_x * sinrot + arc_y * cosrot + Linear(center(Y), center(Y)); + + // ensure that endpoints remain exact + for ( int d = 0 ; d < 2 ; d++ ) { + arc[d][0][0] = initialPoint()[d]; + arc[d][0][1] = finalPoint()[d]; + } + + return arc; +} + +// All operations that do not contain skew can be evaluated +// without passing through the implicit form of the ellipse, +// which preserves precision. + +void EllipticalArc::operator*=(Translate const &tr) +{ + _initial_point *= tr; + _final_point *= tr; + _ellipse *= tr; +} + +void EllipticalArc::operator*=(Scale const &s) +{ + _initial_point *= s; + _final_point *= s; + _ellipse *= s; +} + +void EllipticalArc::operator*=(Rotate const &r) +{ + _initial_point *= r; + _final_point *= r; + _ellipse *= r; +} + +void EllipticalArc::operator*=(Zoom const &z) +{ + _initial_point *= z; + _final_point *= z; + _ellipse *= z; +} + +void EllipticalArc::operator*=(Affine const& m) +{ + if (isChord()) { + _initial_point *= m; + _final_point *= m; + _ellipse.setCenter(middle_point(_initial_point, _final_point)); + _ellipse.setRays(0, 0); + _ellipse.setRotationAngle(0); + return; + } + + _initial_point *= m; + _final_point *= m; + _ellipse *= m; + if (m.det() < 0) { + _angles.setSweep(!sweep()); + } + + // ellipse transformation does not preserve its functional form, + // i.e. e.pointAt(0.5)*m and (e*m).pointAt(0.5) can be different. + // We need to recompute start / end angles. + _angles.setInitial(_ellipse.timeAt(_initial_point)); + _angles.setFinal(_ellipse.timeAt(_final_point)); +} + +bool EllipticalArc::operator==(Curve const &c) const +{ + EllipticalArc const *other = dynamic_cast<EllipticalArc const *>(&c); + if (!other) return false; + if (_initial_point != other->_initial_point) return false; + if (_final_point != other->_final_point) return false; + // TODO: all arcs with ellipse rays which are too small + // and fall back to a line should probably be equal + if (rays() != other->rays()) return false; + if (rotationAngle() != other->rotationAngle()) return false; + if (_large_arc != other->_large_arc) return false; + if (sweep() != other->sweep()) return false; + return true; +} + +bool EllipticalArc::isNear(Curve const &c, Coord precision) const +{ + EllipticalArc const *other = dynamic_cast<EllipticalArc const *>(&c); + if (!other) { + if (isChord()) { + return c.isNear(chord(), precision); + } + return false; + } + + if (!are_near(_initial_point, other->_initial_point, precision)) return false; + if (!are_near(_final_point, other->_final_point, precision)) return false; + if (isChord() && other->isChord()) return true; + + if (sweep() != other->sweep()) return false; + if (!are_near(_ellipse, other->_ellipse, precision)) return false; + return true; +} + +void EllipticalArc::feed(PathSink &sink, bool moveto_initial) const +{ + if (moveto_initial) { + sink.moveTo(_initial_point); + } + sink.arcTo(ray(X), ray(Y), rotationAngle(), _large_arc, sweep(), _final_point); +} + +int EllipticalArc::winding(Point const &p) const +{ + using std::swap; + + double sinrot, cosrot; + sincos(rotationAngle(), sinrot, cosrot); + + Angle ymin_a = std::atan2( ray(Y) * cosrot, ray(X) * sinrot ); + Angle ymax_a = ymin_a + M_PI; + + Point ymin = pointAtAngle(ymin_a); + Point ymax = pointAtAngle(ymax_a); + if (ymin[Y] > ymax[Y]) { + swap(ymin, ymax); + swap(ymin_a, ymax_a); + } + + if (!Interval(ymin[Y], ymax[Y]).lowerContains(p[Y])) { + return 0; + } + + bool const left = cross(ymax - ymin, p - ymin) > 0; + bool const inside = _ellipse.contains(p); + if (_angles.isFull()) { + if (inside) { + return sweep() ? 1 : -1; + } + return 0; + } + bool const includes_ymin = _angles.contains(ymin_a); + bool const includes_ymax = _angles.contains(ymax_a); + + AngleInterval rarc(ymin_a, ymax_a, true), + larc(ymax_a, ymin_a, true); + + // we'll compute the result for an arc in the direction of increasing angles + // and then negate if necessary + Angle ia = initialAngle(), fa = finalAngle(); + Point ip = _initial_point, fp = _final_point; + if (!sweep()) { + swap(ia, fa); + swap(ip, fp); + } + + bool const initial_left = larc.contains(ia); + bool const final_left = larc.contains(fa); + + bool intersects_left = false, intersects_right = false; + if (inside || left) { + // The point is inside the ellipse or to the left of it, so the rightwards horizontal ray + // may intersect the part of the arc contained in the right half of the ellipse. + // There are four ways in which this can happen. + + intersects_right = + // Possiblity 1: the arc extends into the right half through the min-Y point + // and the ray intersects this extension: + (includes_ymin && !final_left && Interval(ymin[Y], fp[Y]).lowerContains(p[Y])) + || + // Possibility 2: the arc starts and ends within the right half (hence, it cannot be the + // "large arc") and the ray's Y-coordinate is within the Y-coordinate range of the arc: + (!initial_left && !final_left && !largeArc() && Interval(ip[Y], fp[Y]).lowerContains(p[Y])) + || + // Possibility 3: the arc starts in the right half and continues through the max-Y + // point into the left half: + (!initial_left && includes_ymax && Interval(ip[Y], ymax[Y]).lowerContains(p[Y])) + || + // Possibility 4: the entire right half of the ellipse is contained in the arc. + (initial_left && final_left && includes_ymin && includes_ymax); + } + if (left && !inside) { + // The point is to the left of the ellipse, so the rightwards horizontal ray + // may intersect the part of the arc contained in the left half of the ellipse. + // There are four ways in which this can happen. + + intersects_left = + // Possibility 1: the arc starts in the left half and continues through the min-Y + // point into the right half: + (includes_ymin && initial_left && Interval(ymin[Y], ip[Y]).lowerContains(p[Y])) + || + // Possibility 2: the arc starts and ends within the left half (hence, it cannot be the + // "large arc") and the ray's Y-coordinate is within the Y-coordinate range of the arc: + (initial_left && final_left && !largeArc() && Interval(ip[Y], fp[Y]).lowerContains(p[Y])) + || + // Possibility 3: the arc extends into the left half through the max-Y point + // and the ray intersects this extension: + (final_left && includes_ymax && Interval(fp[Y], ymax[Y]).lowerContains(p[Y])) + || + // Possibility 4: the entire left half of the ellipse is contained in the arc. + (!initial_left && !final_left && includes_ymin && includes_ymax); + + } + int const winding_assuming_increasing_angles = (int)intersects_right - (int)intersects_left; + return sweep() ? winding_assuming_increasing_angles : -winding_assuming_increasing_angles; +} + +std::ostream &operator<<(std::ostream &out, EllipticalArc const &ea) +{ + out << "EllipticalArc(" + << ea.initialPoint() << ", " + << format_coord_nice(ea.ray(X)) << ", " << format_coord_nice(ea.ray(Y)) << ", " + << format_coord_nice(ea.rotationAngle()) << ", " + << "large_arc=" << (ea.largeArc() ? "true" : "false") << ", " + << "sweep=" << (ea.sweep() ? "true" : "false") << ", " + << ea.finalPoint() << ")"; + return out; +} + +} // 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:fileencoding=utf-8:textwidth=99 : + |