1 /* -*- Mode: C++; tab-width: 20; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 * This Source Code Form is subject to the terms of the Mozilla Public
3 * License, v. 2.0. If a copy of the MPL was not distributed with this
4 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
7 #include "PathAnalysis.h"
8 #include "PathHelpers.h"
13 static float CubicRoot(float aValue
) {
15 return -CubicRoot(-aValue
);
18 return powf(aValue
, 1.0f
/ 3.0f
);
22 struct BezierControlPoints
24 BezierControlPoints() {}
25 BezierControlPoints(const Point
&aCP1
, const Point
&aCP2
,
26 const Point
&aCP3
, const Point
&aCP4
)
27 : mCP1(aCP1
), mCP2(aCP2
), mCP3(aCP3
), mCP4(aCP4
)
31 Point mCP1
, mCP2
, mCP3
, mCP4
;
35 FlattenBezier(const BezierControlPoints
&aPoints
,
36 PathSink
*aSink
, Float aTolerance
);
50 EnsureFlattenedPath();
51 return mFlattenedPath
->ComputeLength();
55 Path::ComputePointAtLength(Float aLength
, Point
* aTangent
)
57 EnsureFlattenedPath();
58 return mFlattenedPath
->ComputePointAtLength(aLength
, aTangent
);
62 Path::EnsureFlattenedPath()
64 if (!mFlattenedPath
) {
65 mFlattenedPath
= new FlattenedPath();
66 StreamToSink(mFlattenedPath
);
70 // This is the maximum deviation we allow (with an additional ~20% margin of
71 // error) of the approximation from the actual Bezier curve.
72 const Float kFlatteningTolerance
= 0.0001f
;
75 FlattenedPath::MoveTo(const Point
&aPoint
)
77 MOZ_ASSERT(!mCalculatedLength
);
79 op
.mType
= FlatPathOp::OP_MOVETO
;
81 mPathOps
.push_back(op
);
87 FlattenedPath::LineTo(const Point
&aPoint
)
89 MOZ_ASSERT(!mCalculatedLength
);
91 op
.mType
= FlatPathOp::OP_LINETO
;
93 mPathOps
.push_back(op
);
97 FlattenedPath::BezierTo(const Point
&aCP1
,
101 MOZ_ASSERT(!mCalculatedLength
);
102 FlattenBezier(BezierControlPoints(CurrentPoint(), aCP1
, aCP2
, aCP3
), this, kFlatteningTolerance
);
106 FlattenedPath::QuadraticBezierTo(const Point
&aCP1
,
109 MOZ_ASSERT(!mCalculatedLength
);
110 // We need to elevate the degree of this quadratic B�zier to cubic, so we're
111 // going to add an intermediate control point, and recompute control point 1.
112 // The first and last control points remain the same.
113 // This formula can be found on http://fontforge.sourceforge.net/bezier.html
114 Point CP0
= CurrentPoint();
115 Point CP1
= (CP0
+ aCP1
* 2.0) / 3.0;
116 Point CP2
= (aCP2
+ aCP1
* 2.0) / 3.0;
119 BezierTo(CP1
, CP2
, CP3
);
123 FlattenedPath::Close()
125 MOZ_ASSERT(!mCalculatedLength
);
130 FlattenedPath::Arc(const Point
&aOrigin
, float aRadius
, float aStartAngle
,
131 float aEndAngle
, bool aAntiClockwise
)
133 ArcToBezier(this, aOrigin
, Size(aRadius
, aRadius
), aStartAngle
, aEndAngle
, aAntiClockwise
);
137 FlattenedPath::ComputeLength()
139 if (!mCalculatedLength
) {
142 for (uint32_t i
= 0; i
< mPathOps
.size(); i
++) {
143 if (mPathOps
[i
].mType
== FlatPathOp::OP_MOVETO
) {
144 currentPoint
= mPathOps
[i
].mPoint
;
146 mCachedLength
+= Distance(currentPoint
, mPathOps
[i
].mPoint
);
147 currentPoint
= mPathOps
[i
].mPoint
;
151 mCalculatedLength
= true;
154 return mCachedLength
;
158 FlattenedPath::ComputePointAtLength(Float aLength
, Point
*aTangent
)
160 // We track the last point that -wasn't- in the same place as the current
161 // point so if we pass the edge of the path with a bunch of zero length
162 // paths we still get the correct tangent vector.
163 Point lastPointSinceMove
;
165 for (uint32_t i
= 0; i
< mPathOps
.size(); i
++) {
166 if (mPathOps
[i
].mType
== FlatPathOp::OP_MOVETO
) {
167 if (Distance(currentPoint
, mPathOps
[i
].mPoint
)) {
168 lastPointSinceMove
= currentPoint
;
170 currentPoint
= mPathOps
[i
].mPoint
;
172 Float segmentLength
= Distance(currentPoint
, mPathOps
[i
].mPoint
);
175 lastPointSinceMove
= currentPoint
;
176 if (segmentLength
> aLength
) {
177 Point currentVector
= mPathOps
[i
].mPoint
- currentPoint
;
178 Point tangent
= currentVector
/ segmentLength
;
182 return currentPoint
+ tangent
* aLength
;
186 aLength
-= segmentLength
;
187 currentPoint
= mPathOps
[i
].mPoint
;
191 Point currentVector
= currentPoint
- lastPointSinceMove
;
193 if (hypotf(currentVector
.x
, currentVector
.y
)) {
194 *aTangent
= currentVector
/ hypotf(currentVector
.x
, currentVector
.y
);
202 // This function explicitly permits aControlPoints to refer to the same object
203 // as either of the other arguments.
205 SplitBezier(const BezierControlPoints
&aControlPoints
,
206 BezierControlPoints
*aFirstSegmentControlPoints
,
207 BezierControlPoints
*aSecondSegmentControlPoints
,
210 MOZ_ASSERT(aSecondSegmentControlPoints
);
212 *aSecondSegmentControlPoints
= aControlPoints
;
214 Point cp1a
= aControlPoints
.mCP1
+ (aControlPoints
.mCP2
- aControlPoints
.mCP1
) * t
;
215 Point cp2a
= aControlPoints
.mCP2
+ (aControlPoints
.mCP3
- aControlPoints
.mCP2
) * t
;
216 Point cp1aa
= cp1a
+ (cp2a
- cp1a
) * t
;
217 Point cp3a
= aControlPoints
.mCP3
+ (aControlPoints
.mCP4
- aControlPoints
.mCP3
) * t
;
218 Point cp2aa
= cp2a
+ (cp3a
- cp2a
) * t
;
219 Point cp1aaa
= cp1aa
+ (cp2aa
- cp1aa
) * t
;
220 aSecondSegmentControlPoints
->mCP4
= aControlPoints
.mCP4
;
222 if(aFirstSegmentControlPoints
) {
223 aFirstSegmentControlPoints
->mCP1
= aControlPoints
.mCP1
;
224 aFirstSegmentControlPoints
->mCP2
= cp1a
;
225 aFirstSegmentControlPoints
->mCP3
= cp1aa
;
226 aFirstSegmentControlPoints
->mCP4
= cp1aaa
;
228 aSecondSegmentControlPoints
->mCP1
= cp1aaa
;
229 aSecondSegmentControlPoints
->mCP2
= cp2aa
;
230 aSecondSegmentControlPoints
->mCP3
= cp3a
;
234 FlattenBezierCurveSegment(const BezierControlPoints
&aControlPoints
,
238 /* The algorithm implemented here is based on:
239 * http://cis.usouthal.edu/~hain/general/Publications/Bezier/Bezier%20Offset%20Curves.pdf
241 * The basic premise is that for a small t the third order term in the
242 * equation of a cubic bezier curve is insignificantly small. This can
243 * then be approximated by a quadratic equation for which the maximum
244 * difference from a linear approximation can be much more easily determined.
246 BezierControlPoints currentCP
= aControlPoints
;
250 Point cp21
= currentCP
.mCP2
- currentCP
.mCP3
;
251 Point cp31
= currentCP
.mCP3
- currentCP
.mCP1
;
253 Float s3
= (cp31
.x
* cp21
.y
- cp31
.y
* cp21
.x
) / hypotf(cp21
.x
, cp21
.y
);
255 t
= 2 * Float(sqrt(aTolerance
/ (3. * abs(s3
))));
258 aSink
->LineTo(aControlPoints
.mCP4
);
262 Point prevCP2
, prevCP3
, nextCP1
, nextCP2
, nextCP3
;
263 SplitBezier(currentCP
, nullptr, ¤tCP
, t
);
265 aSink
->LineTo(currentCP
.mCP1
);
270 FindInflectionApproximationRange(BezierControlPoints aControlPoints
,
271 Float
*aMin
, Float
*aMax
, Float aT
,
274 SplitBezier(aControlPoints
, nullptr, &aControlPoints
, aT
);
276 Point cp21
= aControlPoints
.mCP2
- aControlPoints
.mCP1
;
277 Point cp41
= aControlPoints
.mCP4
- aControlPoints
.mCP1
;
279 if (cp21
.x
== 0.f
&& cp21
.y
== 0.f
) {
280 // In this case s3 becomes lim[n->0] (cp41.x * n) / n - (cp41.y * n) / n = cp41.x - cp41.y.
282 // Use the absolute value so that Min and Max will correspond with the
283 // minimum and maximum of the range.
284 *aMin
= aT
- CubicRoot(abs(aTolerance
/ (cp41
.x
- cp41
.y
)));
285 *aMax
= aT
+ CubicRoot(abs(aTolerance
/ (cp41
.x
- cp41
.y
)));
289 Float s3
= (cp41
.x
* cp21
.y
- cp41
.y
* cp21
.x
) / hypotf(cp21
.x
, cp21
.y
);
292 // This means within the precision we have it can be approximated
293 // infinitely by a linear segment. Deal with this by specifying the
294 // approximation range as extending beyond the entire curve.
300 Float tf
= CubicRoot(abs(aTolerance
/ s3
));
302 *aMin
= aT
- tf
* (1 - aT
);
303 *aMax
= aT
+ tf
* (1 - aT
);
306 /* Find the inflection points of a bezier curve. Will return false if the
307 * curve is degenerate in such a way that it is best approximated by a straight
310 * The below algorithm was written by Jeff Muizelaar <jmuizelaar@mozilla.com>, explanation follows:
312 * The lower inflection point is returned in aT1, the higher one in aT2. In the
313 * case of a single inflection point this will be in aT1.
315 * The method is inspired by the algorithm in "analysis of in?ection points for planar cubic bezier curve"
317 * Here are some differences between this algorithm and versions discussed elsewhere in the literature:
319 * zhang et. al compute a0, d0 and e0 incrementally using the follow formula:
321 * Point a0 = CP2 - CP1
322 * Point a1 = CP3 - CP2
323 * Point a2 = CP4 - CP1
330 * this avoids any multiplications and may or may not be faster than the approach take below.
332 * "fast, precise flattening of cubic bezier path and ofset curves" by hain et. al
333 * Point a = CP1 + 3 * CP2 - 3 * CP3 + CP4
334 * Point b = 3 * CP1 - 6 * CP2 + 3 * CP3
335 * Point c = -3 * CP1 + 3 * CP2
337 * the a, b, c, d can be expressed in terms of a0, d0 and e0 defined above as:
343 * a = 3a = a.y * b.x - a.x * b.y
344 * b = 3b = a.y * c.x - a.x * c.y
345 * c = 9c = b.y * c.x - b.x * c.y
347 * The additional multiples of 3 cancel each other out as show below:
349 * x = (-b + sqrt(b * b - 4 * a * c)) / (2 * a)
350 * x = (-3 * b + sqrt(3 * b * 3 * b - 4 * a * 3 * 9 * c / 3)) / (2 * 3 * a)
351 * x = 3 * (-b + sqrt(b * b - 4 * a * c)) / (2 * 3 * a)
352 * x = (-b + sqrt(b * b - 4 * a * c)) / (2 * a)
354 * I haven't looked into whether the formulation of the quadratic formula in
355 * hain has any numerical advantages over the one used below.
358 FindInflectionPoints(const BezierControlPoints
&aControlPoints
,
359 Float
*aT1
, Float
*aT2
, uint32_t *aCount
)
361 // Find inflection points.
362 // See www.faculty.idc.ac.il/arik/quality/appendixa.html for an explanation
364 Point A
= aControlPoints
.mCP2
- aControlPoints
.mCP1
;
365 Point B
= aControlPoints
.mCP3
- (aControlPoints
.mCP2
* 2) + aControlPoints
.mCP1
;
366 Point C
= aControlPoints
.mCP4
- (aControlPoints
.mCP3
* 3) + (aControlPoints
.mCP2
* 3) - aControlPoints
.mCP1
;
368 Float a
= Float(B
.x
) * C
.y
- Float(B
.y
) * C
.x
;
369 Float b
= Float(A
.x
) * C
.y
- Float(A
.y
) * C
.x
;
370 Float c
= Float(A
.x
) * B
.y
- Float(A
.y
) * B
.x
;
373 // Not a quadratic equation.
375 // Instead of a linear acceleration change we have a constant
376 // acceleration change. This means the equation has no solution
377 // and there are no inflection points, unless the constant is 0.
378 // In that case the curve is a straight line, essentially that means
379 // the easiest way to deal with is is by saying there's an inflection
380 // point at t == 0. The inflection point approximation range found will
381 // automatically extend into infinity.
394 Float discriminant
= b
* b
- 4 * a
* c
;
396 if (discriminant
< 0) {
397 // No inflection points.
399 } else if (discriminant
== 0) {
403 /* Use the following formula for computing the roots:
405 * q = -1/2 * (b + sign(b) * sqrt(b^2 - 4ac))
409 Float q
= sqrtf(discriminant
);
420 std::swap(*aT1
, *aT2
);
430 FlattenBezier(const BezierControlPoints
&aControlPoints
,
431 PathSink
*aSink
, Float aTolerance
)
437 FindInflectionPoints(aControlPoints
, &t1
, &t2
, &count
);
439 // Check that at least one of the inflection points is inside [0..1]
440 if (count
== 0 || ((t1
< 0 || t1
> 1.0) && ((t2
< 0 || t2
> 1.0) || count
== 1)) ) {
441 FlattenBezierCurveSegment(aControlPoints
, aSink
, aTolerance
);
445 Float t1min
= t1
, t1max
= t1
, t2min
= t2
, t2max
= t2
;
447 BezierControlPoints remainingCP
= aControlPoints
;
449 // For both inflection points, calulate the range where they can be linearly
450 // approximated if they are positioned within [0,1]
451 if (count
> 0 && t1
>= 0 && t1
< 1.0) {
452 FindInflectionApproximationRange(aControlPoints
, &t1min
, &t1max
, t1
, aTolerance
);
454 if (count
> 1 && t2
>= 0 && t2
< 1.0) {
455 FindInflectionApproximationRange(aControlPoints
, &t2min
, &t2max
, t2
, aTolerance
);
457 BezierControlPoints nextCPs
= aControlPoints
;
458 BezierControlPoints prevCPs
;
460 // Process ranges. [t1min, t1max] and [t2min, t2max] are approximated by line
463 // Flatten the Bezier up until the first inflection point's approximation
465 SplitBezier(aControlPoints
, &prevCPs
,
466 &remainingCP
, t1min
);
467 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
469 if (t1max
>= 0 && t1max
< 1.0 && (count
== 1 || t2min
> t1max
)) {
470 // The second inflection point's approximation range begins after the end
471 // of the first, approximate the first inflection point by a line and
472 // subsequently flatten up until the end or the next inflection point.
473 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
475 aSink
->LineTo(nextCPs
.mCP1
);
477 if (count
== 1 || (count
> 1 && t2min
>= 1.0)) {
478 // No more inflection points to deal with, flatten the rest of the curve.
479 FlattenBezierCurveSegment(nextCPs
, aSink
, aTolerance
);
481 } else if (count
> 1 && t2min
> 1.0) {
482 // We've already concluded t2min <= t1max, so if this is true the
483 // approximation range for the first inflection point runs past the
484 // end of the curve, draw a line to the end and we're done.
485 aSink
->LineTo(aControlPoints
.mCP4
);
489 if (count
> 1 && t2min
< 1.0 && t2max
> 0) {
490 if (t2min
> 0 && t2min
< t1max
) {
491 // In this case the t2 approximation range starts inside the t1
492 // approximation range.
493 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
494 aSink
->LineTo(nextCPs
.mCP1
);
495 } else if (t2min
> 0 && t1max
> 0) {
496 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t1max
);
498 // Find a control points describing the portion of the curve between t1max and t2min.
499 Float t2mina
= (t2min
- t1max
) / (1 - t1max
);
500 SplitBezier(nextCPs
, &prevCPs
, &nextCPs
, t2mina
);
501 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
502 } else if (t2min
> 0) {
503 // We have nothing interesting before t2min, find that bit and flatten it.
504 SplitBezier(aControlPoints
, &prevCPs
, &nextCPs
, t2min
);
505 FlattenBezierCurveSegment(prevCPs
, aSink
, aTolerance
);
508 // Flatten the portion of the curve after t2max
509 SplitBezier(aControlPoints
, nullptr, &nextCPs
, t2max
);
511 // Draw a line to the start, this is the approximation between t2min and
513 aSink
->LineTo(nextCPs
.mCP1
);
514 FlattenBezierCurveSegment(nextCPs
, aSink
, aTolerance
);
516 // Our approximation range extends beyond the end of the curve.
517 aSink
->LineTo(aControlPoints
.mCP4
);