michael@0: /* -*- Mode: C++; tab-width: 20; indent-tabs-mode: nil; c-basic-offset: 2 -*- michael@0: * This Source Code Form is subject to the terms of the Mozilla Public michael@0: * License, v. 2.0. If a copy of the MPL was not distributed with this michael@0: * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ michael@0: michael@0: #include "2D.h" michael@0: #include "PathAnalysis.h" michael@0: #include "PathHelpers.h" michael@0: michael@0: namespace mozilla { michael@0: namespace gfx { michael@0: michael@0: static float CubicRoot(float aValue) { michael@0: if (aValue < 0.0) { michael@0: return -CubicRoot(-aValue); michael@0: } michael@0: else { michael@0: return powf(aValue, 1.0f / 3.0f); michael@0: } michael@0: } michael@0: michael@0: struct BezierControlPoints michael@0: { michael@0: BezierControlPoints() {} michael@0: BezierControlPoints(const Point &aCP1, const Point &aCP2, michael@0: const Point &aCP3, const Point &aCP4) michael@0: : mCP1(aCP1), mCP2(aCP2), mCP3(aCP3), mCP4(aCP4) michael@0: { michael@0: } michael@0: michael@0: Point mCP1, mCP2, mCP3, mCP4; michael@0: }; michael@0: michael@0: void michael@0: FlattenBezier(const BezierControlPoints &aPoints, michael@0: PathSink *aSink, Float aTolerance); michael@0: michael@0: michael@0: Path::Path() michael@0: { michael@0: } michael@0: michael@0: Path::~Path() michael@0: { michael@0: } michael@0: michael@0: Float michael@0: Path::ComputeLength() michael@0: { michael@0: EnsureFlattenedPath(); michael@0: return mFlattenedPath->ComputeLength(); michael@0: } michael@0: michael@0: Point michael@0: Path::ComputePointAtLength(Float aLength, Point* aTangent) michael@0: { michael@0: EnsureFlattenedPath(); michael@0: return mFlattenedPath->ComputePointAtLength(aLength, aTangent); michael@0: } michael@0: michael@0: void michael@0: Path::EnsureFlattenedPath() michael@0: { michael@0: if (!mFlattenedPath) { michael@0: mFlattenedPath = new FlattenedPath(); michael@0: StreamToSink(mFlattenedPath); michael@0: } michael@0: } michael@0: michael@0: // This is the maximum deviation we allow (with an additional ~20% margin of michael@0: // error) of the approximation from the actual Bezier curve. michael@0: const Float kFlatteningTolerance = 0.0001f; michael@0: michael@0: void michael@0: FlattenedPath::MoveTo(const Point &aPoint) michael@0: { michael@0: MOZ_ASSERT(!mCalculatedLength); michael@0: FlatPathOp op; michael@0: op.mType = FlatPathOp::OP_MOVETO; michael@0: op.mPoint = aPoint; michael@0: mPathOps.push_back(op); michael@0: michael@0: mLastMove = aPoint; michael@0: } michael@0: michael@0: void michael@0: FlattenedPath::LineTo(const Point &aPoint) michael@0: { michael@0: MOZ_ASSERT(!mCalculatedLength); michael@0: FlatPathOp op; michael@0: op.mType = FlatPathOp::OP_LINETO; michael@0: op.mPoint = aPoint; michael@0: mPathOps.push_back(op); michael@0: } michael@0: michael@0: void michael@0: FlattenedPath::BezierTo(const Point &aCP1, michael@0: const Point &aCP2, michael@0: const Point &aCP3) michael@0: { michael@0: MOZ_ASSERT(!mCalculatedLength); michael@0: FlattenBezier(BezierControlPoints(CurrentPoint(), aCP1, aCP2, aCP3), this, kFlatteningTolerance); michael@0: } michael@0: michael@0: void michael@0: FlattenedPath::QuadraticBezierTo(const Point &aCP1, michael@0: const Point &aCP2) michael@0: { michael@0: MOZ_ASSERT(!mCalculatedLength); michael@0: // We need to elevate the degree of this quadratic Bézier to cubic, so we're michael@0: // going to add an intermediate control point, and recompute control point 1. michael@0: // The first and last control points remain the same. michael@0: // This formula can be found on http://fontforge.sourceforge.net/bezier.html michael@0: Point CP0 = CurrentPoint(); michael@0: Point CP1 = (CP0 + aCP1 * 2.0) / 3.0; michael@0: Point CP2 = (aCP2 + aCP1 * 2.0) / 3.0; michael@0: Point CP3 = aCP2; michael@0: michael@0: BezierTo(CP1, CP2, CP3); michael@0: } michael@0: michael@0: void michael@0: FlattenedPath::Close() michael@0: { michael@0: MOZ_ASSERT(!mCalculatedLength); michael@0: LineTo(mLastMove); michael@0: } michael@0: michael@0: void michael@0: FlattenedPath::Arc(const Point &aOrigin, float aRadius, float aStartAngle, michael@0: float aEndAngle, bool aAntiClockwise) michael@0: { michael@0: ArcToBezier(this, aOrigin, Size(aRadius, aRadius), aStartAngle, aEndAngle, aAntiClockwise); michael@0: } michael@0: michael@0: Float michael@0: FlattenedPath::ComputeLength() michael@0: { michael@0: if (!mCalculatedLength) { michael@0: Point currentPoint; michael@0: michael@0: for (uint32_t i = 0; i < mPathOps.size(); i++) { michael@0: if (mPathOps[i].mType == FlatPathOp::OP_MOVETO) { michael@0: currentPoint = mPathOps[i].mPoint; michael@0: } else { michael@0: mCachedLength += Distance(currentPoint, mPathOps[i].mPoint); michael@0: currentPoint = mPathOps[i].mPoint; michael@0: } michael@0: } michael@0: michael@0: mCalculatedLength = true; michael@0: } michael@0: michael@0: return mCachedLength; michael@0: } michael@0: michael@0: Point michael@0: FlattenedPath::ComputePointAtLength(Float aLength, Point *aTangent) michael@0: { michael@0: // We track the last point that -wasn't- in the same place as the current michael@0: // point so if we pass the edge of the path with a bunch of zero length michael@0: // paths we still get the correct tangent vector. michael@0: Point lastPointSinceMove; michael@0: Point currentPoint; michael@0: for (uint32_t i = 0; i < mPathOps.size(); i++) { michael@0: if (mPathOps[i].mType == FlatPathOp::OP_MOVETO) { michael@0: if (Distance(currentPoint, mPathOps[i].mPoint)) { michael@0: lastPointSinceMove = currentPoint; michael@0: } michael@0: currentPoint = mPathOps[i].mPoint; michael@0: } else { michael@0: Float segmentLength = Distance(currentPoint, mPathOps[i].mPoint); michael@0: michael@0: if (segmentLength) { michael@0: lastPointSinceMove = currentPoint; michael@0: if (segmentLength > aLength) { michael@0: Point currentVector = mPathOps[i].mPoint - currentPoint; michael@0: Point tangent = currentVector / segmentLength; michael@0: if (aTangent) { michael@0: *aTangent = tangent; michael@0: } michael@0: return currentPoint + tangent * aLength; michael@0: } michael@0: } michael@0: michael@0: aLength -= segmentLength; michael@0: currentPoint = mPathOps[i].mPoint; michael@0: } michael@0: } michael@0: michael@0: Point currentVector = currentPoint - lastPointSinceMove; michael@0: if (aTangent) { michael@0: if (hypotf(currentVector.x, currentVector.y)) { michael@0: *aTangent = currentVector / hypotf(currentVector.x, currentVector.y); michael@0: } else { michael@0: *aTangent = Point(); michael@0: } michael@0: } michael@0: return currentPoint; michael@0: } michael@0: michael@0: // This function explicitly permits aControlPoints to refer to the same object michael@0: // as either of the other arguments. michael@0: static void michael@0: SplitBezier(const BezierControlPoints &aControlPoints, michael@0: BezierControlPoints *aFirstSegmentControlPoints, michael@0: BezierControlPoints *aSecondSegmentControlPoints, michael@0: Float t) michael@0: { michael@0: MOZ_ASSERT(aSecondSegmentControlPoints); michael@0: michael@0: *aSecondSegmentControlPoints = aControlPoints; michael@0: michael@0: Point cp1a = aControlPoints.mCP1 + (aControlPoints.mCP2 - aControlPoints.mCP1) * t; michael@0: Point cp2a = aControlPoints.mCP2 + (aControlPoints.mCP3 - aControlPoints.mCP2) * t; michael@0: Point cp1aa = cp1a + (cp2a - cp1a) * t; michael@0: Point cp3a = aControlPoints.mCP3 + (aControlPoints.mCP4 - aControlPoints.mCP3) * t; michael@0: Point cp2aa = cp2a + (cp3a - cp2a) * t; michael@0: Point cp1aaa = cp1aa + (cp2aa - cp1aa) * t; michael@0: aSecondSegmentControlPoints->mCP4 = aControlPoints.mCP4; michael@0: michael@0: if(aFirstSegmentControlPoints) { michael@0: aFirstSegmentControlPoints->mCP1 = aControlPoints.mCP1; michael@0: aFirstSegmentControlPoints->mCP2 = cp1a; michael@0: aFirstSegmentControlPoints->mCP3 = cp1aa; michael@0: aFirstSegmentControlPoints->mCP4 = cp1aaa; michael@0: } michael@0: aSecondSegmentControlPoints->mCP1 = cp1aaa; michael@0: aSecondSegmentControlPoints->mCP2 = cp2aa; michael@0: aSecondSegmentControlPoints->mCP3 = cp3a; michael@0: } michael@0: michael@0: static void michael@0: FlattenBezierCurveSegment(const BezierControlPoints &aControlPoints, michael@0: PathSink *aSink, michael@0: Float aTolerance) michael@0: { michael@0: /* The algorithm implemented here is based on: michael@0: * http://cis.usouthal.edu/~hain/general/Publications/Bezier/Bezier%20Offset%20Curves.pdf michael@0: * michael@0: * The basic premise is that for a small t the third order term in the michael@0: * equation of a cubic bezier curve is insignificantly small. This can michael@0: * then be approximated by a quadratic equation for which the maximum michael@0: * difference from a linear approximation can be much more easily determined. michael@0: */ michael@0: BezierControlPoints currentCP = aControlPoints; michael@0: michael@0: Float t = 0; michael@0: while (t < 1.0f) { michael@0: Point cp21 = currentCP.mCP2 - currentCP.mCP3; michael@0: Point cp31 = currentCP.mCP3 - currentCP.mCP1; michael@0: michael@0: Float s3 = (cp31.x * cp21.y - cp31.y * cp21.x) / hypotf(cp21.x, cp21.y); michael@0: michael@0: t = 2 * Float(sqrt(aTolerance / (3. * abs(s3)))); michael@0: michael@0: if (t >= 1.0f) { michael@0: aSink->LineTo(aControlPoints.mCP4); michael@0: break; michael@0: } michael@0: michael@0: Point prevCP2, prevCP3, nextCP1, nextCP2, nextCP3; michael@0: SplitBezier(currentCP, nullptr, ¤tCP, t); michael@0: michael@0: aSink->LineTo(currentCP.mCP1); michael@0: } michael@0: } michael@0: michael@0: static inline void michael@0: FindInflectionApproximationRange(BezierControlPoints aControlPoints, michael@0: Float *aMin, Float *aMax, Float aT, michael@0: Float aTolerance) michael@0: { michael@0: SplitBezier(aControlPoints, nullptr, &aControlPoints, aT); michael@0: michael@0: Point cp21 = aControlPoints.mCP2 - aControlPoints.mCP1; michael@0: Point cp41 = aControlPoints.mCP4 - aControlPoints.mCP1; michael@0: michael@0: if (cp21.x == 0 && cp21.y == 0) { michael@0: // In this case s3 becomes lim[n->0] (cp41.x * n) / n - (cp41.y * n) / n = cp41.x - cp41.y. michael@0: michael@0: // Use the absolute value so that Min and Max will correspond with the michael@0: // minimum and maximum of the range. michael@0: *aMin = aT - CubicRoot(abs(aTolerance / (cp41.x - cp41.y))); michael@0: *aMax = aT + CubicRoot(abs(aTolerance / (cp41.x - cp41.y))); michael@0: return; michael@0: } michael@0: michael@0: Float s3 = (cp41.x * cp21.y - cp41.y * cp21.x) / hypotf(cp21.x, cp21.y); michael@0: michael@0: if (s3 == 0) { michael@0: // This means within the precision we have it can be approximated michael@0: // infinitely by a linear segment. Deal with this by specifying the michael@0: // approximation range as extending beyond the entire curve. michael@0: *aMin = -1.0f; michael@0: *aMax = 2.0f; michael@0: return; michael@0: } michael@0: michael@0: Float tf = CubicRoot(abs(aTolerance / s3)); michael@0: michael@0: *aMin = aT - tf * (1 - aT); michael@0: *aMax = aT + tf * (1 - aT); michael@0: } michael@0: michael@0: /* Find the inflection points of a bezier curve. Will return false if the michael@0: * curve is degenerate in such a way that it is best approximated by a straight michael@0: * line. michael@0: * michael@0: * The below algorithm was written by Jeff Muizelaar , explanation follows: michael@0: * michael@0: * The lower inflection point is returned in aT1, the higher one in aT2. In the michael@0: * case of a single inflection point this will be in aT1. michael@0: * michael@0: * The method is inspired by the algorithm in "analysis of in?ection points for planar cubic bezier curve" michael@0: * michael@0: * Here are some differences between this algorithm and versions discussed elsewhere in the literature: michael@0: * michael@0: * zhang et. al compute a0, d0 and e0 incrementally using the follow formula: michael@0: * michael@0: * Point a0 = CP2 - CP1 michael@0: * Point a1 = CP3 - CP2 michael@0: * Point a2 = CP4 - CP1 michael@0: * michael@0: * Point d0 = a1 - a0 michael@0: * Point d1 = a2 - a1 michael@0: michael@0: * Point e0 = d1 - d0 michael@0: * michael@0: * this avoids any multiplications and may or may not be faster than the approach take below. michael@0: * michael@0: * "fast, precise flattening of cubic bezier path and ofset curves" by hain et. al michael@0: * Point a = CP1 + 3 * CP2 - 3 * CP3 + CP4 michael@0: * Point b = 3 * CP1 - 6 * CP2 + 3 * CP3 michael@0: * Point c = -3 * CP1 + 3 * CP2 michael@0: * Point d = CP1 michael@0: * the a, b, c, d can be expressed in terms of a0, d0 and e0 defined above as: michael@0: * c = 3 * a0 michael@0: * b = 3 * d0 michael@0: * a = e0 michael@0: * michael@0: * michael@0: * a = 3a = a.y * b.x - a.x * b.y michael@0: * b = 3b = a.y * c.x - a.x * c.y michael@0: * c = 9c = b.y * c.x - b.x * c.y michael@0: * michael@0: * The additional multiples of 3 cancel each other out as show below: michael@0: * michael@0: * x = (-b + sqrt(b * b - 4 * a * c)) / (2 * a) michael@0: * x = (-3 * b + sqrt(3 * b * 3 * b - 4 * a * 3 * 9 * c / 3)) / (2 * 3 * a) michael@0: * x = 3 * (-b + sqrt(b * b - 4 * a * c)) / (2 * 3 * a) michael@0: * x = (-b + sqrt(b * b - 4 * a * c)) / (2 * a) michael@0: * michael@0: * I haven't looked into whether the formulation of the quadratic formula in michael@0: * hain has any numerical advantages over the one used below. michael@0: */ michael@0: static inline void michael@0: FindInflectionPoints(const BezierControlPoints &aControlPoints, michael@0: Float *aT1, Float *aT2, uint32_t *aCount) michael@0: { michael@0: // Find inflection points. michael@0: // See www.faculty.idc.ac.il/arik/quality/appendixa.html for an explanation michael@0: // of this approach. michael@0: Point A = aControlPoints.mCP2 - aControlPoints.mCP1; michael@0: Point B = aControlPoints.mCP3 - (aControlPoints.mCP2 * 2) + aControlPoints.mCP1; michael@0: Point C = aControlPoints.mCP4 - (aControlPoints.mCP3 * 3) + (aControlPoints.mCP2 * 3) - aControlPoints.mCP1; michael@0: michael@0: Float a = Float(B.x) * C.y - Float(B.y) * C.x; michael@0: Float b = Float(A.x) * C.y - Float(A.y) * C.x; michael@0: Float c = Float(A.x) * B.y - Float(A.y) * B.x; michael@0: michael@0: if (a == 0) { michael@0: // Not a quadratic equation. michael@0: if (b == 0) { michael@0: // Instead of a linear acceleration change we have a constant michael@0: // acceleration change. This means the equation has no solution michael@0: // and there are no inflection points, unless the constant is 0. michael@0: // In that case the curve is a straight line, but we'll let michael@0: // FlattenBezierCurveSegment deal with this. michael@0: *aCount = 0; michael@0: return; michael@0: } michael@0: *aT1 = -c / b; michael@0: *aCount = 1; michael@0: return; michael@0: } else { michael@0: Float discriminant = b * b - 4 * a * c; michael@0: michael@0: if (discriminant < 0) { michael@0: // No inflection points. michael@0: *aCount = 0; michael@0: } else if (discriminant == 0) { michael@0: *aCount = 1; michael@0: *aT1 = -b / (2 * a); michael@0: } else { michael@0: /* Use the following formula for computing the roots: michael@0: * michael@0: * q = -1/2 * (b + sign(b) * sqrt(b^2 - 4ac)) michael@0: * t1 = q / a michael@0: * t2 = c / q michael@0: */ michael@0: Float q = sqrtf(discriminant); michael@0: if (b < 0) { michael@0: q = b - q; michael@0: } else { michael@0: q = b + q; michael@0: } michael@0: q *= Float(-1./2); michael@0: michael@0: *aT1 = q / a; michael@0: *aT2 = c / q; michael@0: if (*aT1 > *aT2) { michael@0: std::swap(*aT1, *aT2); michael@0: } michael@0: *aCount = 2; michael@0: } michael@0: } michael@0: michael@0: return; michael@0: } michael@0: michael@0: void michael@0: FlattenBezier(const BezierControlPoints &aControlPoints, michael@0: PathSink *aSink, Float aTolerance) michael@0: { michael@0: Float t1; michael@0: Float t2; michael@0: uint32_t count; michael@0: michael@0: FindInflectionPoints(aControlPoints, &t1, &t2, &count); michael@0: michael@0: // Check that at least one of the inflection points is inside [0..1] michael@0: if (count == 0 || ((t1 < 0 || t1 > 1.0) && ((t2 < 0 || t2 > 1.0) || count == 1)) ) { michael@0: FlattenBezierCurveSegment(aControlPoints, aSink, aTolerance); michael@0: return; michael@0: } michael@0: michael@0: Float t1min = t1, t1max = t1, t2min = t2, t2max = t2; michael@0: michael@0: BezierControlPoints remainingCP = aControlPoints; michael@0: michael@0: // For both inflection points, calulate the range where they can be linearly michael@0: // approximated if they are positioned within [0,1] michael@0: if (count > 0 && t1 >= 0 && t1 < 1.0) { michael@0: FindInflectionApproximationRange(aControlPoints, &t1min, &t1max, t1, aTolerance); michael@0: } michael@0: if (count > 1 && t2 >= 0 && t2 < 1.0) { michael@0: FindInflectionApproximationRange(aControlPoints, &t2min, &t2max, t2, aTolerance); michael@0: } michael@0: BezierControlPoints nextCPs = aControlPoints; michael@0: BezierControlPoints prevCPs; michael@0: michael@0: // Process ranges. [t1min, t1max] and [t2min, t2max] are approximated by line michael@0: // segments. michael@0: if (t1min > 0) { michael@0: // Flatten the Bezier up until the first inflection point's approximation michael@0: // point. michael@0: SplitBezier(aControlPoints, &prevCPs, michael@0: &remainingCP, t1min); michael@0: FlattenBezierCurveSegment(prevCPs, aSink, aTolerance); michael@0: } michael@0: if (t1max >= 0 && t1max < 1.0 && (count == 1 || t2min > t1max)) { michael@0: // The second inflection point's approximation range begins after the end michael@0: // of the first, approximate the first inflection point by a line and michael@0: // subsequently flatten up until the end or the next inflection point. michael@0: SplitBezier(aControlPoints, nullptr, &nextCPs, t1max); michael@0: michael@0: aSink->LineTo(nextCPs.mCP1); michael@0: michael@0: if (count == 1 || (count > 1 && t2min >= 1.0)) { michael@0: // No more inflection points to deal with, flatten the rest of the curve. michael@0: FlattenBezierCurveSegment(nextCPs, aSink, aTolerance); michael@0: } michael@0: } else if (count > 1 && t2min > 1.0) { michael@0: // We've already concluded t2min <= t1max, so if this is true the michael@0: // approximation range for the first inflection point runs past the michael@0: // end of the curve, draw a line to the end and we're done. michael@0: aSink->LineTo(aControlPoints.mCP4); michael@0: return; michael@0: } michael@0: michael@0: if (count > 1 && t2min < 1.0 && t2max > 0) { michael@0: if (t2min > 0 && t2min < t1max) { michael@0: // In this case the t2 approximation range starts inside the t1 michael@0: // approximation range. michael@0: SplitBezier(aControlPoints, nullptr, &nextCPs, t1max); michael@0: aSink->LineTo(nextCPs.mCP1); michael@0: } else if (t2min > 0 && t1max > 0) { michael@0: SplitBezier(aControlPoints, nullptr, &nextCPs, t1max); michael@0: michael@0: // Find a control points describing the portion of the curve between t1max and t2min. michael@0: Float t2mina = (t2min - t1max) / (1 - t1max); michael@0: SplitBezier(nextCPs, &prevCPs, &nextCPs, t2mina); michael@0: FlattenBezierCurveSegment(prevCPs, aSink, aTolerance); michael@0: } else if (t2min > 0) { michael@0: // We have nothing interesting before t2min, find that bit and flatten it. michael@0: SplitBezier(aControlPoints, &prevCPs, &nextCPs, t2min); michael@0: FlattenBezierCurveSegment(prevCPs, aSink, aTolerance); michael@0: } michael@0: if (t2max < 1.0) { michael@0: // Flatten the portion of the curve after t2max michael@0: SplitBezier(aControlPoints, nullptr, &nextCPs, t2max); michael@0: michael@0: // Draw a line to the start, this is the approximation between t2min and michael@0: // t2max. michael@0: aSink->LineTo(nextCPs.mCP1); michael@0: FlattenBezierCurveSegment(nextCPs, aSink, aTolerance); michael@0: } else { michael@0: // Our approximation range extends beyond the end of the curve. michael@0: aSink->LineTo(aControlPoints.mCP4); michael@0: return; michael@0: } michael@0: } michael@0: } michael@0: michael@0: } michael@0: }