1 /*
2  * Copyright 2011 Google Inc.
3  *
4  * Use of this source code is governed by a BSD-style license that can be
5  * found in the LICENSE file.
6  */
7 
8 #include "GrPathUtils.h"
9 
10 #include "GrTypes.h"
11 #include "SkMathPriv.h"
12 #include "SkPointPriv.h"
13 
14 static const SkScalar gMinCurveTol = 0.0001f;
15 
16 SkScalar GrPathUtils::scaleToleranceToSrc(SkScalar devTol,
17                                           const SkMatrix& viewM,
18                                           const SkRect& pathBounds) {
19     // In order to tesselate the path we get a bound on how much the matrix can
20     // scale when mapping to screen coordinates.
21     SkScalar stretch = viewM.getMaxScale();
22 
23     if (stretch < 0) {
24         // take worst case mapRadius amoung four corners.
25         // (less than perfect)
26         for (int i = 0; i < 4; ++i) {
27             SkMatrix mat;
28             mat.setTranslate((i % 2) ? pathBounds.fLeft : pathBounds.fRight,
29                              (i < 2) ? pathBounds.fTop : pathBounds.fBottom);
30             mat.postConcat(viewM);
31             stretch = SkMaxScalar(stretch, mat.mapRadius(SK_Scalar1));
32         }
33     }
34     SkScalar srcTol = 0;
35     if (stretch <= 0) {
36         // We have degenerate bounds or some degenerate matrix. Thus we set the tolerance to be the
37         // max of the path pathBounds width and height.
38         srcTol = SkTMax(pathBounds.width(), pathBounds.height());
39     } else {
40         srcTol = devTol / stretch;
41     }
42     if (srcTol < gMinCurveTol) {
43         srcTol = gMinCurveTol;
44     }
45     return srcTol;
46 }
47 
48 uint32_t GrPathUtils::quadraticPointCount(const SkPoint points[], SkScalar tol) {
49     // You should have called scaleToleranceToSrc, which guarantees this
50     SkASSERT(tol >= gMinCurveTol);
51 
52     SkScalar d = SkPointPriv::DistanceToLineSegmentBetween(points[1], points[0], points[2]);
53     if (!SkScalarIsFinite(d)) {
54         return kMaxPointsPerCurve;
55     } else if (d <= tol) {
56         return 1;
57     } else {
58         // Each time we subdivide, d should be cut in 4. So we need to
59         // subdivide x = log4(d/tol) times. x subdivisions creates 2^(x)
60         // points.
61         // 2^(log4(x)) = sqrt(x);
62         SkScalar divSqrt = SkScalarSqrt(d / tol);
63         if (((SkScalar)SK_MaxS32) <= divSqrt) {
64             return kMaxPointsPerCurve;
65         } else {
66             int temp = SkScalarCeilToInt(divSqrt);
67             int pow2 = GrNextPow2(temp);
68             // Because of NaNs & INFs we can wind up with a degenerate temp
69             // such that pow2 comes out negative. Also, our point generator
70             // will always output at least one pt.
71             if (pow2 < 1) {
72                 pow2 = 1;
73             }
74             return SkTMin(pow2, kMaxPointsPerCurve);
75         }
76     }
77 }
78 
79 uint32_t GrPathUtils::generateQuadraticPoints(const SkPoint& p0,
80                                               const SkPoint& p1,
81                                               const SkPoint& p2,
82                                               SkScalar tolSqd,
83                                               SkPoint** points,
84                                               uint32_t pointsLeft) {
85     if (pointsLeft < 2 ||
86         (SkPointPriv::DistanceToLineSegmentBetweenSqd(p1, p0, p2)) < tolSqd) {
87         (*points)[0] = p2;
88         *points += 1;
89         return 1;
90     }
91 
92     SkPoint q[] = {
93         { SkScalarAve(p0.fX, p1.fX), SkScalarAve(p0.fY, p1.fY) },
94         { SkScalarAve(p1.fX, p2.fX), SkScalarAve(p1.fY, p2.fY) },
95     };
96     SkPoint r = { SkScalarAve(q[0].fX, q[1].fX), SkScalarAve(q[0].fY, q[1].fY) };
97 
98     pointsLeft >>= 1;
99     uint32_t a = generateQuadraticPoints(p0, q[0], r, tolSqd, points, pointsLeft);
100     uint32_t b = generateQuadraticPoints(r, q[1], p2, tolSqd, points, pointsLeft);
101     return a + b;
102 }
103 
104 uint32_t GrPathUtils::cubicPointCount(const SkPoint points[],
105                                            SkScalar tol) {
106     // You should have called scaleToleranceToSrc, which guarantees this
107     SkASSERT(tol >= gMinCurveTol);
108 
109     SkScalar d = SkTMax(
110         SkPointPriv::DistanceToLineSegmentBetweenSqd(points[1], points[0], points[3]),
111         SkPointPriv::DistanceToLineSegmentBetweenSqd(points[2], points[0], points[3]));
112     d = SkScalarSqrt(d);
113     if (!SkScalarIsFinite(d)) {
114         return kMaxPointsPerCurve;
115     } else if (d <= tol) {
116         return 1;
117     } else {
118         SkScalar divSqrt = SkScalarSqrt(d / tol);
119         if (((SkScalar)SK_MaxS32) <= divSqrt) {
120             return kMaxPointsPerCurve;
121         } else {
122             int temp = SkScalarCeilToInt(SkScalarSqrt(d / tol));
123             int pow2 = GrNextPow2(temp);
124             // Because of NaNs & INFs we can wind up with a degenerate temp
125             // such that pow2 comes out negative. Also, our point generator
126             // will always output at least one pt.
127             if (pow2 < 1) {
128                 pow2 = 1;
129             }
130             return SkTMin(pow2, kMaxPointsPerCurve);
131         }
132     }
133 }
134 
135 uint32_t GrPathUtils::generateCubicPoints(const SkPoint& p0,
136                                           const SkPoint& p1,
137                                           const SkPoint& p2,
138                                           const SkPoint& p3,
139                                           SkScalar tolSqd,
140                                           SkPoint** points,
141                                           uint32_t pointsLeft) {
142     if (pointsLeft < 2 ||
143         (SkPointPriv::DistanceToLineSegmentBetweenSqd(p1, p0, p3) < tolSqd &&
144          SkPointPriv::DistanceToLineSegmentBetweenSqd(p2, p0, p3) < tolSqd)) {
145         (*points)[0] = p3;
146         *points += 1;
147         return 1;
148     }
149     SkPoint q[] = {
150         { SkScalarAve(p0.fX, p1.fX), SkScalarAve(p0.fY, p1.fY) },
151         { SkScalarAve(p1.fX, p2.fX), SkScalarAve(p1.fY, p2.fY) },
152         { SkScalarAve(p2.fX, p3.fX), SkScalarAve(p2.fY, p3.fY) }
153     };
154     SkPoint r[] = {
155         { SkScalarAve(q[0].fX, q[1].fX), SkScalarAve(q[0].fY, q[1].fY) },
156         { SkScalarAve(q[1].fX, q[2].fX), SkScalarAve(q[1].fY, q[2].fY) }
157     };
158     SkPoint s = { SkScalarAve(r[0].fX, r[1].fX), SkScalarAve(r[0].fY, r[1].fY) };
159     pointsLeft >>= 1;
160     uint32_t a = generateCubicPoints(p0, q[0], r[0], s, tolSqd, points, pointsLeft);
161     uint32_t b = generateCubicPoints(s, r[1], q[2], p3, tolSqd, points, pointsLeft);
162     return a + b;
163 }
164 
165 int GrPathUtils::worstCasePointCount(const SkPath& path, int* subpaths, SkScalar tol) {
166     // You should have called scaleToleranceToSrc, which guarantees this
167     SkASSERT(tol >= gMinCurveTol);
168 
169     int pointCount = 0;
170     *subpaths = 1;
171 
172     bool first = true;
173 
174     SkPath::Iter iter(path, false);
175     SkPath::Verb verb;
176 
177     SkPoint pts[4];
178     while ((verb = iter.next(pts, false)) != SkPath::kDone_Verb) {
179 
180         switch (verb) {
181             case SkPath::kLine_Verb:
182                 pointCount += 1;
183                 break;
184             case SkPath::kConic_Verb: {
185                 SkScalar weight = iter.conicWeight();
186                 SkAutoConicToQuads converter;
187                 const SkPoint* quadPts = converter.computeQuads(pts, weight, tol);
188                 for (int i = 0; i < converter.countQuads(); ++i) {
189                     pointCount += quadraticPointCount(quadPts + 2*i, tol);
190                 }
191             }
192             case SkPath::kQuad_Verb:
193                 pointCount += quadraticPointCount(pts, tol);
194                 break;
195             case SkPath::kCubic_Verb:
196                 pointCount += cubicPointCount(pts, tol);
197                 break;
198             case SkPath::kMove_Verb:
199                 pointCount += 1;
200                 if (!first) {
201                     ++(*subpaths);
202                 }
203                 break;
204             default:
205                 break;
206         }
207         first = false;
208     }
209     return pointCount;
210 }
211 
212 void GrPathUtils::QuadUVMatrix::set(const SkPoint qPts[3]) {
213     SkMatrix m;
214     // We want M such that M * xy_pt = uv_pt
215     // We know M * control_pts = [0  1/2 1]
216     //                           [0  0   1]
217     //                           [1  1   1]
218     // And control_pts = [x0 x1 x2]
219     //                   [y0 y1 y2]
220     //                   [1  1  1 ]
221     // We invert the control pt matrix and post concat to both sides to get M.
222     // Using the known form of the control point matrix and the result, we can
223     // optimize and improve precision.
224 
225     double x0 = qPts[0].fX;
226     double y0 = qPts[0].fY;
227     double x1 = qPts[1].fX;
228     double y1 = qPts[1].fY;
229     double x2 = qPts[2].fX;
230     double y2 = qPts[2].fY;
231     double det = x0*y1 - y0*x1 + x2*y0 - y2*x0 + x1*y2 - y1*x2;
232 
233     if (!sk_float_isfinite(det)
234         || SkScalarNearlyZero((float)det, SK_ScalarNearlyZero * SK_ScalarNearlyZero)) {
235         // The quad is degenerate. Hopefully this is rare. Find the pts that are
236         // farthest apart to compute a line (unless it is really a pt).
237         SkScalar maxD = SkPointPriv::DistanceToSqd(qPts[0], qPts[1]);
238         int maxEdge = 0;
239         SkScalar d = SkPointPriv::DistanceToSqd(qPts[1], qPts[2]);
240         if (d > maxD) {
241             maxD = d;
242             maxEdge = 1;
243         }
244         d = SkPointPriv::DistanceToSqd(qPts[2], qPts[0]);
245         if (d > maxD) {
246             maxD = d;
247             maxEdge = 2;
248         }
249         // We could have a tolerance here, not sure if it would improve anything
250         if (maxD > 0) {
251             // Set the matrix to give (u = 0, v = distance_to_line)
252             SkVector lineVec = qPts[(maxEdge + 1)%3] - qPts[maxEdge];
253             // when looking from the point 0 down the line we want positive
254             // distances to be to the left. This matches the non-degenerate
255             // case.
256             lineVec = SkPointPriv::MakeOrthog(lineVec, SkPointPriv::kLeft_Side);
257             // first row
258             fM[0] = 0;
259             fM[1] = 0;
260             fM[2] = 0;
261             // second row
262             fM[3] = lineVec.fX;
263             fM[4] = lineVec.fY;
264             fM[5] = -lineVec.dot(qPts[maxEdge]);
265         } else {
266             // It's a point. It should cover zero area. Just set the matrix such
267             // that (u, v) will always be far away from the quad.
268             fM[0] = 0; fM[1] = 0; fM[2] = 100.f;
269             fM[3] = 0; fM[4] = 0; fM[5] = 100.f;
270         }
271     } else {
272         double scale = 1.0/det;
273 
274         // compute adjugate matrix
275         double a2, a3, a4, a5, a6, a7, a8;
276         a2 = x1*y2-x2*y1;
277 
278         a3 = y2-y0;
279         a4 = x0-x2;
280         a5 = x2*y0-x0*y2;
281 
282         a6 = y0-y1;
283         a7 = x1-x0;
284         a8 = x0*y1-x1*y0;
285 
286         // this performs the uv_pts*adjugate(control_pts) multiply,
287         // then does the scale by 1/det afterwards to improve precision
288         m[SkMatrix::kMScaleX] = (float)((0.5*a3 + a6)*scale);
289         m[SkMatrix::kMSkewX]  = (float)((0.5*a4 + a7)*scale);
290         m[SkMatrix::kMTransX] = (float)((0.5*a5 + a8)*scale);
291 
292         m[SkMatrix::kMSkewY]  = (float)(a6*scale);
293         m[SkMatrix::kMScaleY] = (float)(a7*scale);
294         m[SkMatrix::kMTransY] = (float)(a8*scale);
295 
296         // kMPersp0 & kMPersp1 should algebraically be zero
297         m[SkMatrix::kMPersp0] = 0.0f;
298         m[SkMatrix::kMPersp1] = 0.0f;
299         m[SkMatrix::kMPersp2] = (float)((a2 + a5 + a8)*scale);
300 
301         // It may not be normalized to have 1.0 in the bottom right
302         float m33 = m.get(SkMatrix::kMPersp2);
303         if (1.f != m33) {
304             m33 = 1.f / m33;
305             fM[0] = m33 * m.get(SkMatrix::kMScaleX);
306             fM[1] = m33 * m.get(SkMatrix::kMSkewX);
307             fM[2] = m33 * m.get(SkMatrix::kMTransX);
308             fM[3] = m33 * m.get(SkMatrix::kMSkewY);
309             fM[4] = m33 * m.get(SkMatrix::kMScaleY);
310             fM[5] = m33 * m.get(SkMatrix::kMTransY);
311         } else {
312             fM[0] = m.get(SkMatrix::kMScaleX);
313             fM[1] = m.get(SkMatrix::kMSkewX);
314             fM[2] = m.get(SkMatrix::kMTransX);
315             fM[3] = m.get(SkMatrix::kMSkewY);
316             fM[4] = m.get(SkMatrix::kMScaleY);
317             fM[5] = m.get(SkMatrix::kMTransY);
318         }
319     }
320 }
321 
322 ////////////////////////////////////////////////////////////////////////////////
323 
324 // k = (y2 - y0, x0 - x2, x2*y0 - x0*y2)
325 // l = (y1 - y0, x0 - x1, x1*y0 - x0*y1) * 2*w
326 // m = (y2 - y1, x1 - x2, x2*y1 - x1*y2) * 2*w
327 void GrPathUtils::getConicKLM(const SkPoint p[3], const SkScalar weight, SkMatrix* out) {
328     SkMatrix& klm = *out;
329     const SkScalar w2 = 2.f * weight;
330     klm[0] = p[2].fY - p[0].fY;
331     klm[1] = p[0].fX - p[2].fX;
332     klm[2] = p[2].fX * p[0].fY - p[0].fX * p[2].fY;
333 
334     klm[3] = w2 * (p[1].fY - p[0].fY);
335     klm[4] = w2 * (p[0].fX - p[1].fX);
336     klm[5] = w2 * (p[1].fX * p[0].fY - p[0].fX * p[1].fY);
337 
338     klm[6] = w2 * (p[2].fY - p[1].fY);
339     klm[7] = w2 * (p[1].fX - p[2].fX);
340     klm[8] = w2 * (p[2].fX * p[1].fY - p[1].fX * p[2].fY);
341 
342     // scale the max absolute value of coeffs to 10
343     SkScalar scale = 0.f;
344     for (int i = 0; i < 9; ++i) {
345        scale = SkMaxScalar(scale, SkScalarAbs(klm[i]));
346     }
347     SkASSERT(scale > 0.f);
348     scale = 10.f / scale;
349     for (int i = 0; i < 9; ++i) {
350         klm[i] *= scale;
351     }
352 }
353 
354 ////////////////////////////////////////////////////////////////////////////////
355 
356 namespace {
357 
358 // a is the first control point of the cubic.
359 // ab is the vector from a to the second control point.
360 // dc is the vector from the fourth to the third control point.
361 // d is the fourth control point.
362 // p is the candidate quadratic control point.
363 // this assumes that the cubic doesn't inflect and is simple
364 bool is_point_within_cubic_tangents(const SkPoint& a,
365                                     const SkVector& ab,
366                                     const SkVector& dc,
367                                     const SkPoint& d,
368                                     SkPathPriv::FirstDirection dir,
369                                     const SkPoint p) {
370     SkVector ap = p - a;
371     SkScalar apXab = ap.cross(ab);
372     if (SkPathPriv::kCW_FirstDirection == dir) {
373         if (apXab > 0) {
374             return false;
375         }
376     } else {
377         SkASSERT(SkPathPriv::kCCW_FirstDirection == dir);
378         if (apXab < 0) {
379             return false;
380         }
381     }
382 
383     SkVector dp = p - d;
384     SkScalar dpXdc = dp.cross(dc);
385     if (SkPathPriv::kCW_FirstDirection == dir) {
386         if (dpXdc < 0) {
387             return false;
388         }
389     } else {
390         SkASSERT(SkPathPriv::kCCW_FirstDirection == dir);
391         if (dpXdc > 0) {
392             return false;
393         }
394     }
395     return true;
396 }
397 
398 void convert_noninflect_cubic_to_quads(const SkPoint p[4],
399                                        SkScalar toleranceSqd,
400                                        SkTArray<SkPoint, true>* quads,
401                                        int sublevel = 0,
402                                        bool preserveFirstTangent = true,
403                                        bool preserveLastTangent = true) {
404     // Notation: Point a is always p[0]. Point b is p[1] unless p[1] == p[0], in which case it is
405     // p[2]. Point d is always p[3]. Point c is p[2] unless p[2] == p[3], in which case it is p[1].
406     SkVector ab = p[1] - p[0];
407     SkVector dc = p[2] - p[3];
408 
409     if (SkPointPriv::LengthSqd(ab) < SK_ScalarNearlyZero) {
410         if (SkPointPriv::LengthSqd(dc) < SK_ScalarNearlyZero) {
411             SkPoint* degQuad = quads->push_back_n(3);
412             degQuad[0] = p[0];
413             degQuad[1] = p[0];
414             degQuad[2] = p[3];
415             return;
416         }
417         ab = p[2] - p[0];
418     }
419     if (SkPointPriv::LengthSqd(dc) < SK_ScalarNearlyZero) {
420         dc = p[1] - p[3];
421     }
422 
423     static const SkScalar kLengthScale = 3 * SK_Scalar1 / 2;
424     static const int kMaxSubdivs = 10;
425 
426     ab.scale(kLengthScale);
427     dc.scale(kLengthScale);
428 
429     // c0 and c1 are extrapolations along vectors ab and dc.
430     SkPoint c0 = p[0] + ab;
431     SkPoint c1 = p[3] + dc;
432 
433     SkScalar dSqd = sublevel > kMaxSubdivs ? 0 : SkPointPriv::DistanceToSqd(c0, c1);
434     if (dSqd < toleranceSqd) {
435         SkPoint newC;
436         if (preserveFirstTangent == preserveLastTangent) {
437             // We used to force a split when both tangents need to be preserved and c0 != c1.
438             // This introduced a large performance regression for tiny paths for no noticeable
439             // quality improvement. However, we aren't quite fulfilling our contract of guaranteeing
440             // the two tangent vectors and this could introduce a missed pixel in
441             // GrAAHairlinePathRenderer.
442             newC = (c0 + c1) * 0.5f;
443         } else if (preserveFirstTangent) {
444             newC = c0;
445         } else {
446             newC = c1;
447         }
448 
449         SkPoint* pts = quads->push_back_n(3);
450         pts[0] = p[0];
451         pts[1] = newC;
452         pts[2] = p[3];
453         return;
454     }
455     SkPoint choppedPts[7];
456     SkChopCubicAtHalf(p, choppedPts);
457     convert_noninflect_cubic_to_quads(
458             choppedPts + 0, toleranceSqd, quads, sublevel + 1, preserveFirstTangent, false);
459     convert_noninflect_cubic_to_quads(
460             choppedPts + 3, toleranceSqd, quads, sublevel + 1, false, preserveLastTangent);
461 }
462 
463 void convert_noninflect_cubic_to_quads_with_constraint(const SkPoint p[4],
464                                                        SkScalar toleranceSqd,
465                                                        SkPathPriv::FirstDirection dir,
466                                                        SkTArray<SkPoint, true>* quads,
467                                                        int sublevel = 0) {
468     // Notation: Point a is always p[0]. Point b is p[1] unless p[1] == p[0], in which case it is
469     // p[2]. Point d is always p[3]. Point c is p[2] unless p[2] == p[3], in which case it is p[1].
470 
471     SkVector ab = p[1] - p[0];
472     SkVector dc = p[2] - p[3];
473 
474     if (SkPointPriv::LengthSqd(ab) < SK_ScalarNearlyZero) {
475         if (SkPointPriv::LengthSqd(dc) < SK_ScalarNearlyZero) {
476             SkPoint* degQuad = quads->push_back_n(3);
477             degQuad[0] = p[0];
478             degQuad[1] = p[0];
479             degQuad[2] = p[3];
480             return;
481         }
482         ab = p[2] - p[0];
483     }
484     if (SkPointPriv::LengthSqd(dc) < SK_ScalarNearlyZero) {
485         dc = p[1] - p[3];
486     }
487 
488     // When the ab and cd tangents are degenerate or nearly parallel with vector from d to a the
489     // constraint that the quad point falls between the tangents becomes hard to enforce and we are
490     // likely to hit the max subdivision count. However, in this case the cubic is approaching a
491     // line and the accuracy of the quad point isn't so important. We check if the two middle cubic
492     // control points are very close to the baseline vector. If so then we just pick quadratic
493     // points on the control polygon.
494 
495     SkVector da = p[0] - p[3];
496     bool doQuads = SkPointPriv::LengthSqd(dc) < SK_ScalarNearlyZero ||
497                    SkPointPriv::LengthSqd(ab) < SK_ScalarNearlyZero;
498     if (!doQuads) {
499         SkScalar invDALengthSqd = SkPointPriv::LengthSqd(da);
500         if (invDALengthSqd > SK_ScalarNearlyZero) {
501             invDALengthSqd = SkScalarInvert(invDALengthSqd);
502             // cross(ab, da)^2/length(da)^2 == sqd distance from b to line from d to a.
503             // same goes for point c using vector cd.
504             SkScalar detABSqd = ab.cross(da);
505             detABSqd = SkScalarSquare(detABSqd);
506             SkScalar detDCSqd = dc.cross(da);
507             detDCSqd = SkScalarSquare(detDCSqd);
508             if (detABSqd * invDALengthSqd < toleranceSqd &&
509                 detDCSqd * invDALengthSqd < toleranceSqd) {
510                 doQuads = true;
511             }
512         }
513     }
514     if (doQuads) {
515         SkPoint b = p[0] + ab;
516         SkPoint c = p[3] + dc;
517         SkPoint mid = b + c;
518         mid.scale(SK_ScalarHalf);
519         // Insert two quadratics to cover the case when ab points away from d and/or dc
520         // points away from a.
521         if (SkVector::DotProduct(da, dc) < 0 || SkVector::DotProduct(ab, da) > 0) {
522             SkPoint* qpts = quads->push_back_n(6);
523             qpts[0] = p[0];
524             qpts[1] = b;
525             qpts[2] = mid;
526             qpts[3] = mid;
527             qpts[4] = c;
528             qpts[5] = p[3];
529         } else {
530             SkPoint* qpts = quads->push_back_n(3);
531             qpts[0] = p[0];
532             qpts[1] = mid;
533             qpts[2] = p[3];
534         }
535         return;
536     }
537 
538     static const SkScalar kLengthScale = 3 * SK_Scalar1 / 2;
539     static const int kMaxSubdivs = 10;
540 
541     ab.scale(kLengthScale);
542     dc.scale(kLengthScale);
543 
544     // c0 and c1 are extrapolations along vectors ab and dc.
545     SkVector c0 = p[0] + ab;
546     SkVector c1 = p[3] + dc;
547 
548     SkScalar dSqd = sublevel > kMaxSubdivs ? 0 : SkPointPriv::DistanceToSqd(c0, c1);
549     if (dSqd < toleranceSqd) {
550         SkPoint cAvg = (c0 + c1) * 0.5f;
551         bool subdivide = false;
552 
553         if (!is_point_within_cubic_tangents(p[0], ab, dc, p[3], dir, cAvg)) {
554             // choose a new cAvg that is the intersection of the two tangent lines.
555             ab = SkPointPriv::MakeOrthog(ab);
556             SkScalar z0 = -ab.dot(p[0]);
557             dc = SkPointPriv::MakeOrthog(dc);
558             SkScalar z1 = -dc.dot(p[3]);
559             cAvg.fX = ab.fY * z1 - z0 * dc.fY;
560             cAvg.fY = z0 * dc.fX - ab.fX * z1;
561             SkScalar z = ab.fX * dc.fY - ab.fY * dc.fX;
562             z = SkScalarInvert(z);
563             cAvg.fX *= z;
564             cAvg.fY *= z;
565             if (sublevel <= kMaxSubdivs) {
566                 SkScalar d0Sqd = SkPointPriv::DistanceToSqd(c0, cAvg);
567                 SkScalar d1Sqd = SkPointPriv::DistanceToSqd(c1, cAvg);
568                 // We need to subdivide if d0 + d1 > tolerance but we have the sqd values. We know
569                 // the distances and tolerance can't be negative.
570                 // (d0 + d1)^2 > toleranceSqd
571                 // d0Sqd + 2*d0*d1 + d1Sqd > toleranceSqd
572                 SkScalar d0d1 = SkScalarSqrt(d0Sqd * d1Sqd);
573                 subdivide = 2 * d0d1 + d0Sqd + d1Sqd > toleranceSqd;
574             }
575         }
576         if (!subdivide) {
577             SkPoint* pts = quads->push_back_n(3);
578             pts[0] = p[0];
579             pts[1] = cAvg;
580             pts[2] = p[3];
581             return;
582         }
583     }
584     SkPoint choppedPts[7];
585     SkChopCubicAtHalf(p, choppedPts);
586     convert_noninflect_cubic_to_quads_with_constraint(
587             choppedPts + 0, toleranceSqd, dir, quads, sublevel + 1);
588     convert_noninflect_cubic_to_quads_with_constraint(
589             choppedPts + 3, toleranceSqd, dir, quads, sublevel + 1);
590 }
591 }
592 
593 void GrPathUtils::convertCubicToQuads(const SkPoint p[4],
594                                       SkScalar tolScale,
595                                       SkTArray<SkPoint, true>* quads) {
596     if (!p[0].isFinite() || !p[1].isFinite() || !p[2].isFinite() || !p[3].isFinite()) {
597         return;
598     }
599     if (!SkScalarIsFinite(tolScale)) {
600         return;
601     }
602     SkPoint chopped[10];
603     int count = SkChopCubicAtInflections(p, chopped);
604 
605     const SkScalar tolSqd = SkScalarSquare(tolScale);
606 
607     for (int i = 0; i < count; ++i) {
608         SkPoint* cubic = chopped + 3*i;
609         convert_noninflect_cubic_to_quads(cubic, tolSqd, quads);
610     }
611 }
612 
613 void GrPathUtils::convertCubicToQuadsConstrainToTangents(const SkPoint p[4],
614                                                          SkScalar tolScale,
615                                                          SkPathPriv::FirstDirection dir,
616                                                          SkTArray<SkPoint, true>* quads) {
617     if (!p[0].isFinite() || !p[1].isFinite() || !p[2].isFinite() || !p[3].isFinite()) {
618         return;
619     }
620     if (!SkScalarIsFinite(tolScale)) {
621         return;
622     }
623     SkPoint chopped[10];
624     int count = SkChopCubicAtInflections(p, chopped);
625 
626     const SkScalar tolSqd = SkScalarSquare(tolScale);
627 
628     for (int i = 0; i < count; ++i) {
629         SkPoint* cubic = chopped + 3*i;
630         convert_noninflect_cubic_to_quads_with_constraint(cubic, tolSqd, dir, quads);
631     }
632 }
633 
634 ////////////////////////////////////////////////////////////////////////////////
635 
636 using ExcludedTerm = GrPathUtils::ExcludedTerm;
637 
638 ExcludedTerm GrPathUtils::calcCubicInverseTransposePowerBasisMatrix(const SkPoint p[4],
639                                                                     SkMatrix* out) {
640     GR_STATIC_ASSERT(SK_SCALAR_IS_FLOAT);
641 
642     // First convert the bezier coordinates p[0..3] to power basis coefficients X,Y(,W=[0 0 0 1]).
643     // M3 is the matrix that does this conversion. The homogeneous equation for the cubic becomes:
644     //
645     //                                     | X   Y   0 |
646     // C(t,s) = [t^3  t^2*s  t*s^2  s^3] * | .   .   0 |
647     //                                     | .   .   0 |
648     //                                     | .   .   1 |
649     //
650     const Sk4f M3[3] = {Sk4f(-1, 3, -3, 1),
651                         Sk4f(3, -6, 3, 0),
652                         Sk4f(-3, 3, 0, 0)};
653     // 4th col of M3 =  Sk4f(1, 0, 0, 0)};
654     Sk4f X(p[3].x(), 0, 0, 0);
655     Sk4f Y(p[3].y(), 0, 0, 0);
656     for (int i = 2; i >= 0; --i) {
657         X += M3[i] * p[i].x();
658         Y += M3[i] * p[i].y();
659     }
660 
661     // The matrix is 3x4. In order to invert it, we first need to make it square by throwing out one
662     // of the middle two rows. We toss the row that leaves us with the largest absolute determinant.
663     // Since the right column will be [0 0 1], the respective determinants reduce to x0*y2 - y0*x2
664     // and x0*y1 - y0*x1.
665     SkScalar dets[4];
666     Sk4f D = SkNx_shuffle<0,0,2,1>(X) * SkNx_shuffle<2,1,0,0>(Y);
667     D -= SkNx_shuffle<2,3,0,1>(D);
668     D.store(dets);
669     ExcludedTerm skipTerm = SkScalarAbs(dets[0]) > SkScalarAbs(dets[1]) ?
670                             ExcludedTerm::kQuadraticTerm : ExcludedTerm::kLinearTerm;
671     SkScalar det = dets[ExcludedTerm::kQuadraticTerm == skipTerm ? 0 : 1];
672     if (0 == det) {
673         return ExcludedTerm::kNonInvertible;
674     }
675     SkScalar rdet = 1 / det;
676 
677     // Compute the inverse-transpose of the power basis matrix with the 'skipRow'th row removed.
678     // Since W=[0 0 0 1], it follows that our corresponding solution will be equal to:
679     //
680     //             |  y1  -x1   x1*y2 - y1*x2 |
681     //     1/det * | -y0   x0  -x0*y2 + y0*x2 |
682     //             |   0    0             det |
683     //
684     SkScalar x[4], y[4], z[4];
685     X.store(x);
686     Y.store(y);
687     (X * SkNx_shuffle<3,3,3,3>(Y) - Y * SkNx_shuffle<3,3,3,3>(X)).store(z);
688 
689     int middleRow = ExcludedTerm::kQuadraticTerm == skipTerm ? 2 : 1;
690     out->setAll( y[middleRow] * rdet, -x[middleRow] * rdet,  z[middleRow] * rdet,
691                         -y[0] * rdet,          x[0] * rdet,         -z[0] * rdet,
692                                    0,                    0,                    1);
693 
694     return skipTerm;
695 }
696 
697 inline static void calc_serp_kcoeffs(SkScalar tl, SkScalar sl, SkScalar tm, SkScalar sm,
698                                      ExcludedTerm skipTerm, SkScalar outCoeffs[3]) {
699     SkASSERT(ExcludedTerm::kQuadraticTerm == skipTerm || ExcludedTerm::kLinearTerm == skipTerm);
700     outCoeffs[0] = 0;
701     outCoeffs[1] = (ExcludedTerm::kLinearTerm == skipTerm) ? sl*sm : -tl*sm - tm*sl;
702     outCoeffs[2] = tl*tm;
703 }
704 
705 inline static void calc_serp_lmcoeffs(SkScalar t, SkScalar s, ExcludedTerm skipTerm,
706                                       SkScalar outCoeffs[3]) {
707     SkASSERT(ExcludedTerm::kQuadraticTerm == skipTerm || ExcludedTerm::kLinearTerm == skipTerm);
708     outCoeffs[0] = -s*s*s;
709     outCoeffs[1] = (ExcludedTerm::kLinearTerm == skipTerm) ? 3*s*s*t : -3*s*t*t;
710     outCoeffs[2] = t*t*t;
711 }
712 
713 inline static void calc_loop_kcoeffs(SkScalar td, SkScalar sd, SkScalar te, SkScalar se,
714                                      SkScalar tdse, SkScalar tesd, ExcludedTerm skipTerm,
715                                      SkScalar outCoeffs[3]) {
716     SkASSERT(ExcludedTerm::kQuadraticTerm == skipTerm || ExcludedTerm::kLinearTerm == skipTerm);
717     outCoeffs[0] = 0;
718     outCoeffs[1] = (ExcludedTerm::kLinearTerm == skipTerm) ? sd*se : -tdse - tesd;
719     outCoeffs[2] = td*te;
720 }
721 
722 inline static void calc_loop_lmcoeffs(SkScalar t2, SkScalar s2, SkScalar t1, SkScalar s1,
723                                       SkScalar t2s1, SkScalar t1s2, ExcludedTerm skipTerm,
724                                       SkScalar outCoeffs[3]) {
725     SkASSERT(ExcludedTerm::kQuadraticTerm == skipTerm || ExcludedTerm::kLinearTerm == skipTerm);
726     outCoeffs[0] = -s2*s2*s1;
727     outCoeffs[1] = (ExcludedTerm::kLinearTerm == skipTerm) ? s2 * (2*t2s1 + t1s2)
728                                                            : -t2 * (t2s1 + 2*t1s2);
729     outCoeffs[2] = t2*t2*t1;
730 }
731 
732 // For the case when a cubic bezier is actually a quadratic. We duplicate k in l so that the
733 // implicit becomes:
734 //
735 //     k^3 - l*m == k^3 - l*k == k * (k^2 - l)
736 //
737 // In the quadratic case we can simply assign fixed values at each control point:
738 //
739 //     | ..K.. |     | pts[0]  pts[1]  pts[2]  pts[3] |      | 0   1/3  2/3  1 |
740 //     | ..L.. |  *  |   .       .       .       .    |  ==  | 0     0  1/3  1 |
741 //     | ..K.. |     |   1       1       1       1    |      | 0   1/3  2/3  1 |
742 //
743 static void calc_quadratic_klm(const SkPoint pts[4], double d3, SkMatrix* klm) {
744     SkMatrix klmAtPts;
745     klmAtPts.setAll(0,  1.f/3,  1,
746                     0,      0,  1,
747                     0,  1.f/3,  1);
748 
749     SkMatrix inversePts;
750     inversePts.setAll(pts[0].x(),  pts[1].x(),  pts[3].x(),
751                       pts[0].y(),  pts[1].y(),  pts[3].y(),
752                                1,           1,           1);
753     SkAssertResult(inversePts.invert(&inversePts));
754 
755     klm->setConcat(klmAtPts, inversePts);
756 
757     // If d3 > 0 we need to flip the orientation of our curve
758     // This is done by negating the k and l values
759     if (d3 > 0) {
760         klm->postScale(-1, -1);
761     }
762 }
763 
764 // For the case when a cubic bezier is actually a line. We set K=0, L=1, M=-line, which results in
765 // the following implicit:
766 //
767 //     k^3 - l*m == 0^3 - 1*(-line) == -(-line) == line
768 //
769 static void calc_line_klm(const SkPoint pts[4], SkMatrix* klm) {
770     SkScalar ny = pts[0].x() - pts[3].x();
771     SkScalar nx = pts[3].y() - pts[0].y();
772     SkScalar k = nx * pts[0].x() + ny * pts[0].y();
773     klm->setAll(  0,   0, 0,
774                   0,   0, 1,
775                 -nx, -ny, k);
776 }
777 
778 SkCubicType GrPathUtils::getCubicKLM(const SkPoint src[4], SkMatrix* klm, double tt[2],
779                                      double ss[2]) {
780     double d[4];
781     SkCubicType type = SkClassifyCubic(src, tt, ss, d);
782 
783     if (SkCubicType::kLineOrPoint == type) {
784         calc_line_klm(src, klm);
785         return SkCubicType::kLineOrPoint;
786     }
787 
788     if (SkCubicType::kQuadratic == type) {
789         calc_quadratic_klm(src, d[3], klm);
790         return SkCubicType::kQuadratic;
791     }
792 
793     SkMatrix CIT;
794     ExcludedTerm skipTerm = calcCubicInverseTransposePowerBasisMatrix(src, &CIT);
795     if (ExcludedTerm::kNonInvertible == skipTerm) {
796         // This could technically also happen if the curve were quadratic, but SkClassifyCubic
797         // should have detected that case already with tolerance.
798         calc_line_klm(src, klm);
799         return SkCubicType::kLineOrPoint;
800     }
801 
802     const SkScalar t0 = static_cast<SkScalar>(tt[0]), t1 = static_cast<SkScalar>(tt[1]),
803                    s0 = static_cast<SkScalar>(ss[0]), s1 = static_cast<SkScalar>(ss[1]);
804 
805     SkMatrix klmCoeffs;
806     switch (type) {
807         case SkCubicType::kCuspAtInfinity:
808             SkASSERT(1 == t1 && 0 == s1); // Infinity.
809             // fallthru.
810         case SkCubicType::kLocalCusp:
811         case SkCubicType::kSerpentine:
812             calc_serp_kcoeffs(t0, s0, t1, s1, skipTerm, &klmCoeffs[0]);
813             calc_serp_lmcoeffs(t0, s0, skipTerm, &klmCoeffs[3]);
814             calc_serp_lmcoeffs(t1, s1, skipTerm, &klmCoeffs[6]);
815             break;
816         case SkCubicType::kLoop: {
817             const SkScalar tdse = t0 * s1;
818             const SkScalar tesd = t1 * s0;
819             calc_loop_kcoeffs(t0, s0, t1, s1, tdse, tesd, skipTerm, &klmCoeffs[0]);
820             calc_loop_lmcoeffs(t0, s0, t1, s1, tdse, tesd, skipTerm, &klmCoeffs[3]);
821             calc_loop_lmcoeffs(t1, s1, t0, s0, tesd, tdse, skipTerm, &klmCoeffs[6]);
822             break;
823         }
824         default:
825             SK_ABORT("Unexpected cubic type.");
826             break;
827     }
828 
829     klm->setConcat(klmCoeffs, CIT);
830     return type;
831 }
832 
833 int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[10], SkMatrix* klm,
834                                              int* loopIndex) {
835     SkSTArray<2, SkScalar> chops;
836     *loopIndex = -1;
837 
838     double t[2], s[2];
839     if (SkCubicType::kLoop == GrPathUtils::getCubicKLM(src, klm, t, s)) {
840         SkScalar t0 = static_cast<SkScalar>(t[0] / s[0]);
841         SkScalar t1 = static_cast<SkScalar>(t[1] / s[1]);
842         SkASSERT(t0 <= t1); // Technically t0 != t1 in a loop, but there may be FP error.
843 
844         if (t0 < 1 && t1 > 0) {
845             *loopIndex = 0;
846             if (t0 > 0) {
847                 chops.push_back(t0);
848                 *loopIndex = 1;
849             }
850             if (t1 < 1) {
851                 chops.push_back(t1);
852                 *loopIndex = chops.count() - 1;
853             }
854         }
855     }
856 
857     SkChopCubicAt(src, dst, chops.begin(), chops.count());
858     return chops.count() + 1;
859 }
860