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 // We track the last point that -wasn't- in the same place as the current
142 // point so if we pass the edge of the path with a bunch of zero length
143 // paths we still get the correct tangent vector.
144 Point lastPointSinceMove
;
146 for (uint32_t i
= 0; i
< mPathOps
.size(); i
++) {
147 if (mPathOps
[i
].mType
== FlatPathOp::OP_MOVETO
) {
148 if (Distance(currentPoint
, mPathOps
[i
].mPoint
)) {
149 lastPointSinceMove
= currentPoint
;
151 currentPoint
= mPathOps
[i
].mPoint
;
153 Float segmentLength
= Distance(currentPoint
, mPathOps
[i
].mPoint
);
156 lastPointSinceMove
= currentPoint
;
157 if (segmentLength
> aLength
) {
158 Point currentVector
= mPathOps
[i
].mPoint
- currentPoint
;
159 Point tangent
= currentVector
/ segmentLength
;
163 return currentPoint
+ tangent
* aLength
;
167 aLength
-= segmentLength
;
168 currentPoint
= mPathOps
[i
].mPoint
;
172 Point currentVector
= currentPoint
- lastPointSinceMove
;
174 if (hypotf(currentVector
.x
, currentVector
.y
)) {
175 *aTangent
= currentVector
/ hypotf(currentVector
.x
, currentVector
.y
);
183 // This function explicitly permits aControlPoints to refer to the same object
184 // as either of the other arguments.
185 static void SplitBezier(const BezierControlPoints
& aControlPoints
,
186 BezierControlPoints
* aFirstSegmentControlPoints
,
187 BezierControlPoints
* aSecondSegmentControlPoints
,
189 MOZ_ASSERT(aSecondSegmentControlPoints
);
191 *aSecondSegmentControlPoints
= aControlPoints
;
194 aControlPoints
.mCP1
+ (aControlPoints
.mCP2
- aControlPoints
.mCP1
) * t
;
196 aControlPoints
.mCP2
+ (aControlPoints
.mCP3
- aControlPoints
.mCP2
) * t
;
197 PointD cp1aa
= cp1a
+ (cp2a
- cp1a
) * t
;
199 aControlPoints
.mCP3
+ (aControlPoints
.mCP4
- aControlPoints
.mCP3
) * t
;
200 PointD cp2aa
= cp2a
+ (cp3a
- cp2a
) * t
;
201 PointD cp1aaa
= cp1aa
+ (cp2aa
- cp1aa
) * t
;
202 aSecondSegmentControlPoints
->mCP4
= aControlPoints
.mCP4
;
204 if (aFirstSegmentControlPoints
) {
205 aFirstSegmentControlPoints
->mCP1
= aControlPoints
.mCP1
;
206 aFirstSegmentControlPoints
->mCP2
= cp1a
;
207 aFirstSegmentControlPoints
->mCP3
= cp1aa
;
208 aFirstSegmentControlPoints
->mCP4
= cp1aaa
;
210 aSecondSegmentControlPoints
->mCP1
= cp1aaa
;
211 aSecondSegmentControlPoints
->mCP2
= cp2aa
;
212 aSecondSegmentControlPoints
->mCP3
= cp3a
;
215 static void FlattenBezierCurveSegment(const BezierControlPoints
& aControlPoints
,
216 PathSink
* aSink
, double aTolerance
) {
217 /* The algorithm implemented here is based on:
218 * http://cis.usouthal.edu/~hain/general/Publications/Bezier/Bezier%20Offset%20Curves.pdf
220 * The basic premise is that for a small t the third order term in the
221 * equation of a cubic bezier curve is insignificantly small. This can
222 * then be approximated by a quadratic equation for which the maximum
223 * difference from a linear approximation can be much more easily determined.
225 BezierControlPoints currentCP
= aControlPoints
;
228 double currentTolerance
= aTolerance
;
230 PointD cp21
= currentCP
.mCP2
- currentCP
.mCP1
;
231 PointD cp31
= currentCP
.mCP3
- currentCP
.mCP1
;
233 /* To remove divisions and check for divide-by-zero, this is optimized from:
234 * Float s3 = (cp31.x * cp21.y - cp31.y * cp21.x) / hypotf(cp21.x, cp21.y);
235 * t = 2 * Float(sqrt(aTolerance / (3. * std::abs(s3))));
237 double cp21x31
= cp31
.x
* cp21
.y
- cp31
.y
* cp21
.x
;
238 double h
= hypot(cp21
.x
, cp21
.y
);
239 if (cp21x31
* h
== 0) {
243 double s3inv
= h
/ cp21x31
;
244 t
= 2 * sqrt(currentTolerance
* std::abs(s3inv
) / 3.);
245 currentTolerance
*= 1 + aTolerance
;
246 // Increase tolerance every iteration to prevent this loop from executing
247 // too many times. This approximates the length of large curves more
248 // roughly. In practice, aTolerance is the constant kFlatteningTolerance
249 // which has value 0.0001. With this value, it takes 6,932 splits to double
250 // currentTolerance (to 0.0002) and 23,028 splits to increase
251 // currentTolerance by an order of magnitude (to 0.001).
256 SplitBezier(currentCP
, nullptr, ¤tCP
, t
);
258 aSink
->LineTo(currentCP
.mCP1
.ToPoint());
261 aSink
->LineTo(currentCP
.mCP4
.ToPoint());
264 static inline void FindInflectionApproximationRange(
265 BezierControlPoints aControlPoints
, double* aMin
, double* aMax
, double aT
,
267 SplitBezier(aControlPoints
, nullptr, &aControlPoints
, aT
);
269 PointD cp21
= aControlPoints
.mCP2
- aControlPoints
.mCP1
;
270 PointD cp41
= aControlPoints
.mCP4
- aControlPoints
.mCP1
;
272 if (cp21
.x
== 0. && cp21
.y
== 0.) {
273 cp21
= aControlPoints
.mCP3
- aControlPoints
.mCP1
;
276 if (cp21
.x
== 0. && cp21
.y
== 0.) {
277 // In this case s3 becomes lim[n->0] (cp41.x * n) / n - (cp41.y * n) / n =
279 double s3
= cp41
.x
- cp41
.y
;
281 // Use the absolute value so that Min and Max will correspond with the
282 // minimum and maximum of the range.
287 double r
= CubicRoot(std::abs(aTolerance
/ s3
));
294 double s3
= (cp41
.x
* cp21
.y
- cp41
.y
* cp21
.x
) / hypot(cp21
.x
, cp21
.y
);
297 // This means within the precision we have it can be approximated
298 // infinitely by a linear segment. Deal with this by specifying the
299 // approximation range as extending beyond the entire curve.
305 double tf
= CubicRoot(std::abs(aTolerance
/ s3
));
307 *aMin
= aT
- tf
* (1 - aT
);
308 *aMax
= aT
+ tf
* (1 - aT
);
311 /* Find the inflection points of a bezier curve. Will return false if the
312 * curve is degenerate in such a way that it is best approximated by a straight
315 * The below algorithm was written by Jeff Muizelaar <jmuizelaar@mozilla.com>,
316 * explanation follows:
318 * The lower inflection point is returned in aT1, the higher one in aT2. In the
319 * case of a single inflection point this will be in aT1.
321 * The method is inspired by the algorithm in "analysis of in?ection points for
322 * planar cubic bezier curve"
324 * Here are some differences between this algorithm and versions discussed
325 * elsewhere in the literature:
327 * zhang et. al compute a0, d0 and e0 incrementally using the follow formula:
329 * Point a0 = CP2 - CP1
330 * Point a1 = CP3 - CP2
331 * Point a2 = CP4 - CP1
338 * this avoids any multiplications and may or may not be faster than the
339 * approach take below.
341 * "fast, precise flattening of cubic bezier path and ofset curves" by hain et.
343 * Point a = CP1 + 3 * CP2 - 3 * CP3 + CP4
344 * Point b = 3 * CP1 - 6 * CP2 + 3 * CP3
345 * Point c = -3 * CP1 + 3 * CP2
347 * the a, b, c, d can be expressed in terms of a0, d0 and e0 defined above as:
353 * a = 3a = a.y * b.x - a.x * b.y
354 * b = 3b = a.y * c.x - a.x * c.y
355 * c = 9c = b.y * c.x - b.x * c.y
357 * The additional multiples of 3 cancel each other out as show below:
359 * x = (-b + sqrt(b * b - 4 * a * c)) / (2 * a)
360 * x = (-3 * b + sqrt(3 * b * 3 * b - 4 * a * 3 * 9 * c / 3)) / (2 * 3 * a)
361 * x = 3 * (-b + sqrt(b * b - 4 * a * c)) / (2 * 3 * a)
362 * x = (-b + sqrt(b * b - 4 * a * c)) / (2 * a)
364 * I haven't looked into whether the formulation of the quadratic formula in
365 * hain has any numerical advantages over the one used below.
367 static inline void FindInflectionPoints(
368 const BezierControlPoints
& aControlPoints
, double* aT1
, double* aT2
,
370 // Find inflection points.
371 // See www.faculty.idc.ac.il/arik/quality/appendixa.html for an explanation
373 PointD A
= aControlPoints
.mCP2
- aControlPoints
.mCP1
;
375 aControlPoints
.mCP3
- (aControlPoints
.mCP2
* 2) + aControlPoints
.mCP1
;
376 PointD C
= aControlPoints
.mCP4
- (aControlPoints
.mCP3
* 3) +
377 (aControlPoints
.mCP2
* 3) - aControlPoints
.mCP1
;
379 double a
= B
.x
* C
.y
- B
.y
* C
.x
;
380 double b
= A
.x
* C
.y
- A
.y
* C
.x
;
381 double c
= A
.x
* B
.y
- A
.y
* B
.x
;
384 // Not a quadratic equation.
386 // Instead of a linear acceleration change we have a constant
387 // acceleration change. This means the equation has no solution
388 // and there are no inflection points, unless the constant is 0.
389 // In that case the curve is a straight line, essentially that means
390 // the easiest way to deal with is is by saying there's an inflection
391 // point at t == 0. The inflection point approximation range found will
392 // automatically extend into infinity.
405 double discriminant
= b
* b
- 4 * a
* c
;
407 if (discriminant
< 0) {
408 // No inflection points.
410 } else if (discriminant
== 0) {
414 /* Use the following formula for computing the roots:
416 * q = -1/2 * (b + sign(b) * sqrt(b^2 - 4ac))
420 double q
= sqrt(discriminant
);
431 std::swap(*aT1
, *aT2
);
438 void FlattenBezier(const BezierControlPoints
& aControlPoints
, PathSink
* aSink
,
444 FindInflectionPoints(aControlPoints
, &t1
, &t2
, &count
);
446 // Check that at least one of the inflection points is inside [0..1]
448 ((t1
< 0.0 || t1
>= 1.0) && (count
== 1 || (t2
< 0.0 || t2
>= 1.0)))) {
449 FlattenBezierCurveSegment(aControlPoints
, aSink
, aTolerance
);
453 double t1min
= t1
, t1max
= t1
, t2min
= t2
, t2max
= t2
;
455 BezierControlPoints remainingCP
= aControlPoints
;
457 // For both inflection points, calulate the range where they can be linearly
458 // approximated if they are positioned within [0,1]
459 if (count
> 0 && t1
>= 0 && t1
< 1.0) {
460 FindInflectionApproximationRange(aControlPoints
, &t1min
, &t1max
, t1
,
463 if (count
> 1 && t2
>= 0 && t2
< 1.0) {
464 FindInflectionApproximationRange(aControlPoints
, &t2min
, &t2max
, t2
,
467 BezierControlPoints nextCPs
= aControlPoints
;
468 BezierControlPoints prevCPs
;
470 // Process ranges. [t1min, t1max] and [t2min, t2max] are approximated by line
472 if (count
== 1 && t1min
<= 0 && t1max
>= 1.0) {
473 // The whole range can be approximated by a line segment.
474 aSink
->LineTo(aControlPoints
.mCP4
.ToPoint());
479 // Flatten the Bezier up until the first inflection point's approximation
481 SplitBezier(aControlPoints
, &prevCPs
, &remainingCP
, t1min
);
482 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
484 if (t1max
>= 0 && t1max
< 1.0 && (count
== 1 || t2min
> t1max
)) {
485 // The second inflection point's approximation range begins after the end
486 // of the first, approximate the first inflection point by a line and
487 // subsequently flatten up until the end or the next inflection point.
488 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
490 aSink
->LineTo(nextCPs
.mCP1
.ToPoint());
492 if (count
== 1 || (count
> 1 && t2min
>= 1.0)) {
493 // No more inflection points to deal with, flatten the rest of the curve.
494 FlattenBezierCurveSegment(nextCPs
, aSink
, aTolerance
);
496 } else if (count
> 1 && t2min
> 1.0) {
497 // We've already concluded t2min <= t1max, so if this is true the
498 // approximation range for the first inflection point runs past the
499 // end of the curve, draw a line to the end and we're done.
500 aSink
->LineTo(aControlPoints
.mCP4
.ToPoint());
504 if (count
> 1 && t2min
< 1.0 && t2max
> 0) {
505 if (t2min
> 0 && t2min
< t1max
) {
506 // In this case the t2 approximation range starts inside the t1
507 // approximation range.
508 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
509 aSink
->LineTo(nextCPs
.mCP1
.ToPoint());
510 } else if (t2min
> 0 && t1max
> 0) {
511 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
513 // Find a control points describing the portion of the curve between t1max
515 double t2mina
= (t2min
- t1max
) / (1 - t1max
);
516 SplitBezier(nextCPs
, &prevCPs
, &nextCPs
, t2mina
);
517 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
518 } else if (t2min
> 0) {
519 // We have nothing interesting before t2min, find that bit and flatten it.
520 SplitBezier(aControlPoints
, &prevCPs
, &nextCPs
, t2min
);
521 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
524 // Flatten the portion of the curve after t2max
525 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t2max
);
527 // Draw a line to the start, this is the approximation between t2min and
529 aSink
->LineTo(nextCPs
.mCP1
.ToPoint());
530 FlattenBezierCurveSegment(nextCPs
, aSink
, aTolerance
);
532 // Our approximation range extends beyond the end of the curve.
533 aSink
->LineTo(aControlPoints
.mCP4
.ToPoint());
540 } // namespace mozilla