From c883d7627a479c8c5b6a9f77b9841fa5631572ad Mon Sep 17 00:00:00 2001 From: Krzysztof Kosi??ski Date: Mon, 27 Apr 2015 19:39:29 -0400 Subject: 2Geom sync - initial commit (bzr r14059.2.1) --- src/2geom/ellipse.cpp | 231 ++++++++++++++++++++++++-------------------------- 1 file changed, 112 insertions(+), 119 deletions(-) (limited to 'src/2geom/ellipse.cpp') diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index 2686844b2..4e26707ef 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -1,10 +1,11 @@ -/* - * Ellipse Curve - * +/** @file + * @brief Ellipse shape + *//* * Authors: - * Marco Cecchetti + * Marco Cecchetti + * Krzysztof Kosiński * - * Copyright 2008 authors + * Copyright 2008-2014 Authors * * This library is free software; you can redistribute it and/or * modify it either under the terms of the GNU Lesser General Public @@ -30,33 +31,28 @@ * the specific language governing rights and limitations. */ - #include <2geom/ellipse.h> #include <2geom/svg-elliptical-arc.h> #include <2geom/numeric/fitting-tool.h> #include <2geom/numeric/fitting-model.h> -using std::swap; +namespace Geom { -namespace Geom -{ - -void Ellipse::set(double A, double B, double C, double D, double E, double F) +void Ellipse::setCoefficients(double A, double B, double C, double D, double E, double F) { double den = 4*A*C - B*B; - if ( den == 0 ) - { - THROW_LOGICALERROR("den == 0, while computing ellipse centre"); + if (den == 0) { + THROW_RANGEERROR("den == 0, while computing ellipse centre"); } - m_centre[X] = (B*E - 2*C*D) / den; - m_centre[Y] = (B*D - 2*A*E) / den; + _center[X] = (B*E - 2*C*D) / den; + _center[Y] = (B*D - 2*A*E) / den; // evaluate the a coefficient of the ellipse equation in normal form // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1 // where b = a*B , c = a*C, (cx,cy) == centre - double num = A * sqr(m_centre[X]) - + B * m_centre[X] * m_centre[Y] - + C * sqr(m_centre[Y]) + double num = A * sqr(_center[X]) + + B * _center[X] * _center[Y] + + C * sqr(_center[Y]) - F; @@ -64,73 +60,74 @@ void Ellipse::set(double A, double B, double C, double D, double E, double F) double rot = std::atan2( -B, -(A - C) )/2; // std::cerr << "rot = " << rot << std::endl; bool swap_axes = false; - if ( are_near(rot, 0) ) rot = 0; - if ( are_near(rot, M_PI/2) || rot < 0 ) - { + + if (rot >= M_PI/2 || rot < 0) { swap_axes = true; } // evaluate the length of the ellipse rays - double cosrot = std::cos(rot); - double sinrot = std::sin(rot); + double sinrot, cosrot; + sincos(rot, sinrot, cosrot); double cos2 = cosrot * cosrot; double sin2 = sinrot * sinrot; double cossin = cosrot * sinrot; den = A * cos2 + B * cossin + C * sin2; - if ( den == 0 ) - { - THROW_LOGICALERROR("den == 0, while computing 'rx' coefficient"); + if (den == 0) { + THROW_RANGEERROR("den == 0, while computing 'rx' coefficient"); } - double rx2 = num/den; - if ( rx2 < 0 ) - { - THROW_LOGICALERROR("rx2 < 0, while computing 'rx' coefficient"); + double rx2 = num / den; + if (rx2 < 0) { + THROW_RANGEERROR("rx2 < 0, while computing 'rx' coefficient"); } double rx = std::sqrt(rx2); den = C * cos2 - B * cossin + A * sin2; - if ( den == 0 ) - { - THROW_LOGICALERROR("den == 0, while computing 'ry' coefficient"); + if (den == 0) { + THROW_RANGEERROR("den == 0, while computing 'ry' coefficient"); } - double ry2 = num/den; - if ( ry2 < 0 ) - { - THROW_LOGICALERROR("ry2 < 0, while computing 'rx' coefficient"); + double ry2 = num / den; + if (ry2 < 0) { + THROW_RANGEERROR("ry2 < 0, while computing 'rx' coefficient"); } double ry = std::sqrt(ry2); // the solution is not unique so we choose always the ellipse // with a rotation angle between 0 and PI/2 - if ( swap_axes ) swap(rx, ry); - if ( are_near(rot, M_PI/2) - || are_near(rot, -M_PI/2) - || are_near(rx, ry) ) - { + if (swap_axes) { + std::swap(rx, ry); + } + + if (rx == ry) { rot = 0; } - else if ( rot < 0 ) - { + if (rot < 0) { rot += M_PI/2; } - m_ray[X] = rx; - m_ray[Y] = ry; - m_angle = rot; + _rays[X] = rx; + _rays[Y] = ry; + _angle = rot; +} + + +Affine Ellipse::unitCircleTransform() const +{ + Affine ret = Scale(ray(X), ray(Y)) * Rotate(_angle); + ret.setTranslation(center()); + return ret; } -std::vector Ellipse::implicit_form_coefficients() const +std::vector Ellipse::coefficients() const { - if (ray(X) == 0 || ray(Y) == 0) - { - THROW_LOGICALERROR("a degenerate ellipse doesn't own an implicit form"); + if (ray(X) == 0 || ray(Y) == 0) { + THROW_RANGEERROR("a degenerate ellipse doesn't have an implicit form"); } std::vector coeff(6); - double cosrot = std::cos(rot_angle()); - double sinrot = std::sin(rot_angle()); + double cosrot, sinrot; + sincos(_angle, sinrot, cosrot); double cos2 = cosrot * cosrot; double sin2 = sinrot * sinrot; double cossin = cosrot * sinrot; @@ -150,18 +147,16 @@ std::vector Ellipse::implicit_form_coefficients() const } -void Ellipse::set(std::vector const& points) +void Ellipse::fit(std::vector const &points) { size_t sz = points.size(); - if (sz < 5) - { + if (sz < 5) { THROW_RANGEERROR("fitting error: too few points passed"); } NL::LFMEllipse model; NL::least_squeares_fitter fitter(model, sz); - for (size_t i = 0; i < sz; ++i) - { + for (size_t i = 0; i < sz; ++i) { fitter.append(points[i]); } fitter.update(); @@ -172,12 +167,12 @@ void Ellipse::set(std::vector const& points) EllipticalArc * -Ellipse::arc(Point const& initial, Point const& inner, Point const& final, +Ellipse::arc(Point const &ip, Point const &inner, Point const &fp, bool _svg_compliant) { - Point sp_cp = initial - center(); - Point ep_cp = final - center(); - Point ip_cp = inner - center(); + Point sp_cp = ip - center(); + Point ep_cp = fp - center(); + Point ip_cp = inner - center(); double angle1 = angle_between(sp_cp, ep_cp); double angle2 = angle_between(sp_cp, ip_cp); @@ -186,28 +181,19 @@ Ellipse::arc(Point const& initial, Point const& inner, Point const& final, bool large_arc_flag = true; bool sweep_flag = true; - if ( angle1 > 0 ) - { - if ( angle2 > 0 && angle3 > 0 ) - { + if (angle1 > 0) { + if (angle2 > 0 && angle3 > 0) { large_arc_flag = false; sweep_flag = true; - } - else - { + } else { large_arc_flag = true; sweep_flag = false; } - } - else - { - if ( angle2 < 0 && angle3 < 0 ) - { + } else { + if (angle2 < 0 && angle3 < 0) { large_arc_flag = false; sweep_flag = false; - } - else - { + } else { large_arc_flag = true; sweep_flag = true; } @@ -215,66 +201,73 @@ Ellipse::arc(Point const& initial, Point const& inner, Point const& final, EllipticalArc *ret_arc; if (_svg_compliant) { - ret_arc = new SVGEllipticalArc(initial, ray(X), ray(Y), rot_angle(), - large_arc_flag, sweep_flag, final); + ret_arc = new SVGEllipticalArc(ip, ray(X), ray(Y), rotationAngle(), + large_arc_flag, sweep_flag, fp); } else { - ret_arc = new EllipticalArc(initial, ray(X), ray(Y), rot_angle(), - large_arc_flag, sweep_flag, final); + ret_arc = new EllipticalArc(ip, ray(X), ray(Y), rotationAngle(), + large_arc_flag, sweep_flag, fp); } return ret_arc; } -Ellipse Ellipse::transformed(Affine const& m) const +Ellipse &Ellipse::operator*=(Rotate const &r) { - double cosrot = std::cos(rot_angle()); - double sinrot = std::sin(rot_angle()); - Affine A( ray(X) * cosrot, ray(X) * sinrot, - -ray(Y) * sinrot, ray(Y) * cosrot, - 0, 0 ); - Point new_center = center() * m; - Affine M = m.withoutTranslation(); - Affine AM = A * M; - if ( are_near(std::sqrt(fabs(AM.det())), 0) ) - { + _angle += r.angle(); + // keep the angle in the first quadrant + if (_angle < 0) { + _angle += M_PI; + } + if (_angle >= M_PI/2) { + std::swap(_rays[X], _rays[Y]); + _angle -= M_PI/2; + } + _center *= r; + return *this; +} + +Ellipse &Ellipse::operator*=(Affine const& m) +{ + Rotate a(_angle); + Affine mwot = m.withoutTranslation(); + Affine am = a * mwot; + if (are_near(am.descrim(), 0)) { double angle; - if (AM[0] != 0) - { - angle = std::atan2(AM[2], AM[0]); - } - else if (AM[1] != 0) - { - angle = std::atan2(AM[3], AM[1]); - } - else - { + if (am[0] != 0) { + angle = std::atan2(am[2], am[0]); + } else if (am[1] != 0) { + angle = std::atan2(am[3], am[1]); + } else { angle = M_PI/2; } - Point V(std::cos(angle), std::sin(angle)); - V *= AM; - double rx = L2(V); - angle = atan2(V); - return Ellipse(new_center[X], new_center[Y], rx, 0, angle); + Point v; + sincos(angle, v[X], v[Y]); + v *= am; + _angle = atan2(v); + _center *= m; + _rays[X] = L2(v); + _rays[Y] = 0; + return *this; } - std::vector coeff = implicit_form_coefficients(); - Affine Q( coeff[0], coeff[1]/2, + std::vector coeff = coefficients(); + Affine q( coeff[0], coeff[1]/2, coeff[1]/2, coeff[2], 0, 0 ); - Affine invm = M.inverse(); - Q = invm * Q ; - swap( invm[1], invm[2] ); - Q *= invm; - Ellipse e(Q[0], 2*Q[1], Q[3], 0, 0, -1); - e.m_centre = new_center; + Affine invm = mwot.inverse(); + q = invm * q ; + std::swap(invm[1], invm[2]); + q *= invm; + setCoefficients(q[0], 2*q[1], q[3], 0, 0, -1); + _center *= m; - return e; + return *this; } Ellipse::Ellipse(Geom::Circle const &c) { - m_centre = c.center(); - m_ray[X] = m_ray[Y] = c.ray(); + _center = c.center(); + _rays[X] = _rays[Y] = c.radius(); } } // end namespace Geom -- cgit v1.2.3 From 6a9762c7603a32c7ec5cc0aaed8048d84daee6e8 Mon Sep 17 00:00:00 2001 From: Krzysztof Kosi??ski Date: Thu, 30 Apr 2015 11:17:07 +0200 Subject: Update 2Geom to r2347 (bzr r14059.2.3) --- src/2geom/ellipse.cpp | 50 ++++++++++++++++++++++++-------------------------- 1 file changed, 24 insertions(+), 26 deletions(-) (limited to 'src/2geom/ellipse.cpp') diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index 4e26707ef..9e927c9e8 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -170,33 +170,31 @@ EllipticalArc * Ellipse::arc(Point const &ip, Point const &inner, Point const &fp, bool _svg_compliant) { - Point sp_cp = ip - center(); - Point ep_cp = fp - center(); - Point ip_cp = inner - center(); - - double angle1 = angle_between(sp_cp, ep_cp); - double angle2 = angle_between(sp_cp, ip_cp); - double angle3 = angle_between(ip_cp, ep_cp); - - bool large_arc_flag = true; - bool sweep_flag = true; + // This is resistant to degenerate ellipses: + // both flags evaluate to false in that case. + + bool large_arc_flag = false; + bool sweep_flag = false; + + // Determination of large arc flag: + // The arc is larger than half of the ellipse if the inner point + // is on the same side of the line going from the initial + // to the final point as the center of the ellipse + Line chord(ip, fp); + Point versor = fp - ip; + double sdist_c = cross(versor, _center - ip); + double sdist_inner = cross(versor, inner - ip); + + // if we have exactly half of an arc, do not set the large flag. + if (sdist_c != 0 && sgn(sdist_c) == sgn(sdist_inner)) { + large_arc_flag = true; + } - if (angle1 > 0) { - if (angle2 > 0 && angle3 > 0) { - large_arc_flag = false; - sweep_flag = true; - } else { - large_arc_flag = true; - sweep_flag = false; - } - } else { - if (angle2 < 0 && angle3 < 0) { - large_arc_flag = false; - sweep_flag = false; - } else { - large_arc_flag = true; - sweep_flag = true; - } + // Determination of sweep flag: + // If the inner point is on the left side of the ip-fp line, + // we go in clockwise direction. + if (sdist_inner < 0) { + sweep_flag = true; } EllipticalArc *ret_arc; -- cgit v1.2.3 From 7b7620b5dc7f3c57be46a5078f0bec2e3868f553 Mon Sep 17 00:00:00 2001 From: Krzysztof Kosi??ski Date: Fri, 8 May 2015 14:09:07 +0200 Subject: Fix ellipse problems (bzr r14059.2.5) --- src/2geom/ellipse.cpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) (limited to 'src/2geom/ellipse.cpp') diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index 9e927c9e8..f23c9b24e 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -180,7 +180,6 @@ Ellipse::arc(Point const &ip, Point const &inner, Point const &fp, // The arc is larger than half of the ellipse if the inner point // is on the same side of the line going from the initial // to the final point as the center of the ellipse - Line chord(ip, fp); Point versor = fp - ip; double sdist_c = cross(versor, _center - ip); double sdist_inner = cross(versor, inner - ip); @@ -225,9 +224,11 @@ Ellipse &Ellipse::operator*=(Rotate const &r) Ellipse &Ellipse::operator*=(Affine const& m) { - Rotate a(_angle); + Affine a = Scale(ray(X), ray(Y)) * Rotate(_angle); Affine mwot = m.withoutTranslation(); Affine am = a * mwot; + Point new_center = _center * m; + if (are_near(am.descrim(), 0)) { double angle; if (am[0] != 0) { @@ -237,13 +238,11 @@ Ellipse &Ellipse::operator*=(Affine const& m) } else { angle = M_PI/2; } - Point v; - sincos(angle, v[X], v[Y]); - v *= am; - _angle = atan2(v); - _center *= m; + Point v = Point::polar(angle) * am; + _center = new_center; _rays[X] = L2(v); _rays[Y] = 0; + _angle = atan2(v); return *this; } @@ -257,7 +256,7 @@ Ellipse &Ellipse::operator*=(Affine const& m) std::swap(invm[1], invm[2]); q *= invm; setCoefficients(q[0], 2*q[1], q[3], 0, 0, -1); - _center *= m; + _center = new_center; return *this; } -- cgit v1.2.3 From 25fa09178b7d0d0befa708e93ea5316ef381caa0 Mon Sep 17 00:00:00 2001 From: Krzysztof Kosi??ski Date: Fri, 22 May 2015 10:23:27 +0200 Subject: Update to 2Geom revision 2396 (bzr r14059.2.16) --- src/2geom/ellipse.cpp | 421 +++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 364 insertions(+), 57 deletions(-) (limited to 'src/2geom/ellipse.cpp') diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index f23c9b24e..46c60d85d 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -32,12 +32,18 @@ */ #include <2geom/ellipse.h> -#include <2geom/svg-elliptical-arc.h> +#include <2geom/elliptical-arc.h> #include <2geom/numeric/fitting-tool.h> #include <2geom/numeric/fitting-model.h> namespace Geom { +Ellipse::Ellipse(Geom::Circle const &c) + : _center(c.center()) + , _rays(c.radius(), c.radius()) + , _angle(0) +{} + void Ellipse::setCoefficients(double A, double B, double C, double D, double E, double F) { double den = 4*A*C - B*B; @@ -57,17 +63,11 @@ void Ellipse::setCoefficients(double A, double B, double C, double D, double E, //evaluate ellipse rotation angle - double rot = std::atan2( -B, -(A - C) )/2; -// std::cerr << "rot = " << rot << std::endl; - bool swap_axes = false; - - if (rot >= M_PI/2 || rot < 0) { - swap_axes = true; - } + _angle = std::atan2( -B, -(A - C) )/2; // evaluate the length of the ellipse rays double sinrot, cosrot; - sincos(rot, sinrot, cosrot); + sincos(_angle, sinrot, cosrot); double cos2 = cosrot * cosrot; double sin2 = sinrot * sinrot; double cossin = cosrot * sinrot; @@ -80,7 +80,7 @@ void Ellipse::setCoefficients(double A, double B, double C, double D, double E, if (rx2 < 0) { THROW_RANGEERROR("rx2 < 0, while computing 'rx' coefficient"); } - double rx = std::sqrt(rx2); + _rays[X] = std::sqrt(rx2); den = C * cos2 - B * cossin + A * sin2; if (den == 0) { @@ -90,24 +90,11 @@ void Ellipse::setCoefficients(double A, double B, double C, double D, double E, if (ry2 < 0) { THROW_RANGEERROR("ry2 < 0, while computing 'rx' coefficient"); } - double ry = std::sqrt(ry2); + _rays[Y] = std::sqrt(ry2); // the solution is not unique so we choose always the ellipse // with a rotation angle between 0 and PI/2 - if (swap_axes) { - std::swap(rx, ry); - } - - if (rx == ry) { - rot = 0; - } - if (rot < 0) { - rot += M_PI/2; - } - - _rays[X] = rx; - _rays[Y] = ry; - _angle = rot; + makeCanonical(); } @@ -118,14 +105,49 @@ Affine Ellipse::unitCircleTransform() const return ret; } +Affine Ellipse::inverseUnitCircleTransform() const +{ + if (ray(X) == 0 || ray(Y) == 0) { + THROW_RANGEERROR("a degenerate ellipse doesn't have an inverse unit circle transform"); + } + Affine ret = Translate(-center()) * Rotate(-_angle) * Scale(1/ray(X), 1/ray(Y)); + return ret; +} + + +LineSegment Ellipse::axis(Dim2 d) const +{ + Point a(0, 0), b(0, 0); + a[d] = -1; + b[d] = 1; + LineSegment ls(a, b); + ls.transform(unitCircleTransform()); + return ls; +} + +LineSegment Ellipse::semiaxis(Dim2 d, int sign) const +{ + Point a(0, 0), b(0, 0); + b[d] = sgn(sign); + LineSegment ls(a, b); + ls.transform(unitCircleTransform()); + return ls; +} + std::vector Ellipse::coefficients() const +{ + std::vector c(6); + coefficients(c[0], c[1], c[2], c[3], c[4], c[5]); + return c; +} + +void Ellipse::coefficients(Coord &A, Coord &B, Coord &C, Coord &D, Coord &E, Coord &F) const { if (ray(X) == 0 || ray(Y) == 0) { THROW_RANGEERROR("a degenerate ellipse doesn't have an implicit form"); } - std::vector coeff(6); double cosrot, sinrot; sincos(_angle, sinrot, cosrot); double cos2 = cosrot * cosrot; @@ -134,16 +156,15 @@ std::vector Ellipse::coefficients() const double invrx2 = 1 / (ray(X) * ray(X)); double invry2 = 1 / (ray(Y) * ray(Y)); - coeff[0] = invrx2 * cos2 + invry2 * sin2; - coeff[1] = 2 * (invrx2 - invry2) * cossin; - coeff[2] = invrx2 * sin2 + invry2 * cos2; - coeff[3] = -(2 * coeff[0] * center(X) + coeff[1] * center(Y)); - coeff[4] = -(2 * coeff[2] * center(Y) + coeff[1] * center(X)); - coeff[5] = coeff[0] * center(X) * center(X) - + coeff[1] * center(X) * center(Y) - + coeff[2] * center(Y) * center(Y) - - 1; - return coeff; + A = invrx2 * cos2 + invry2 * sin2; + B = 2 * (invrx2 - invry2) * cossin; + C = invrx2 * sin2 + invry2 * cos2; + D = -2 * A * center(X) - B * center(Y); + E = -2 * C * center(Y) - B * center(X); + F = A * center(X) * center(X) + + B * center(X) * center(Y) + + C * center(Y) * center(Y) + - 1; } @@ -167,8 +188,7 @@ void Ellipse::fit(std::vector const &points) EllipticalArc * -Ellipse::arc(Point const &ip, Point const &inner, Point const &fp, - bool _svg_compliant) +Ellipse::arc(Point const &ip, Point const &inner, Point const &fp) { // This is resistant to degenerate ellipses: // both flags evaluate to false in that case. @@ -196,28 +216,14 @@ Ellipse::arc(Point const &ip, Point const &inner, Point const &fp, sweep_flag = true; } - EllipticalArc *ret_arc; - if (_svg_compliant) { - ret_arc = new SVGEllipticalArc(ip, ray(X), ray(Y), rotationAngle(), - large_arc_flag, sweep_flag, fp); - } else { - ret_arc = new EllipticalArc(ip, ray(X), ray(Y), rotationAngle(), - large_arc_flag, sweep_flag, fp); - } + EllipticalArc *ret_arc = new EllipticalArc(ip, ray(X), ray(Y), rotationAngle(), + large_arc_flag, sweep_flag, fp); return ret_arc; } Ellipse &Ellipse::operator*=(Rotate const &r) { _angle += r.angle(); - // keep the angle in the first quadrant - if (_angle < 0) { - _angle += M_PI; - } - if (_angle >= M_PI/2) { - std::swap(_rays[X], _rays[Y]); - _angle -= M_PI/2; - } _center *= r; return *this; } @@ -261,10 +267,311 @@ Ellipse &Ellipse::operator*=(Affine const& m) return *this; } -Ellipse::Ellipse(Geom::Circle const &c) +Ellipse Ellipse::canonicalForm() const +{ + Ellipse result(*this); + result.makeCanonical(); + return result; +} + +void Ellipse::makeCanonical() +{ + if (_rays[X] == _rays[Y]) { + _angle = 0; + return; + } + + if (_angle < 0) { + _angle += M_PI; + } + if (_angle >= M_PI/2) { + std::swap(_rays[X], _rays[Y]); + _angle -= M_PI/2; + } +} + +Point Ellipse::pointAt(Coord t) const +{ + Point p = Point::polar(t); + p *= unitCircleTransform(); + return p; +} + +Coord Ellipse::valueAt(Coord t, Dim2 d) const +{ + // TODO: more efficient version. + return pointAt(t)[d]; +} + +Coord Ellipse::timeAt(Point const &p) const +{ + // degenerate ellipse is basically a reparametrized line segment + if (ray(X) == 0 || ray(Y) == 0) { + if (ray(X) != 0) { + return asin(Line(axis(X)).timeAt(p)); + } else if (ray(Y) != 0) { + return acos(Line(axis(Y)).timeAt(p)); + } else { + return 0; + } + } + Affine iuct = inverseUnitCircleTransform(); + return Angle(atan2(p * iuct)).radians0(); // return a value in [0, 2pi) +} + +std::vector Ellipse::intersect(Line const &line) const +{ + + std::vector result; + + if (line.isDegenerate()) return result; + if (ray(X) == 0 || ray(Y) == 0) { + // TODO intersect with line segment. + return result; + } + + // Ax^2 + Bxy + Cy^2 + Dx + Ey + F + Coord A, B, C, D, E, F; + coefficients(A, B, C, D, E, F); + Affine iuct = inverseUnitCircleTransform(); + + if (line.isHorizontal()) { + // substitute y into the ellipse equation and solve + Coord y = line.initialPoint()[Y]; + std::vector xs = solve_quadratic(A, B*y + D, C*y*y + E*y + F); + for (unsigned i = 0; i < xs.size(); ++i) { + Point p(xs[i], y); + result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p)); + } + return result; + } + if (line.isVertical()) { + // substitute y into the ellipse equation and solve + Coord x = line.initialPoint()[X]; + std::vector ys = solve_quadratic(C, B*x + E, A*x*x + D*x + F); + for (unsigned i = 0; i < ys.size(); ++i) { + Point p(x, ys[i]); + result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p)); + } + return result; + } + + // generic case + Coord a, b, c; + line.coefficients(a, b, c); + + // y = -a/b x - C/B + // TODO: when is it better to substitute X? + Coord q = -a/b; + Coord r = -c/b; + + // substitute that into the ellipse equation, making it quadratic in x + Coord I = A + B*q + C*q*q; // x^2 terms + Coord J = B*r + C*2*q*r + D + E*q; // x^1 terms + Coord K = C*r*r + E*r + F; // x^0 terms + std::vector xs = solve_quadratic(I, J, K); + + for (unsigned i = 0; i < xs.size(); ++i) { + Point p(xs[i], q*xs[i] + r); + result.push_back(ShapeIntersection(atan2(p * iuct), line.timeAt(p), p)); + } + return result; +} + +std::vector Ellipse::intersect(LineSegment const &seg) const +{ + // we simply re-use the procedure for lines and filter out + // results where the line time value is outside of the unit interval. + std::vector result = intersect(Line(seg)); + filter_line_segment_intersections(result); + return result; +} + +std::vector Ellipse::intersect(Ellipse const &other) const +{ + // handle degenerate cases first + if (ray(X) == 0 || ray(Y) == 0) { + + } + // intersection of two ellipses can be solved analytically. + // http://maptools.home.comcast.net/~maptools/BivariateQuadratics.pdf + + Coord A, B, C, D, E, F; + Coord a, b, c, d, e, f; + + // NOTE: the order of coefficients is different to match the convention in the PDF above + // Ax^2 + Bx^2 + Cx + Dy + Exy + F + this->coefficients(A, E, B, C, D, F); + other.coefficients(a, e, b, c, d, f); + + // Assume that Q is the ellipse equation given by uppercase letters + // and R is the equation given by lowercase ones. An intersection exists when + // there is a coefficient mu such that + // mu Q + R = 0 + // + // This can be written in the following way: + // + // | ff cc/2 dd/2 | |1| + // mu Q + R = [1 x y] | cc/2 aa ee/2 | |x| = 0 + // | dd/2 ee/2 bb | |y| + // + // where aa = mu A + a and so on. The determinant can be explicitly written out, + // giving an equation which is cubic in mu and can be solved analytically. + + Coord I, J, K, L; + I = (-E*E*F + 4*A*B*F + C*D*E - A*D*D - B*C*C) / 4; + J = -((E*E - 4*A*B) * f + (2*E*F - C*D) * e + (2*A*D - C*E) * d + + (2*B*C - D*E) * c + (C*C - 4*A*F) * b + (D*D - 4*B*F) * a) / 4; + K = -((e*e - 4*a*b) * F + (2*e*f - c*d) * E + (2*a*d - c*e) * D + + (2*b*c - d*e) * C + (c*c - 4*a*f) * B + (d*d - 4*b*f) * A) / 4; + L = (-e*e*f + 4*a*b*f + c*d*e - a*d*d - b*c*c) / 4; + + std::vector mus = solve_cubic(I, J, K, L); + Coord mu = infinity(); + std::vector result; + + // Now that we have solved for mu, we need to check whether the conic + // determined by mu Q + R is reducible to a product of two lines. If it's not, + // it means that there are no intersections. If it is, the intersections of these + // lines with the original ellipses (if there are any) give the coordinates + // of intersections. + + // Prefer middle root if there are three. + // Out of three possible pairs of lines that go through four points of intersection + // of two ellipses, this corresponds to cross-lines. These intersect the ellipses + // at less shallow angles than the other two options. + if (mus.size() == 3) { + std::swap(mus[1], mus[0]); + } + for (unsigned i = 0; i < mus.size(); ++i) { + Coord aa = mus[i] * A + a; + Coord bb = mus[i] * B + b; + Coord ee = mus[i] * E + e; + Coord delta = ee*ee - 4*aa*bb; + if (delta < 0) continue; + mu = mus[i]; + break; + } + + // if no suitable mu was found, there are no intersections + if (mu == infinity()) return result; + + Coord aa = mu * A + a; + Coord bb = mu * B + b; + Coord cc = mu * C + c; + Coord dd = mu * D + d; + Coord ee = mu * E + e; + Coord ff = mu * F + f; + + Line lines[2]; + + if (aa != 0) { + bb /= aa; cc /= aa; dd /= aa; ee /= aa; ff /= aa; + Coord s = (ee + std::sqrt(ee*ee - 4*bb)) / 2; + Coord q = ee - s; + Coord alpha = (dd - cc*q) / (s - q); + Coord beta = cc - alpha; + + lines[0] = Line(1, q, alpha); + lines[1] = Line(1, s, beta); + } else if (bb != 0) { + cc /= bb; dd /= bb; ee /= bb; ff /= bb; + Coord s = ee; + Coord q = 0; + Coord alpha = cc / ee; + Coord beta = ff * ee / cc; + + lines[0] = Line(q, 1, alpha); + lines[1] = Line(s, 1, beta); + } else { + lines[0] = Line(ee, 0, dd); + lines[1] = Line(0, 1, cc/ee); + } + + // intersect with the obtained lines and report intersections + for (unsigned li = 0; li < 2; ++li) { + std::vector as = intersect(lines[li]); + std::vector bs = other.intersect(lines[li]); + + if (!as.empty() && as.size() == bs.size()) { + for (unsigned i = 0; i < as.size(); ++i) { + ShapeIntersection ix(as[i].first, bs[i].first, + middle_point(as[i].point(), bs[i].point())); + result.push_back(ix); + } + } + } + return result; +} + +std::vector Ellipse::intersect(D2 const &b) const +{ + Coord A, B, C, D, E, F; + coefficients(A, B, C, D, E, F); + + Bezier x = A*b[X]*b[X] + B*b[X]*b[Y] + C*b[Y]*b[Y] + D*b[X] + E*b[Y] + F; + std::vector r = x.roots(); + + std::vector result; + for (unsigned i = 0; i < r.size(); ++i) { + Point p = b.valueAt(r[i]); + result.push_back(ShapeIntersection(timeAt(p), r[i], p)); + } + return result; +} + +bool Ellipse::operator==(Ellipse const &other) const +{ + if (_center != other._center) return false; + + Ellipse a = this->canonicalForm(); + Ellipse b = other.canonicalForm(); + + if (a._rays != b._rays) return false; + if (a._angle != b._angle) return false; + + return true; +} + + +bool are_near(Ellipse const &a, Ellipse const &b, Coord precision) +{ + // We want to know whether no point on ellipse a is further than precision + // from the corresponding point on ellipse b. To check this, we compute + // the four extreme points at the end of each ray for each ellipse + // and check whether they are sufficiently close. + + // First, we need to correct the angles on the ellipses, so that they are + // no further than M_PI/4 apart. This can always be done by rotating + // and exchanging axes. + Ellipse ac = a, bc = b; + if (distance(ac.rotationAngle(), bc.rotationAngle()).radians0() >= M_PI/2) { + ac.setRotationAngle(ac.rotationAngle() + M_PI); + } + if (distance(ac.rotationAngle(), bc.rotationAngle()) >= M_PI/4) { + Angle d1 = distance(ac.rotationAngle() + M_PI/2, bc.rotationAngle()); + Angle d2 = distance(ac.rotationAngle() - M_PI/2, bc.rotationAngle()); + Coord adj = d1.radians0() < d2.radians0() ? M_PI/2 : -M_PI/2; + ac.setRotationAngle(ac.rotationAngle() + adj); + ac.setRays(ac.ray(Y), ac.ray(X)); + } + + // Do the actual comparison by computing four points on each ellipse. + Point tps[] = {Point(1,0), Point(0,1), Point(-1,0), Point(0,-1)}; + for (unsigned i = 0; i < 4; ++i) { + if (!are_near(tps[i] * ac.unitCircleTransform(), + tps[i] * bc.unitCircleTransform(), + precision)) + return false; + } + return true; +} + +std::ostream &operator<<(std::ostream &out, Ellipse const &e) { - _center = c.center(); - _rays[X] = _rays[Y] = c.radius(); + out << "Ellipse(" << e.center() << ", " << e.rays() + << ", " << format_coord_nice(e.rotationAngle()) << ")"; + return out; } } // end namespace Geom -- cgit v1.2.3 From 60437ac397d41678daba5daece227240e8ddd364 Mon Sep 17 00:00:00 2001 From: Krzysztof Kosi??ski Date: Sat, 4 Jul 2015 17:25:59 +0200 Subject: Upgrade to 2Geom r2413 (bzr r14059.2.18) --- src/2geom/ellipse.cpp | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) (limited to 'src/2geom/ellipse.cpp') diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index 46c60d85d..0264ab4ae 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -97,6 +97,14 @@ void Ellipse::setCoefficients(double A, double B, double C, double D, double E, makeCanonical(); } +Point Ellipse::initialPoint() const +{ + Coord sinrot, cosrot; + sincos(_angle, sinrot, cosrot); + Point p(ray(X) * cosrot + center(X), ray(X) * sinrot + center(Y)); + return p; +} + Affine Ellipse::unitCircleTransform() const { @@ -299,8 +307,19 @@ Point Ellipse::pointAt(Coord t) const Coord Ellipse::valueAt(Coord t, Dim2 d) const { - // TODO: more efficient version. - return pointAt(t)[d]; + Coord sinrot, cosrot, cost, sint; + sincos(rotationAngle(), sinrot, cosrot); + sincos(t, sint, cost); + + if ( d == X ) { + return ray(X) * cosrot * cost + - ray(Y) * sinrot * sint + + center(X); + } else { + return ray(X) * sinrot * cost + + ray(Y) * cosrot * sint + + center(Y); + } } Coord Ellipse::timeAt(Point const &p) const @@ -319,6 +338,20 @@ Coord Ellipse::timeAt(Point const &p) const return Angle(atan2(p * iuct)).radians0(); // return a value in [0, 2pi) } +Point Ellipse::unitTangentAt(Coord t) const +{ + Point p = Point::polar(t + M_PI/2); + p *= unitCircleTransform().withoutTranslation(); + p.normalize(); + return p; +} + +bool Ellipse::contains(Point const &p) const +{ + Point tp = p * inverseUnitCircleTransform(); + return tp.length() <= 1; +} + std::vector Ellipse::intersect(Line const &line) const { -- cgit v1.2.3 From 02aedb930177c008f1074fa7b17dd6085825ee1b Mon Sep 17 00:00:00 2001 From: Krzysztof Kosi??ski Date: Tue, 7 Jul 2015 01:29:26 +0200 Subject: Update 2Geom to r2417: make scan-build happy (bzr r14235) --- src/2geom/ellipse.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/2geom/ellipse.cpp') diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index 0264ab4ae..fed3bf86f 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -499,7 +499,7 @@ std::vector Ellipse::intersect(Ellipse const &other) const Line lines[2]; if (aa != 0) { - bb /= aa; cc /= aa; dd /= aa; ee /= aa; ff /= aa; + bb /= aa; cc /= aa; dd /= aa; ee /= aa; /*ff /= aa;*/ Coord s = (ee + std::sqrt(ee*ee - 4*bb)) / 2; Coord q = ee - s; Coord alpha = (dd - cc*q) / (s - q); @@ -508,7 +508,7 @@ std::vector Ellipse::intersect(Ellipse const &other) const lines[0] = Line(1, q, alpha); lines[1] = Line(1, s, beta); } else if (bb != 0) { - cc /= bb; dd /= bb; ee /= bb; ff /= bb; + cc /= bb; /*dd /= bb;*/ ee /= bb; ff /= bb; Coord s = ee; Coord q = 0; Coord alpha = cc / ee; -- cgit v1.2.3 From 8d9336326e52ee6a181c658d57c7c3a3b776fcf8 Mon Sep 17 00:00:00 2001 From: Krzysztof Kosi??ski Date: Wed, 15 Jul 2015 20:18:11 +0200 Subject: Sync 2Geom to revision 2420. Should fix build failure on llvm-gcc 4.2 and an issue with powerstroke. (bzr r14247) --- src/2geom/ellipse.cpp | 52 +++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 40 insertions(+), 12 deletions(-) (limited to 'src/2geom/ellipse.cpp') diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index fed3bf86f..4b5ad9762 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -205,22 +205,50 @@ Ellipse::arc(Point const &ip, Point const &inner, Point const &fp) bool sweep_flag = false; // Determination of large arc flag: - // The arc is larger than half of the ellipse if the inner point - // is on the same side of the line going from the initial - // to the final point as the center of the ellipse - Point versor = fp - ip; - double sdist_c = cross(versor, _center - ip); - double sdist_inner = cross(versor, inner - ip); - - // if we have exactly half of an arc, do not set the large flag. - if (sdist_c != 0 && sgn(sdist_c) == sgn(sdist_inner)) { + // large_arc is false when the inner point is on the same side + // of the center---initial point line as the final point, AND + // is on the same side of the center---final point line as the + // initial point. + // Additionally, large_arc is always false when we have exactly + // 1/2 of an arc, i.e. the cross product of the center -> initial point + // and center -> final point vectors is zero. + // Negating the above leads to the condition for large_arc being true. + Point fv = fp - _center; + Point iv = ip - _center; + Point innerv = inner - _center; + double ifcp = cross(fv, iv); + + if (ifcp != 0 && (sgn(cross(fv, innerv)) != sgn(ifcp) || + sgn(cross(iv, innerv)) != sgn(-ifcp))) + { large_arc_flag = true; } + //cross(-iv, fv) && large_arc_flag + + // Determination of sweep flag: - // If the inner point is on the left side of the ip-fp line, - // we go in clockwise direction. - if (sdist_inner < 0) { + // For clarity, let's assume that Y grows up. Then the cross product + // is positive for points on the left side of a vector and negative + // on the right side of a vector. + // + // cross(?, v) > 0 + // o-------------------> + // cross(?, v) < 0 + // + // If the arc is small (large_arc_flag is false) and the final point + // is on the right side of the vector initial point -> center, + // we have to go in the direction of increasing angles + // (counter-clockwise) and the sweep flag is true. + // If the arc is large, the opposite is true, since we have to reach + // the final point going the long way - in the other direction. + // We can express this observation as: + // cross(_center - ip, fp - _center) < 0 xor large_arc flag + // This is equal to: + // cross(-iv, fv) < 0 xor large_arc flag + // But cross(-iv, fv) is equal to cross(fv, iv) due to antisymmetry + // of the cross product, so we end up with the condition below. + if ((ifcp < 0) ^ large_arc_flag) { sweep_flag = true; } -- cgit v1.2.3