1 /* -*- Mode: C++; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
2 /* vim: set ts=8 sts=2 et sw=2 tw=80: */
3 /* This Source Code Form is subject to the terms of the Mozilla Public
4 * License, v. 2.0. If a copy of the MPL was not distributed with this
5 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
8 #include "PathAnalysis.h"
9 #include "PathHelpers.h"
14 static double CubicRoot(double aValue
) {
16 return -CubicRoot(-aValue
);
18 return pow(aValue
, 1.0 / 3.0);
22 struct PointD
: public BasePoint
<double, PointD
> {
23 typedef BasePoint
<double, PointD
> Super
;
26 PointD(double aX
, double aY
) : Super(aX
, aY
) {}
27 MOZ_IMPLICIT
PointD(const Point
& aPoint
) : Super(aPoint
.x
, aPoint
.y
) {}
29 Point
ToPoint() const {
30 return Point(static_cast<Float
>(x
), static_cast<Float
>(y
));
34 struct BezierControlPoints
{
35 BezierControlPoints() = default;
36 BezierControlPoints(const PointD
& aCP1
, const PointD
& aCP2
,
37 const PointD
& aCP3
, const PointD
& aCP4
)
38 : mCP1(aCP1
), mCP2(aCP2
), mCP3(aCP3
), mCP4(aCP4
) {}
40 PointD mCP1
, mCP2
, mCP3
, mCP4
;
43 void FlattenBezier(const BezierControlPoints
& aPoints
, PathSink
* aSink
,
46 Path::Path() = default;
48 Path::~Path() = default;
50 Float
Path::ComputeLength() {
51 EnsureFlattenedPath();
52 return mFlattenedPath
->ComputeLength();
55 Point
Path::ComputePointAtLength(Float aLength
, Point
* aTangent
) {
56 EnsureFlattenedPath();
57 return mFlattenedPath
->ComputePointAtLength(aLength
, aTangent
);
60 void Path::EnsureFlattenedPath() {
61 if (!mFlattenedPath
) {
62 mFlattenedPath
= new FlattenedPath();
63 StreamToSink(mFlattenedPath
);
67 // This is the maximum deviation we allow (with an additional ~20% margin of
68 // error) of the approximation from the actual Bezier curve.
69 const Float kFlatteningTolerance
= 0.0001f
;
71 void FlattenedPath::MoveTo(const Point
& aPoint
) {
72 MOZ_ASSERT(!mCalculatedLength
);
74 op
.mType
= FlatPathOp::OP_MOVETO
;
76 mPathOps
.push_back(op
);
81 void FlattenedPath::LineTo(const Point
& aPoint
) {
82 MOZ_ASSERT(!mCalculatedLength
);
84 op
.mType
= FlatPathOp::OP_LINETO
;
86 mPathOps
.push_back(op
);
89 void FlattenedPath::BezierTo(const Point
& aCP1
, const Point
& aCP2
,
91 MOZ_ASSERT(!mCalculatedLength
);
92 FlattenBezier(BezierControlPoints(CurrentPoint(), aCP1
, aCP2
, aCP3
), this,
93 kFlatteningTolerance
);
96 void FlattenedPath::QuadraticBezierTo(const Point
& aCP1
, const Point
& aCP2
) {
97 MOZ_ASSERT(!mCalculatedLength
);
98 // We need to elevate the degree of this quadratic B�zier to cubic, so we're
99 // going to add an intermediate control point, and recompute control point 1.
100 // The first and last control points remain the same.
101 // This formula can be found on http://fontforge.sourceforge.net/bezier.html
102 Point CP0
= CurrentPoint();
103 Point CP1
= (CP0
+ aCP1
* 2.0) / 3.0;
104 Point CP2
= (aCP2
+ aCP1
* 2.0) / 3.0;
107 BezierTo(CP1
, CP2
, CP3
);
110 void FlattenedPath::Close() {
111 MOZ_ASSERT(!mCalculatedLength
);
115 void FlattenedPath::Arc(const Point
& aOrigin
, float aRadius
, float aStartAngle
,
116 float aEndAngle
, bool aAntiClockwise
) {
117 ArcToBezier(this, aOrigin
, Size(aRadius
, aRadius
), aStartAngle
, aEndAngle
,
121 Float
FlattenedPath::ComputeLength() {
122 if (!mCalculatedLength
) {
125 for (uint32_t i
= 0; i
< mPathOps
.size(); i
++) {
126 if (mPathOps
[i
].mType
== FlatPathOp::OP_MOVETO
) {
127 currentPoint
= mPathOps
[i
].mPoint
;
129 mCachedLength
+= Distance(currentPoint
, mPathOps
[i
].mPoint
);
130 currentPoint
= mPathOps
[i
].mPoint
;
134 mCalculatedLength
= true;
137 return mCachedLength
;
140 Point
FlattenedPath::ComputePointAtLength(Float aLength
, Point
* aTangent
) {
141 if (aLength
< mCursor
.mLength
) {
142 // If cursor is beyond the target length, reset to the beginning.
145 // Adjust aLength to account for the position where we'll start searching.
146 aLength
-= mCursor
.mLength
;
149 while (mCursor
.mIndex
< mPathOps
.size()) {
150 const auto& op
= mPathOps
[mCursor
.mIndex
];
151 if (op
.mType
== FlatPathOp::OP_MOVETO
) {
152 if (Distance(mCursor
.mCurrentPoint
, op
.mPoint
) > 0.0f
) {
153 mCursor
.mLastPointSinceMove
= mCursor
.mCurrentPoint
;
155 mCursor
.mCurrentPoint
= op
.mPoint
;
157 Float segmentLength
= Distance(mCursor
.mCurrentPoint
, op
.mPoint
);
160 mCursor
.mLastPointSinceMove
= mCursor
.mCurrentPoint
;
161 if (segmentLength
> aLength
) {
162 Point currentVector
= op
.mPoint
- mCursor
.mCurrentPoint
;
163 Point tangent
= currentVector
/ segmentLength
;
167 return mCursor
.mCurrentPoint
+ tangent
* aLength
;
171 aLength
-= segmentLength
;
172 mCursor
.mLength
+= segmentLength
;
173 mCursor
.mCurrentPoint
= op
.mPoint
;
179 Point currentVector
= mCursor
.mCurrentPoint
- mCursor
.mLastPointSinceMove
;
180 if (auto h
= hypotf(currentVector
.x
, currentVector
.y
)) {
181 *aTangent
= currentVector
/ h
;
186 return mCursor
.mCurrentPoint
;
189 // This function explicitly permits aControlPoints to refer to the same object
190 // as either of the other arguments.
191 static void SplitBezier(const BezierControlPoints
& aControlPoints
,
192 BezierControlPoints
* aFirstSegmentControlPoints
,
193 BezierControlPoints
* aSecondSegmentControlPoints
,
195 MOZ_ASSERT(aSecondSegmentControlPoints
);
197 *aSecondSegmentControlPoints
= aControlPoints
;
200 aControlPoints
.mCP1
+ (aControlPoints
.mCP2
- aControlPoints
.mCP1
) * t
;
202 aControlPoints
.mCP2
+ (aControlPoints
.mCP3
- aControlPoints
.mCP2
) * t
;
203 PointD cp1aa
= cp1a
+ (cp2a
- cp1a
) * t
;
205 aControlPoints
.mCP3
+ (aControlPoints
.mCP4
- aControlPoints
.mCP3
) * t
;
206 PointD cp2aa
= cp2a
+ (cp3a
- cp2a
) * t
;
207 PointD cp1aaa
= cp1aa
+ (cp2aa
- cp1aa
) * t
;
208 aSecondSegmentControlPoints
->mCP4
= aControlPoints
.mCP4
;
210 if (aFirstSegmentControlPoints
) {
211 aFirstSegmentControlPoints
->mCP1
= aControlPoints
.mCP1
;
212 aFirstSegmentControlPoints
->mCP2
= cp1a
;
213 aFirstSegmentControlPoints
->mCP3
= cp1aa
;
214 aFirstSegmentControlPoints
->mCP4
= cp1aaa
;
216 aSecondSegmentControlPoints
->mCP1
= cp1aaa
;
217 aSecondSegmentControlPoints
->mCP2
= cp2aa
;
218 aSecondSegmentControlPoints
->mCP3
= cp3a
;
221 static void FlattenBezierCurveSegment(const BezierControlPoints
& aControlPoints
,
222 PathSink
* aSink
, double aTolerance
) {
223 /* The algorithm implemented here is based on:
224 * http://cis.usouthal.edu/~hain/general/Publications/Bezier/Bezier%20Offset%20Curves.pdf
226 * The basic premise is that for a small t the third order term in the
227 * equation of a cubic bezier curve is insignificantly small. This can
228 * then be approximated by a quadratic equation for which the maximum
229 * difference from a linear approximation can be much more easily determined.
231 BezierControlPoints currentCP
= aControlPoints
;
234 double currentTolerance
= aTolerance
;
236 PointD cp21
= currentCP
.mCP2
- currentCP
.mCP1
;
237 PointD cp31
= currentCP
.mCP3
- currentCP
.mCP1
;
239 /* To remove divisions and check for divide-by-zero, this is optimized from:
240 * Float s3 = (cp31.x * cp21.y - cp31.y * cp21.x) / hypotf(cp21.x, cp21.y);
241 * t = 2 * Float(sqrt(aTolerance / (3. * std::abs(s3))));
243 double cp21x31
= cp31
.x
* cp21
.y
- cp31
.y
* cp21
.x
;
244 double h
= hypot(cp21
.x
, cp21
.y
);
245 if (cp21x31
* h
== 0) {
249 double s3inv
= h
/ cp21x31
;
250 t
= 2 * sqrt(currentTolerance
* std::abs(s3inv
) / 3.);
251 currentTolerance
*= 1 + aTolerance
;
252 // Increase tolerance every iteration to prevent this loop from executing
253 // too many times. This approximates the length of large curves more
254 // roughly. In practice, aTolerance is the constant kFlatteningTolerance
255 // which has value 0.0001. With this value, it takes 6,932 splits to double
256 // currentTolerance (to 0.0002) and 23,028 splits to increase
257 // currentTolerance by an order of magnitude (to 0.001).
262 SplitBezier(currentCP
, nullptr, ¤tCP
, t
);
264 aSink
->LineTo(currentCP
.mCP1
.ToPoint());
267 aSink
->LineTo(currentCP
.mCP4
.ToPoint());
270 static inline void FindInflectionApproximationRange(
271 BezierControlPoints aControlPoints
, double* aMin
, double* aMax
, double aT
,
273 SplitBezier(aControlPoints
, nullptr, &aControlPoints
, aT
);
275 PointD cp21
= aControlPoints
.mCP2
- aControlPoints
.mCP1
;
276 PointD cp41
= aControlPoints
.mCP4
- aControlPoints
.mCP1
;
278 if (cp21
.x
== 0. && cp21
.y
== 0.) {
279 cp21
= aControlPoints
.mCP3
- aControlPoints
.mCP1
;
282 if (cp21
.x
== 0. && cp21
.y
== 0.) {
283 // In this case s3 becomes lim[n->0] (cp41.x * n) / n - (cp41.y * n) / n =
285 double s3
= cp41
.x
- cp41
.y
;
287 // Use the absolute value so that Min and Max will correspond with the
288 // minimum and maximum of the range.
293 double r
= CubicRoot(std::abs(aTolerance
/ s3
));
300 double s3
= (cp41
.x
* cp21
.y
- cp41
.y
* cp21
.x
) / hypot(cp21
.x
, cp21
.y
);
303 // This means within the precision we have it can be approximated
304 // infinitely by a linear segment. Deal with this by specifying the
305 // approximation range as extending beyond the entire curve.
311 double tf
= CubicRoot(std::abs(aTolerance
/ s3
));
313 *aMin
= aT
- tf
* (1 - aT
);
314 *aMax
= aT
+ tf
* (1 - aT
);
317 /* Find the inflection points of a bezier curve. Will return false if the
318 * curve is degenerate in such a way that it is best approximated by a straight
321 * The below algorithm was written by Jeff Muizelaar <jmuizelaar@mozilla.com>,
322 * explanation follows:
324 * The lower inflection point is returned in aT1, the higher one in aT2. In the
325 * case of a single inflection point this will be in aT1.
327 * The method is inspired by the algorithm in "analysis of in?ection points for
328 * planar cubic bezier curve"
330 * Here are some differences between this algorithm and versions discussed
331 * elsewhere in the literature:
333 * zhang et. al compute a0, d0 and e0 incrementally using the follow formula:
335 * Point a0 = CP2 - CP1
336 * Point a1 = CP3 - CP2
337 * Point a2 = CP4 - CP1
344 * this avoids any multiplications and may or may not be faster than the
345 * approach take below.
347 * "fast, precise flattening of cubic bezier path and ofset curves" by hain et.
349 * Point a = CP1 + 3 * CP2 - 3 * CP3 + CP4
350 * Point b = 3 * CP1 - 6 * CP2 + 3 * CP3
351 * Point c = -3 * CP1 + 3 * CP2
353 * the a, b, c, d can be expressed in terms of a0, d0 and e0 defined above as:
359 * a = 3a = a.y * b.x - a.x * b.y
360 * b = 3b = a.y * c.x - a.x * c.y
361 * c = 9c = b.y * c.x - b.x * c.y
363 * The additional multiples of 3 cancel each other out as show below:
365 * x = (-b + sqrt(b * b - 4 * a * c)) / (2 * a)
366 * x = (-3 * b + sqrt(3 * b * 3 * b - 4 * a * 3 * 9 * c / 3)) / (2 * 3 * a)
367 * x = 3 * (-b + sqrt(b * b - 4 * a * c)) / (2 * 3 * a)
368 * x = (-b + sqrt(b * b - 4 * a * c)) / (2 * a)
370 * I haven't looked into whether the formulation of the quadratic formula in
371 * hain has any numerical advantages over the one used below.
373 static inline void FindInflectionPoints(
374 const BezierControlPoints
& aControlPoints
, double* aT1
, double* aT2
,
376 // Find inflection points.
377 // See www.faculty.idc.ac.il/arik/quality/appendixa.html for an explanation
379 PointD A
= aControlPoints
.mCP2
- aControlPoints
.mCP1
;
381 aControlPoints
.mCP3
- (aControlPoints
.mCP2
* 2) + aControlPoints
.mCP1
;
382 PointD C
= aControlPoints
.mCP4
- (aControlPoints
.mCP3
* 3) +
383 (aControlPoints
.mCP2
* 3) - aControlPoints
.mCP1
;
385 double a
= B
.x
* C
.y
- B
.y
* C
.x
;
386 double b
= A
.x
* C
.y
- A
.y
* C
.x
;
387 double c
= A
.x
* B
.y
- A
.y
* B
.x
;
390 // Not a quadratic equation.
392 // Instead of a linear acceleration change we have a constant
393 // acceleration change. This means the equation has no solution
394 // and there are no inflection points, unless the constant is 0.
395 // In that case the curve is a straight line, essentially that means
396 // the easiest way to deal with is is by saying there's an inflection
397 // point at t == 0. The inflection point approximation range found will
398 // automatically extend into infinity.
412 double discriminant
= b
* b
- 4 * a
* c
;
414 if (discriminant
< 0) {
415 // No inflection points.
417 } else if (discriminant
== 0) {
421 /* Use the following formula for computing the roots:
423 * q = -1/2 * (b + sign(b) * sqrt(b^2 - 4ac))
427 double q
= sqrt(discriminant
);
438 std::swap(*aT1
, *aT2
);
444 void FlattenBezier(const BezierControlPoints
& aControlPoints
, PathSink
* aSink
,
450 FindInflectionPoints(aControlPoints
, &t1
, &t2
, &count
);
452 // Check that at least one of the inflection points is inside [0..1]
454 ((t1
< 0.0 || t1
>= 1.0) && (count
== 1 || (t2
< 0.0 || t2
>= 1.0)))) {
455 FlattenBezierCurveSegment(aControlPoints
, aSink
, aTolerance
);
459 double t1min
= t1
, t1max
= t1
, t2min
= t2
, t2max
= t2
;
461 BezierControlPoints remainingCP
= aControlPoints
;
463 // For both inflection points, calulate the range where they can be linearly
464 // approximated if they are positioned within [0,1]
465 if (count
> 0 && t1
>= 0 && t1
< 1.0) {
466 FindInflectionApproximationRange(aControlPoints
, &t1min
, &t1max
, t1
,
469 if (count
> 1 && t2
>= 0 && t2
< 1.0) {
470 FindInflectionApproximationRange(aControlPoints
, &t2min
, &t2max
, t2
,
473 BezierControlPoints nextCPs
= aControlPoints
;
474 BezierControlPoints prevCPs
;
476 // Process ranges. [t1min, t1max] and [t2min, t2max] are approximated by line
478 if (count
== 1 && t1min
<= 0 && t1max
>= 1.0) {
479 // The whole range can be approximated by a line segment.
480 aSink
->LineTo(aControlPoints
.mCP4
.ToPoint());
485 // Flatten the Bezier up until the first inflection point's approximation
487 SplitBezier(aControlPoints
, &prevCPs
, &remainingCP
, t1min
);
488 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
490 if (t1max
>= 0 && t1max
< 1.0 && (count
== 1 || t2min
> t1max
)) {
491 // The second inflection point's approximation range begins after the end
492 // of the first, approximate the first inflection point by a line and
493 // subsequently flatten up until the end or the next inflection point.
494 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
496 aSink
->LineTo(nextCPs
.mCP1
.ToPoint());
498 if (count
== 1 || (count
> 1 && t2min
>= 1.0)) {
499 // No more inflection points to deal with, flatten the rest of the curve.
500 FlattenBezierCurveSegment(nextCPs
, aSink
, aTolerance
);
502 } else if (count
> 1 && t2min
> 1.0) {
503 // We've already concluded t2min <= t1max, so if this is true the
504 // approximation range for the first inflection point runs past the
505 // end of the curve, draw a line to the end and we're done.
506 aSink
->LineTo(aControlPoints
.mCP4
.ToPoint());
510 if (count
> 1 && t2min
< 1.0 && t2max
> 0) {
511 if (t2min
> 0 && t2min
< t1max
) {
512 // In this case the t2 approximation range starts inside the t1
513 // approximation range.
514 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
515 aSink
->LineTo(nextCPs
.mCP1
.ToPoint());
516 } else if (t2min
> 0 && t1max
> 0) {
517 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
519 // Find a control points describing the portion of the curve between t1max
521 double t2mina
= (t2min
- t1max
) / (1 - t1max
);
522 SplitBezier(nextCPs
, &prevCPs
, &nextCPs
, t2mina
);
523 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
524 } else if (t2min
> 0) {
525 // We have nothing interesting before t2min, find that bit and flatten it.
526 SplitBezier(aControlPoints
, &prevCPs
, &nextCPs
, t2min
);
527 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
530 // Flatten the portion of the curve after t2max
531 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t2max
);
533 // Draw a line to the start, this is the approximation between t2min and
535 aSink
->LineTo(nextCPs
.mCP1
.ToPoint());
536 FlattenBezierCurveSegment(nextCPs
, aSink
, aTolerance
);
538 // Our approximation range extends beyond the end of the curve.
539 aSink
->LineTo(aControlPoints
.mCP4
.ToPoint());
545 Rect
Path::GetFastBounds(const Matrix
& aTransform
,
546 const StrokeOptions
* aStrokeOptions
) const {
547 return aStrokeOptions
? GetStrokedBounds(*aStrokeOptions
, aTransform
)
548 : GetBounds(aTransform
);
552 } // namespace mozilla