diff options
Diffstat (limited to 'src/2geom/elliptical-arc.cpp')
| -rw-r--r-- | src/2geom/elliptical-arc.cpp | 444 |
1 files changed, 209 insertions, 235 deletions
diff --git a/src/2geom/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp index 601f0c8c8..a04b730b3 100644 --- a/src/2geom/elliptical-arc.cpp +++ b/src/2geom/elliptical-arc.cpp @@ -38,7 +38,6 @@ #include <2geom/ellipse.h> #include <2geom/elliptical-arc.h> #include <2geom/path-sink.h> -#include <2geom/poly.h> #include <2geom/sbasis-geometric.h> #include <2geom/transforms.h> #include <2geom/utils.h> @@ -93,11 +92,16 @@ namespace Geom Rect EllipticalArc::boundsExact() const { + if (isChord()) { + return chord().boundsExact(); + } + using std::swap; + // TODO: simplify / document what is going on here. double extremes[4]; double sinrot, cosrot; - sincos(_rot_angle, sinrot, cosrot); + sincos(rotationAngle(), sinrot, cosrot); extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot ); extremes[1] = extremes[0] + M_PI; @@ -132,14 +136,14 @@ Rect EllipticalArc::boundsExact() const Point EllipticalArc::pointAtAngle(Coord t) const { - Point ret = Point::polar(t) * unitCircleTransform(); + Point ret = _ellipse.pointAt(t); return ret; } Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const { Coord sinrot, cosrot, cost, sint; - sincos(_rot_angle, sinrot, cosrot); + sincos(rotationAngle(), sinrot, cosrot); sincos(t, sint, cost); if ( d == X ) { @@ -155,13 +159,22 @@ Coord EllipticalArc::valueAtAngle(Coord t, Dim2 d) const Affine EllipticalArc::unitCircleTransform() const { - Affine ret = Scale(ray(X), ray(Y)) * Rotate(_rot_angle); - ret.setTranslation(center()); + Affine ret = _ellipse.unitCircleTransform(); + return ret; +} + +Affine EllipticalArc::inverseUnitCircleTransform() const +{ + Affine ret = _ellipse.inverseUnitCircleTransform(); return ret; } std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const { + if (isChord()) { + return chord().roots(v, d); + } + std::vector<Coord> sol; if ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) { @@ -190,7 +203,7 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const for ( unsigned int dim = 0; dim < 2; ++dim ) { - if ( are_near(ray((Dim2) dim), 0) ) + if (ray((Dim2) dim) == 0) { if ( initialPoint()[d] == v && finalPoint()[d] == v ) { @@ -212,18 +225,18 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const case X: switch(dim) { - case X: ray_prj = -ray(Y) * std::sin(_rot_angle); + case X: ray_prj = -ray(Y) * std::sin(rotationAngle()); break; - case Y: ray_prj = ray(X) * std::cos(_rot_angle); + case Y: ray_prj = ray(X) * std::cos(rotationAngle()); break; } break; case Y: switch(dim) { - case X: ray_prj = ray(Y) * std::cos(_rot_angle); + case X: ray_prj = ray(Y) * std::cos(rotationAngle()); break; - case Y: ray_prj = ray(X) * std::sin(_rot_angle); + case Y: ray_prj = ray(X) * std::sin(rotationAngle()); break; } break; @@ -269,10 +282,10 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const double rotx, roty; if (d == X) { - sincos(_rot_angle, roty, rotx); + sincos(rotationAngle(), roty, rotx); roty = -roty; } else { - sincos(_rot_angle, rotx, roty); + sincos(rotationAngle(), rotx, roty); } double rxrotx = ray(X) * rotx; @@ -336,8 +349,12 @@ std::vector<Coord> EllipticalArc::roots(Coord v, Dim2 d) const // of such an angle in the cw direction Curve *EllipticalArc::derivative() const { + if (isChord()) { + return chord().derivative(); + } + EllipticalArc *result = static_cast<EllipticalArc*>(duplicate()); - result->_center[X] = result->_center[Y] = 0; + result->_ellipse.setCenter(0, 0); result->_start_angle += M_PI/2; if( !( result->_start_angle < 2*M_PI ) ) { @@ -357,13 +374,17 @@ Curve *EllipticalArc::derivative() const 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 = map_unit_interval_on_circular_arc(t, initialAngle(), finalAngle(), _sweep); std::auto_ptr<EllipticalArc> ea( static_cast<EllipticalArc*>(duplicate()) ); - ea->_center = Point(0,0); + ea->_ellipse.setCenter(0, 0); unsigned int m = std::min(nn, 4u); for ( unsigned int i = 0; i < m; ++i ) { @@ -392,6 +413,12 @@ bool EllipticalArc::containsAngle(Coord angle) const return AngleInterval::contains(angle); } +Point EllipticalArc::pointAt(Coord t) const +{ + if (isChord()) return chord().pointAt(t); + return _ellipse.pointAt(angleAt(t)); +} + Curve* EllipticalArc::portion(double f, double t) const { // fix input arguments @@ -400,14 +427,14 @@ Curve* EllipticalArc::portion(double f, double t) const if (t < 0) t = 0; if (t > 1) t = 1; - if ( are_near(f, t) ) + if (f == t) { EllipticalArc *arc = static_cast<EllipticalArc*>(duplicate()); - arc->_center = arc->_initial_point = arc->_final_point = pointAt(f); - arc->_start_angle = arc->_end_angle = _start_angle; - arc->_rot_angle = _rot_angle; - arc->_sweep = _sweep; - arc->_large_arc = _large_arc; + arc->_initial_point = arc->_final_point = pointAt(f); + arc->_start_angle = arc->_end_angle = 0; + arc->_ellipse = Ellipse(); + arc->_sweep = false; + arc->_large_arc = false; return arc; } @@ -415,21 +442,24 @@ Curve* EllipticalArc::portion(double f, double t) const arc->_initial_point = pointAt(f); arc->_final_point = pointAt(t); if ( f > t ) arc->_sweep = !_sweep; - if ( _large_arc && fabs(sweepAngle() * (t-f)) < M_PI) + if ( _large_arc && fabs(sweepAngle() * (t-f)) < M_PI) { arc->_large_arc = false; - arc->_updateCenterAndAngles(arc->isSVGCompliant()); //TODO: be more clever + } + //arc->_start_angle = angleAt(f); + //arc->_end_angle = angleAt(t); + arc->_updateCenterAndAngles(); //TODO: be more clever return arc; } // the arc is the same but traversed in the opposite direction -Curve *EllipticalArc::reverse() const { +Curve *EllipticalArc::reverse() const +{ EllipticalArc *rarc = static_cast<EllipticalArc*>(duplicate()); rarc->_sweep = !_sweep; rarc->_initial_point = _final_point; rarc->_final_point = _initial_point; rarc->_start_angle = _end_angle; rarc->_end_angle = _start_angle; - rarc->_updateCenterAndAngles(rarc->isSVGCompliant()); return rarc; } @@ -455,8 +485,8 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from, Point np = seg.pointAt( seg.nearestTime(p) ); if ( are_near(ray(Y), 0) ) { - if ( are_near(_rot_angle, M_PI/2) - || are_near(_rot_angle, 3*M_PI/2) ) + if ( are_near(rotationAngle(), M_PI/2) + || are_near(rotationAngle(), 3*M_PI/2) ) { result = roots(np[Y], Y); } @@ -467,8 +497,8 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from, } else { - if ( are_near(_rot_angle, M_PI/2) - || are_near(_rot_angle, 3*M_PI/2) ) + if ( are_near(rotationAngle(), M_PI/2) + || are_near(rotationAngle(), 3*M_PI/2) ) { result = roots(np[X], X); } @@ -527,7 +557,7 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from, Point p_c = p - center(); double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y)); double sinrot, cosrot; - sincos(_rot_angle, sinrot, cosrot); + sincos(rotationAngle(), sinrot, cosrot); double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot); Poly coeff; coeff.resize(5); @@ -662,227 +692,156 @@ std::vector<double> EllipticalArc::allNearestTimes( Point const& p, double from, } #endif -/* - * NOTE: this implementation follows Standard SVG 1.1 implementation guidelines - * for elliptical arc curves. See Appendix F.6. - */ -void EllipticalArc::_updateCenterAndAngles(bool svg) -{ - Point d = initialPoint() - finalPoint(); - // TODO move this to SVGElipticalArc? - if (svg) - { - if ( initialPoint() == finalPoint() ) - { - _rot_angle = _start_angle = _end_angle = 0; - _center = initialPoint(); - _rays = Geom::Point(0,0); - _large_arc = _sweep = false; - return; +void EllipticalArc::_filterIntersections(std::vector<ShapeIntersection> &xs, bool is_first) const +{ + Interval unit(0, 1); + std::vector<ShapeIntersection>::reverse_iterator i = xs.rbegin(), last = xs.rend(); + while (i != last) { + Coord &t = is_first ? i->first : i->second; + assert(are_near(_ellipse.pointAt(t), i->point(), 1e-6)); + t = timeAtAngle(t); + if (!unit.contains(t)) { + xs.erase((++i).base()); + continue; + } else { + assert(are_near(pointAt(t), i->point(), 1e-6)); + ++i; } + } +} - _rays[X] = std::fabs(_rays[X]); - _rays[Y] = std::fabs(_rays[Y]); - - if ( are_near(ray(X), 0) || are_near(ray(Y), 0) ) - { - _rays[X] = L2(d) / 2; - _rays[Y] = 0; - _rot_angle = std::atan2(d[Y], d[X]); - _start_angle = 0; - _end_angle = M_PI; - _center = middle_point(initialPoint(), finalPoint()); - _large_arc = false; - _sweep = false; - return; - } +std::vector<CurveIntersection> EllipticalArc::intersect(Curve const &other, Coord eps) const +{ + if (isLineSegment()) { + LineSegment ls(_initial_point, _final_point); + return ls.intersect(other, eps); } - else - { - if ( are_near(initialPoint(), finalPoint()) ) - { - if ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) - { - _start_angle = _end_angle = 0; - _center = initialPoint(); - return; - } - else - { - THROW_RANGEERROR("initial and final point are the same"); - } - } - if ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) - { // but initialPoint != finalPoint - THROW_RANGEERROR( - "there is no ellipse that satisfies the given constraints: " - "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint" - ); - } - if ( are_near(ray(Y), 0) ) - { - Point v = initialPoint() - finalPoint(); - if ( are_near(L2sq(v), 4*ray(X)*ray(X)) ) - { - Angle angle(v); - if ( are_near( angle, _rot_angle ) ) - { - _start_angle = 0; - _end_angle = M_PI; - _center = v/2 + finalPoint(); - return; - } - angle -= M_PI; - if ( are_near( angle, _rot_angle ) ) - { - _start_angle = M_PI; - _end_angle = 0; - _center = v/2 + finalPoint(); - return; - } - THROW_RANGEERROR( - "there is no ellipse that satisfies the given constraints: " - "ray(Y) == 0 " - "and slope(initialPoint - finalPoint) != rotation_angle " - "and != rotation_angle + PI" - ); - } - if ( L2sq(v) > 4*ray(X)*ray(X) ) - { - THROW_RANGEERROR( - "there is no ellipse that satisfies the given constraints: " - "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)" - ); - } - else - { - THROW_RANGEERROR( - "there is infinite ellipses that satisfy the given constraints: " - "ray(Y) == 0 and distance(initialPoint, finalPoint) < 2*ray(X)" - ); - } - } + std::vector<CurveIntersection> result; - if ( are_near(ray(X), 0) ) - { - Point v = initialPoint() - finalPoint(); - if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) ) - { - double angle = std::atan2(v[Y], v[X]); - if (angle < 0) angle += 2*M_PI; - double rot_angle = _rot_angle.radians() + M_PI/2; - if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI; - if ( are_near( angle, rot_angle ) ) - { - _start_angle = M_PI/2; - _end_angle = 3*M_PI/2; - _center = v/2 + finalPoint(); - return; - } - angle -= M_PI; - if ( angle < 0 ) angle += 2*M_PI; - if ( are_near( angle, rot_angle ) ) - { - _start_angle = 3*M_PI/2; - _end_angle = M_PI/2; - _center = v/2 + finalPoint(); - return; - } - THROW_RANGEERROR( - "there is no ellipse that satisfies the given constraints: " - "ray(X) == 0 " - "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 " - "and != rotation_angle + (3/2)*PI" - ); - } - if ( L2sq(v) > 4*ray(Y)*ray(Y) ) - { - THROW_RANGEERROR( - "there is no ellipse that satisfies the given constraints: " - "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)" - ); - } - else - { - THROW_RANGEERROR( - "there is infinite ellipses that satisfy the given constraints: " - "ray(X) == 0 and distance(initialPoint, finalPoint) < 2*ray(Y)" - ); - } + if (other.isLineSegment()) { + LineSegment ls(other.initialPoint(), other.finalPoint()); + result = _ellipse.intersect(ls); + _filterIntersections(result, true); + return result; + } - } + BezierCurve const *bez = dynamic_cast<BezierCurve const *>(&other); + if (bez) { + result = _ellipse.intersect(bez->fragment()); + _filterIntersections(result, true); + return result; + } + EllipticalArc const *arc = dynamic_cast<EllipticalArc const *>(&other); + if (arc) { + result = _ellipse.intersect(arc->_ellipse); + _filterIntersections(result, true); + arc->_filterIntersections(result, false); + return result; } - Rotate rm(_rot_angle); - Affine m(rm); - m[1] = -m[1]; - m[2] = -m[2]; - - Point p = (d / 2) * m; - double rx2 = _rays[X] * _rays[X]; - double ry2 = _rays[Y] * _rays[Y]; - double rxpy = _rays[X] * p[Y]; - double rypx = _rays[Y] * p[X]; - double rx2py2 = rxpy * rxpy; - double ry2px2 = rypx * rypx; - double num = rx2 * ry2; - double den = rx2py2 + ry2px2; - assert(den != 0); - double rad = num / den; - Point c(0,0); - if (rad > 1) - { - rad -= 1; - rad = std::sqrt(rad); + // in case someone wants to make a custom curve type + result = other.intersect(*this, eps); + transpose_in_place(result); + return result; +} - if (_large_arc == _sweep) rad = -rad; - c = rad * Point(rxpy / ray(Y), -rypx / ray(X)); - _center = c * rm + middle_point(initialPoint(), finalPoint()); +void EllipticalArc::_updateCenterAndAngles() +{ + // See: http://www.w3.org/TR/SVG/implnote.html#ArcImplementationNotes + Point d = initialPoint() - finalPoint(); + Point mid = middle_point(initialPoint(), finalPoint()); + + // if ip = sp, the arc contains no other points + if (initialPoint() == finalPoint()) { + _start_angle = _end_angle = 0; + _ellipse.setRotationAngle(0); + _ellipse.setCenter(initialPoint()); + _ellipse.setRays(0, 0); + _large_arc = _sweep = false; + return; } - else if (rad == 1 || svg) - { - double lamda = std::sqrt(1 / rad); - _rays[X] *= lamda; - _rays[Y] *= lamda; - _center = middle_point(initialPoint(), finalPoint()); + + // rays should be positive + _ellipse.setRays(std::fabs(ray(X)), std::fabs(ray(Y))); + + if (ray(X) == 0 || ray(Y) == 0) { + _start_angle = 0; + _end_angle = M_PI; + _ellipse.setRays(L2(d) / 2, 0); + _ellipse.setRotationAngle(atan2(d)); + _ellipse.setCenter(mid); + _large_arc = false; + _sweep = false; + return; } - else - { - THROW_RANGEERROR( - "there is no ellipse that satisfies the given constraints" - ); + + 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 = (initialPoint() - mid) * 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 rad = (rxry - pxry - rxpy)/(rxpy + pxry); + // 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); + } } - Point sp((p[X] - c[X]) / ray(X), (p[Y] - c[Y]) / ray(Y)); - Point ep((-p[X] - c[X]) / ray(X), (-p[Y] - c[Y]) / ray(Y)); + // 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); _start_angle = angle_between(v, sp); 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; - _end_angle = _start_angle; - _end_angle += sweep_angle; + _end_angle = _start_angle + sweep_angle; } 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() + ( _sweep ? sweepAngle() : -sweepAngle() ); - Linear param(initialAngle(), et); - Coord cos_rot_angle, sin_rot_angle; - sincos(_rot_angle, sin_rot_angle, cos_rot_angle); + 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 * 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)); + 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++ ) { @@ -898,22 +857,24 @@ void EllipticalArc::transform(Affine const& m) if (isChord()) { _initial_point *= m; _final_point *= m; - _center = middle_point(_initial_point, _final_point); - _rays = Point(0,0); - _rot_angle = 0; + _ellipse.setCenter(middle_point(_initial_point, _final_point)); + _ellipse.setRays(0, 0); + _ellipse.setRotationAngle(0); return; } - // TODO avoid allocating a new arc here - Ellipse e(center(X), center(Y), ray(X), ray(Y), _rot_angle); - e *= m; - Point inner_point = pointAt(0.5); - EllipticalArc *arc = e.arc(initialPoint() * m, - inner_point * m, - finalPoint() * m, - isSVGCompliant()); - *this = *arc; - delete arc; + _initial_point *= m; + _final_point *= m; + _ellipse *= m; + if (m.det() < 0) { + _sweep = !_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. + _start_angle = _ellipse.timeAt(_initial_point); + _end_angle = _ellipse.timeAt(_final_point); } bool EllipticalArc::operator==(Curve const &c) const @@ -924,8 +885,8 @@ bool EllipticalArc::operator==(Curve const &c) const 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 (_rot_angle != other->_rot_angle) return false; + 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; @@ -936,7 +897,7 @@ void EllipticalArc::feed(PathSink &sink, bool moveto_initial) const if (moveto_initial) { sink.moveTo(_initial_point); } - sink.arcTo(_rays[X], _rays[Y], _rot_angle, _large_arc, _sweep, _final_point); + sink.arcTo(ray(X), ray(Y), rotationAngle(), _large_arc, _sweep, _final_point); } Coord EllipticalArc::map_to_01(Coord angle) const @@ -945,6 +906,19 @@ Coord EllipticalArc::map_to_01(Coord angle) const finalAngle(), _sweep); } + +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 /* |