• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * Copyright 2017 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 "GrCCFillGeometry.h"
9 
10 #include "GrTypes.h"
11 #include "SkGeometry.h"
12 #include <algorithm>
13 #include <cmath>
14 #include <cstdlib>
15 
16 static constexpr float kFlatnessThreshold = 1/16.f; // 1/16 of a pixel.
17 
beginPath()18 void GrCCFillGeometry::beginPath() {
19     SkASSERT(!fBuildingContour);
20     fVerbs.push_back(Verb::kBeginPath);
21 }
22 
beginContour(const SkPoint & pt)23 void GrCCFillGeometry::beginContour(const SkPoint& pt) {
24     SkASSERT(!fBuildingContour);
25     // Store the current verb count in the fTriangles field for now. When we close the contour we
26     // will use this value to calculate the actual number of triangles in its fan.
27     fCurrContourTallies = {fVerbs.count(), 0, 0, 0, 0};
28 
29     fPoints.push_back(pt);
30     fVerbs.push_back(Verb::kBeginContour);
31     fCurrAnchorPoint = pt;
32 
33     SkDEBUGCODE(fBuildingContour = true);
34 }
35 
lineTo(const SkPoint P[2])36 void GrCCFillGeometry::lineTo(const SkPoint P[2]) {
37     SkASSERT(fBuildingContour);
38     SkASSERT(P[0] == fPoints.back());
39     Sk2f p0 = Sk2f::Load(P);
40     Sk2f p1 = Sk2f::Load(P+1);
41     this->appendLine(p0, p1);
42 }
43 
appendLine(const Sk2f & p0,const Sk2f & p1)44 inline void GrCCFillGeometry::appendLine(const Sk2f& p0, const Sk2f& p1) {
45     SkASSERT(fPoints.back() == SkPoint::Make(p0[0], p0[1]));
46     if ((p0 == p1).allTrue()) {
47         return;
48     }
49     p1.store(&fPoints.push_back());
50     fVerbs.push_back(Verb::kLineTo);
51 }
52 
normalize(const Sk2f & n)53 static inline Sk2f normalize(const Sk2f& n) {
54     Sk2f nn = n*n;
55     return n * (nn + SkNx_shuffle<1,0>(nn)).rsqrt();
56 }
57 
dot(const Sk2f & a,const Sk2f & b)58 static inline float dot(const Sk2f& a, const Sk2f& b) {
59     float product[2];
60     (a * b).store(product);
61     return product[0] + product[1];
62 }
63 
are_collinear(const Sk2f & p0,const Sk2f & p1,const Sk2f & p2,float tolerance=kFlatnessThreshold)64 static inline bool are_collinear(const Sk2f& p0, const Sk2f& p1, const Sk2f& p2,
65                                  float tolerance = kFlatnessThreshold) {
66     Sk2f l = p2 - p0; // Line from p0 -> p2.
67 
68     // lwidth = Manhattan width of l.
69     Sk2f labs = l.abs();
70     float lwidth = labs[0] + labs[1];
71 
72     // d = |p1 - p0| dot | l.y|
73     //                   |-l.x| = distance from p1 to l.
74     Sk2f dd = (p1 - p0) * SkNx_shuffle<1,0>(l);
75     float d = dd[0] - dd[1];
76 
77     // We are collinear if a box with radius "tolerance", centered on p1, touches the line l.
78     // To decide this, we check if the distance from p1 to the line is less than the distance from
79     // p1 to the far corner of this imaginary box, along that same normal vector.
80     // The far corner of the box can be found at "p1 + sign(n) * tolerance", where n is normal to l:
81     //
82     //   abs(dot(p1 - p0, n)) <= dot(sign(n) * tolerance, n)
83     //
84     // Which reduces to:
85     //
86     //   abs(d) <= (n.x * sign(n.x) + n.y * sign(n.y)) * tolerance
87     //   abs(d) <= (abs(n.x) + abs(n.y)) * tolerance
88     //
89     // Use "<=" in case l == 0.
90     return std::abs(d) <= lwidth * tolerance;
91 }
92 
are_collinear(const SkPoint P[4],float tolerance=kFlatnessThreshold)93 static inline bool are_collinear(const SkPoint P[4], float tolerance = kFlatnessThreshold) {
94     Sk4f Px, Py;               // |Px  Py|   |p0 - p3|
95     Sk4f::Load2(P, &Px, &Py);  // |.   . | = |p1 - p3|
96     Px -= Px[3];               // |.   . |   |p2 - p3|
97     Py -= Py[3];               // |.   . |   |   0   |
98 
99     // Find [lx, ly] = the line from p3 to the furthest-away point from p3.
100     Sk4f Pwidth = Px.abs() + Py.abs(); // Pwidth = Manhattan width of each point.
101     int lidx = Pwidth[0] > Pwidth[1] ? 0 : 1;
102     lidx = Pwidth[lidx] > Pwidth[2] ? lidx : 2;
103     float lx = Px[lidx], ly = Py[lidx];
104     float lwidth = Pwidth[lidx]; // lwidth = Manhattan width of [lx, ly].
105 
106     //     |Px  Py|
107     // d = |.   . | * | ly| = distances from each point to l (two of the distances will be zero).
108     //     |.   . |   |-lx|
109     //     |.   . |
110     Sk4f d = Px*ly - Py*lx;
111 
112     // We are collinear if boxes with radius "tolerance", centered on all 4 points all touch line l.
113     // (See the rationale for this formula in the above, 3-point version of this function.)
114     // Use "<=" in case l == 0.
115     return (d.abs() <= lwidth * tolerance).allTrue();
116 }
117 
118 // Returns whether the (convex) curve segment is monotonic with respect to [endPt - startPt].
is_convex_curve_monotonic(const Sk2f & startPt,const Sk2f & tan0,const Sk2f & endPt,const Sk2f & tan1)119 static inline bool is_convex_curve_monotonic(const Sk2f& startPt, const Sk2f& tan0,
120                                              const Sk2f& endPt, const Sk2f& tan1) {
121     Sk2f v = endPt - startPt;
122     float dot0 = dot(tan0, v);
123     float dot1 = dot(tan1, v);
124 
125     // A small, negative tolerance handles floating-point error in the case when one tangent
126     // approaches 0 length, meaning the (convex) curve segment is effectively a flat line.
127     float tolerance = -std::max(std::abs(dot0), std::abs(dot1)) * SK_ScalarNearlyZero;
128     return dot0 >= tolerance && dot1 >= tolerance;
129 }
130 
lerp(const SkNx<N,float> & a,const SkNx<N,float> & b,const SkNx<N,float> & t)131 template<int N> static inline SkNx<N,float> lerp(const SkNx<N,float>& a, const SkNx<N,float>& b,
132                                                  const SkNx<N,float>& t) {
133     return SkNx_fma(t, b - a, a);
134 }
135 
quadraticTo(const SkPoint P[3])136 void GrCCFillGeometry::quadraticTo(const SkPoint P[3]) {
137     SkASSERT(fBuildingContour);
138     SkASSERT(P[0] == fPoints.back());
139     Sk2f p0 = Sk2f::Load(P);
140     Sk2f p1 = Sk2f::Load(P+1);
141     Sk2f p2 = Sk2f::Load(P+2);
142 
143     // Don't crunch on the curve if it is nearly flat (or just very small). Flat curves can break
144     // The monotonic chopping math.
145     if (are_collinear(p0, p1, p2)) {
146         this->appendLine(p0, p2);
147         return;
148     }
149 
150     this->appendQuadratics(p0, p1, p2);
151 }
152 
appendQuadratics(const Sk2f & p0,const Sk2f & p1,const Sk2f & p2)153 inline void GrCCFillGeometry::appendQuadratics(const Sk2f& p0, const Sk2f& p1, const Sk2f& p2) {
154     Sk2f tan0 = p1 - p0;
155     Sk2f tan1 = p2 - p1;
156 
157     // This should almost always be this case for well-behaved curves in the real world.
158     if (is_convex_curve_monotonic(p0, tan0, p2, tan1)) {
159         this->appendMonotonicQuadratic(p0, p1, p2);
160         return;
161     }
162 
163     // Chop the curve into two segments with equal curvature. To do this we find the T value whose
164     // tangent angle is halfway between tan0 and tan1.
165     Sk2f n = normalize(tan0) - normalize(tan1);
166 
167     // The midtangent can be found where (dQ(t) dot n) = 0:
168     //
169     //   0 = (dQ(t) dot n) = | 2*t  1 | * | p0 - 2*p1 + p2 | * | n |
170     //                                    | -2*p0 + 2*p1   |   | . |
171     //
172     //                     = | 2*t  1 | * | tan1 - tan0 | * | n |
173     //                                    | 2*tan0      |   | . |
174     //
175     //                     = 2*t * ((tan1 - tan0) dot n) + (2*tan0 dot n)
176     //
177     //   t = (tan0 dot n) / ((tan0 - tan1) dot n)
178     Sk2f dQ1n = (tan0 - tan1) * n;
179     Sk2f dQ0n = tan0 * n;
180     Sk2f t = (dQ0n + SkNx_shuffle<1,0>(dQ0n)) / (dQ1n + SkNx_shuffle<1,0>(dQ1n));
181     t = Sk2f::Min(Sk2f::Max(t, 0), 1); // Clamp for FP error.
182 
183     Sk2f p01 = SkNx_fma(t, tan0, p0);
184     Sk2f p12 = SkNx_fma(t, tan1, p1);
185     Sk2f p012 = lerp(p01, p12, t);
186 
187     this->appendMonotonicQuadratic(p0, p01, p012);
188     this->appendMonotonicQuadratic(p012, p12, p2);
189 }
190 
appendMonotonicQuadratic(const Sk2f & p0,const Sk2f & p1,const Sk2f & p2)191 inline void GrCCFillGeometry::appendMonotonicQuadratic(const Sk2f& p0, const Sk2f& p1,
192                                                        const Sk2f& p2) {
193     // Don't send curves to the GPU if we know they are nearly flat (or just very small).
194     if (are_collinear(p0, p1, p2)) {
195         this->appendLine(p0, p2);
196         return;
197     }
198 
199     SkASSERT(fPoints.back() == SkPoint::Make(p0[0], p0[1]));
200     SkASSERT((p0 != p2).anyTrue());
201     p1.store(&fPoints.push_back());
202     p2.store(&fPoints.push_back());
203     fVerbs.push_back(Verb::kMonotonicQuadraticTo);
204     ++fCurrContourTallies.fQuadratics;
205 }
206 
first_unless_nearly_zero(const Sk2f & a,const Sk2f & b)207 static inline Sk2f first_unless_nearly_zero(const Sk2f& a, const Sk2f& b) {
208     Sk2f aa = a*a;
209     aa += SkNx_shuffle<1,0>(aa);
210     SkASSERT(aa[0] == aa[1]);
211 
212     Sk2f bb = b*b;
213     bb += SkNx_shuffle<1,0>(bb);
214     SkASSERT(bb[0] == bb[1]);
215 
216     return (aa > bb * SK_ScalarNearlyZero).thenElse(a, b);
217 }
218 
get_cubic_tangents(const Sk2f & p0,const Sk2f & p1,const Sk2f & p2,const Sk2f & p3,Sk2f * tan0,Sk2f * tan1)219 static inline void get_cubic_tangents(const Sk2f& p0, const Sk2f& p1, const Sk2f& p2,
220                                       const Sk2f& p3, Sk2f* tan0, Sk2f* tan1) {
221     *tan0 = first_unless_nearly_zero(p1 - p0, p2 - p0);
222     *tan1 = first_unless_nearly_zero(p3 - p2, p3 - p1);
223 }
224 
is_cubic_nearly_quadratic(const Sk2f & p0,const Sk2f & p1,const Sk2f & p2,const Sk2f & p3,const Sk2f & tan0,const Sk2f & tan1,Sk2f * c)225 static inline bool is_cubic_nearly_quadratic(const Sk2f& p0, const Sk2f& p1, const Sk2f& p2,
226                                              const Sk2f& p3, const Sk2f& tan0, const Sk2f& tan1,
227                                              Sk2f* c) {
228     Sk2f c1 = SkNx_fma(Sk2f(1.5f), tan0, p0);
229     Sk2f c2 = SkNx_fma(Sk2f(-1.5f), tan1, p3);
230     *c = (c1 + c2) * .5f; // Hopefully optimized out if not used?
231     return ((c1 - c2).abs() <= 1).allTrue();
232 }
233 
234 enum class ExcludedTerm : bool {
235     kQuadraticTerm,
236     kLinearTerm
237 };
238 
239 // Finds where to chop a non-loop around its inflection points. The resulting cubic segments will be
240 // chopped such that a box of radius 'padRadius', centered at any point along the curve segment, is
241 // guaranteed to not cross the tangent lines at the inflection points (a.k.a lines L & M).
242 //
243 // 'chops' will be filled with 0, 2, or 4 T values. The segments between T0..T1 and T2..T3 must be
244 // drawn with flat lines instead of cubics.
245 //
246 // A serpentine cubic has two inflection points, so this method takes Sk2f and computes the padding
247 // for both in SIMD.
find_chops_around_inflection_points(float padRadius,Sk2f tl,Sk2f sl,const Sk2f & C0,const Sk2f & C1,ExcludedTerm skipTerm,float Cdet,SkSTArray<4,float> * chops)248 static inline void find_chops_around_inflection_points(float padRadius, Sk2f tl, Sk2f sl,
249                                                        const Sk2f& C0, const Sk2f& C1,
250                                                        ExcludedTerm skipTerm, float Cdet,
251                                                        SkSTArray<4, float>* chops) {
252     SkASSERT(chops->empty());
253     SkASSERT(padRadius >= 0);
254 
255     padRadius /= std::abs(Cdet); // Scale this single value rather than all of C^-1 later on.
256 
257     // The homogeneous parametric functions for distance from lines L & M are:
258     //
259     //     l(t,s) = (t*sl - s*tl)^3
260     //     m(t,s) = (t*sm - s*tm)^3
261     //
262     // See "Resolution Independent Curve Rendering using Programmable Graphics Hardware",
263     // 4.3 Finding klmn:
264     //
265     // https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
266     //
267     // From here on we use Sk2f with "L" names, but the second lane will be for line M.
268     tl = (sl > 0).thenElse(tl, -tl); // Tl=tl/sl is the triple root of l(t,s). Normalize so s >= 0.
269     sl = sl.abs();
270 
271     // Convert l(t,s), m(t,s) to power-basis form:
272     //
273     //                                                  | l3  m3 |
274     //    |l(t,s)  m(t,s)| = |t^3  t^2*s  t*s^2  s^3| * | l2  m2 |
275     //                                                  | l1  m1 |
276     //                                                  | l0  m0 |
277     //
278     Sk2f l3 = sl*sl*sl;
279     Sk2f l2or1 = (ExcludedTerm::kLinearTerm == skipTerm) ? sl*sl*tl*-3 : sl*tl*tl*3;
280 
281     // The equation for line L can be found as follows:
282     //
283     //     L = C^-1 * (l excluding skipTerm)
284     //
285     // (See comments for GrPathUtils::calcCubicInverseTransposePowerBasisMatrix.)
286     // We are only interested in the normal to L, so only need the upper 2x2 of C^-1. And rather
287     // than divide by determinant(C) here, we have already performed this divide on padRadius.
288     Sk2f Lx =  C1[1]*l3 - C0[1]*l2or1;
289     Sk2f Ly = -C1[0]*l3 + C0[0]*l2or1;
290 
291     // A box of radius "padRadius" is touching line L if "center dot L" is less than the Manhattan
292     // with of L. (See rationale in are_collinear.)
293     Sk2f Lwidth = Lx.abs() + Ly.abs();
294     Sk2f pad = Lwidth * padRadius;
295 
296     // Will T=(t + cbrt(pad))/s be greater than 0? No need to solve roots outside T=0..1.
297     Sk2f insideLeftPad = pad + tl*tl*tl;
298 
299     // Will T=(t - cbrt(pad))/s be less than 1? No need to solve roots outside T=0..1.
300     Sk2f tms = tl - sl;
301     Sk2f insideRightPad = pad - tms*tms*tms;
302 
303     // Solve for the T values where abs(l(T)) = pad.
304     if (insideLeftPad[0] > 0 && insideRightPad[0] > 0) {
305         float padT = cbrtf(pad[0]);
306         Sk2f pts = (tl[0] + Sk2f(-padT, +padT)) / sl[0];
307         pts.store(chops->push_back_n(2));
308     }
309 
310     // Solve for the T values where abs(m(T)) = pad.
311     if (insideLeftPad[1] > 0 && insideRightPad[1] > 0) {
312         float padT = cbrtf(pad[1]);
313         Sk2f pts = (tl[1] + Sk2f(-padT, +padT)) / sl[1];
314         pts.store(chops->push_back_n(2));
315     }
316 }
317 
swap_if_greater(float & a,float & b)318 static inline void swap_if_greater(float& a, float& b) {
319     if (a > b) {
320         std::swap(a, b);
321     }
322 }
323 
324 // Finds where to chop a non-loop around its intersection point. The resulting cubic segments will
325 // be chopped such that a box of radius 'padRadius', centered at any point along the curve segment,
326 // is guaranteed to not cross the tangent lines at the intersection point (a.k.a lines L & M).
327 //
328 // 'chops' will be filled with 0, 2, or 4 T values. The segments between T0..T1 and T2..T3 must be
329 // drawn with quadratic splines instead of cubics.
330 //
331 // A loop intersection falls at two different T values, so this method takes Sk2f and computes the
332 // padding for both in SIMD.
find_chops_around_loop_intersection(float padRadius,Sk2f t2,Sk2f s2,const Sk2f & C0,const Sk2f & C1,ExcludedTerm skipTerm,float Cdet,SkSTArray<4,float> * chops)333 static inline void find_chops_around_loop_intersection(float padRadius, Sk2f t2, Sk2f s2,
334                                                        const Sk2f& C0, const Sk2f& C1,
335                                                        ExcludedTerm skipTerm, float Cdet,
336                                                        SkSTArray<4, float>* chops) {
337     SkASSERT(chops->empty());
338     SkASSERT(padRadius >= 0);
339 
340     padRadius /= std::abs(Cdet); // Scale this single value rather than all of C^-1 later on.
341 
342     // The parametric functions for distance from lines L & M are:
343     //
344     //     l(T) = (T - Td)^2 * (T - Te)
345     //     m(T) = (T - Td) * (T - Te)^2
346     //
347     // See "Resolution Independent Curve Rendering using Programmable Graphics Hardware",
348     // 4.3 Finding klmn:
349     //
350     // https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
351     Sk2f T2 = t2/s2; // T2 is the double root of l(T).
352     Sk2f T1 = SkNx_shuffle<1,0>(T2); // T1 is the other root of l(T).
353 
354     // Convert l(T), m(T) to power-basis form:
355     //
356     //                                      |  1   1 |
357     //    |l(T)  m(T)| = |T^3  T^2  T  1| * | l2  m2 |
358     //                                      | l1  m1 |
359     //                                      | l0  m0 |
360     //
361     // From here on we use Sk2f with "L" names, but the second lane will be for line M.
362     Sk2f l2 = SkNx_fma(Sk2f(-2), T2, -T1);
363     Sk2f l1 = T2 * SkNx_fma(Sk2f(2), T1, T2);
364     Sk2f l0 = -T2*T2*T1;
365 
366     // The equation for line L can be found as follows:
367     //
368     //     L = C^-1 * (l excluding skipTerm)
369     //
370     // (See comments for GrPathUtils::calcCubicInverseTransposePowerBasisMatrix.)
371     // We are only interested in the normal to L, so only need the upper 2x2 of C^-1. And rather
372     // than divide by determinant(C) here, we have already performed this divide on padRadius.
373     Sk2f l2or1 = (ExcludedTerm::kLinearTerm == skipTerm) ? l2 : l1;
374     Sk2f Lx = -C0[1]*l2or1 + C1[1]; // l3 is always 1.
375     Sk2f Ly =  C0[0]*l2or1 - C1[0];
376 
377     // A box of radius "padRadius" is touching line L if "center dot L" is less than the Manhattan
378     // with of L. (See rationale in are_collinear.)
379     Sk2f Lwidth = Lx.abs() + Ly.abs();
380     Sk2f pad = Lwidth * padRadius;
381 
382     // Is l(T=0) outside the padding around line L?
383     Sk2f lT0 = l0; // l(T=0) = |0  0  0  1| dot |1  l2  l1  l0| = l0
384     Sk2f outsideT0 = lT0.abs() - pad;
385 
386     // Is l(T=1) outside the padding around line L?
387     Sk2f lT1 = (Sk2f(1) + l2 + l1 + l0).abs(); // l(T=1) = |1  1  1  1| dot |1  l2  l1  l0|
388     Sk2f outsideT1 = lT1.abs() - pad;
389 
390     // Values for solving the cubic.
391     Sk2f p, q, qqq, discr, numRoots, D;
392     bool hasDiscr = false;
393 
394     // Values for calculating one root (rarely needed).
395     Sk2f R, QQ;
396     bool hasOneRootVals = false;
397 
398     // Values for calculating three roots.
399     Sk2f P, cosTheta3;
400     bool hasThreeRootVals = false;
401 
402     // Solve for the T values where l(T) = +pad and m(T) = -pad.
403     for (int i = 0; i < 2; ++i) {
404         float T = T2[i]; // T is the point we are chopping around.
405         if ((T < 0 && outsideT0[i] >= 0) || (T > 1 && outsideT1[i] >= 0)) {
406             // The padding around T is completely out of range. No point solving for it.
407             continue;
408         }
409 
410         if (!hasDiscr) {
411             p = Sk2f(+.5f, -.5f) * pad;
412             q = (1.f/3) * (T2 - T1);
413             qqq = q*q*q;
414             discr = qqq*p*2 + p*p;
415             numRoots = (discr < 0).thenElse(3, 1);
416             D = T2 - q;
417             hasDiscr = true;
418         }
419 
420         if (1 == numRoots[i]) {
421             if (!hasOneRootVals) {
422                 Sk2f r = qqq + p;
423                 Sk2f s = r.abs() + discr.sqrt();
424                 R = (r > 0).thenElse(-s, s);
425                 QQ = q*q;
426                 hasOneRootVals = true;
427             }
428 
429             float A = cbrtf(R[i]);
430             float B = A != 0 ? QQ[i]/A : 0;
431             // When there is only one root, ine L chops from root..1, line M chops from 0..root.
432             if (1 == i) {
433                 chops->push_back(0);
434             }
435             chops->push_back(A + B + D[i]);
436             if (0 == i) {
437                 chops->push_back(1);
438             }
439             continue;
440         }
441 
442         if (!hasThreeRootVals) {
443             P = q.abs() * -2;
444             cosTheta3 = (q >= 0).thenElse(1, -1) + p / qqq.abs();
445             hasThreeRootVals = true;
446         }
447 
448         static constexpr float k2PiOver3 = 2 * SK_ScalarPI / 3;
449         float theta = std::acos(cosTheta3[i]) * (1.f/3);
450         float roots[3] = {P[i] * std::cos(theta) + D[i],
451                           P[i] * std::cos(theta + k2PiOver3) + D[i],
452                           P[i] * std::cos(theta - k2PiOver3) + D[i]};
453 
454         // Sort the three roots.
455         swap_if_greater(roots[0], roots[1]);
456         swap_if_greater(roots[1], roots[2]);
457         swap_if_greater(roots[0], roots[1]);
458 
459         // Line L chops around the first 2 roots, line M chops around the second 2.
460         chops->push_back_n(2, &roots[i]);
461     }
462 }
463 
cubicTo(const SkPoint P[4],float inflectPad,float loopIntersectPad)464 void GrCCFillGeometry::cubicTo(const SkPoint P[4], float inflectPad, float loopIntersectPad) {
465     SkASSERT(fBuildingContour);
466     SkASSERT(P[0] == fPoints.back());
467 
468     // Don't crunch on the curve or inflate geometry if it is nearly flat (or just very small).
469     // Flat curves can break the math below.
470     if (are_collinear(P)) {
471         Sk2f p0 = Sk2f::Load(P);
472         Sk2f p3 = Sk2f::Load(P+3);
473         this->appendLine(p0, p3);
474         return;
475     }
476 
477     Sk2f p0 = Sk2f::Load(P);
478     Sk2f p1 = Sk2f::Load(P+1);
479     Sk2f p2 = Sk2f::Load(P+2);
480     Sk2f p3 = Sk2f::Load(P+3);
481 
482     // Also detect near-quadratics ahead of time.
483     Sk2f tan0, tan1, c;
484     get_cubic_tangents(p0, p1, p2, p3, &tan0, &tan1);
485     if (is_cubic_nearly_quadratic(p0, p1, p2, p3, tan0, tan1, &c)) {
486         this->appendQuadratics(p0, c, p3);
487         return;
488     }
489 
490     double tt[2], ss[2], D[4];
491     fCurrCubicType = SkClassifyCubic(P, tt, ss, D);
492     SkASSERT(!SkCubicIsDegenerate(fCurrCubicType));
493     Sk2f t = Sk2f(static_cast<float>(tt[0]), static_cast<float>(tt[1]));
494     Sk2f s = Sk2f(static_cast<float>(ss[0]), static_cast<float>(ss[1]));
495 
496     ExcludedTerm skipTerm = (std::abs(D[2]) > std::abs(D[1]))
497                                     ? ExcludedTerm::kQuadraticTerm
498                                     : ExcludedTerm::kLinearTerm;
499     Sk2f C0 = SkNx_fma(Sk2f(3), p1 - p2, p3 - p0);
500     Sk2f C1 = (ExcludedTerm::kLinearTerm == skipTerm
501                        ? SkNx_fma(Sk2f(-2), p1, p0 + p2)
502                        : p1 - p0) * 3;
503     Sk2f C0x1 = C0 * SkNx_shuffle<1,0>(C1);
504     float Cdet = C0x1[0] - C0x1[1];
505 
506     SkSTArray<4, float> chops;
507     if (SkCubicType::kLoop != fCurrCubicType) {
508         find_chops_around_inflection_points(inflectPad, t, s, C0, C1, skipTerm, Cdet, &chops);
509     } else {
510         find_chops_around_loop_intersection(loopIntersectPad, t, s, C0, C1, skipTerm, Cdet, &chops);
511     }
512     if (4 == chops.count() && chops[1] >= chops[2]) {
513         // This just the means the KLM roots are so close that their paddings overlap. We will
514         // approximate the entire middle section, but still have it chopped midway. For loops this
515         // chop guarantees the append code only sees convex segments. Otherwise, it means we are (at
516         // least almost) a cusp and the chop makes sure we get a sharp point.
517         Sk2f ts = t * SkNx_shuffle<1,0>(s);
518         chops[1] = chops[2] = (ts[0] + ts[1]) / (2*s[0]*s[1]);
519     }
520 
521 #ifdef SK_DEBUG
522     for (int i = 1; i < chops.count(); ++i) {
523         SkASSERT(chops[i] >= chops[i - 1]);
524     }
525 #endif
526     this->appendCubics(AppendCubicMode::kLiteral, p0, p1, p2, p3, chops.begin(), chops.count());
527 }
528 
chop_cubic(const Sk2f & p0,const Sk2f & p1,const Sk2f & p2,const Sk2f & p3,float T,Sk2f * ab,Sk2f * abc,Sk2f * abcd,Sk2f * bcd,Sk2f * cd)529 static inline void chop_cubic(const Sk2f& p0, const Sk2f& p1, const Sk2f& p2, const Sk2f& p3,
530                               float T, Sk2f* ab, Sk2f* abc, Sk2f* abcd, Sk2f* bcd, Sk2f* cd) {
531     Sk2f TT = T;
532     *ab = lerp(p0, p1, TT);
533     Sk2f bc = lerp(p1, p2, TT);
534     *cd = lerp(p2, p3, TT);
535     *abc = lerp(*ab, bc, TT);
536     *bcd = lerp(bc, *cd, TT);
537     *abcd = lerp(*abc, *bcd, TT);
538 }
539 
appendCubics(AppendCubicMode mode,const Sk2f & p0,const Sk2f & p1,const Sk2f & p2,const Sk2f & p3,const float chops[],int numChops,float localT0,float localT1)540 void GrCCFillGeometry::appendCubics(AppendCubicMode mode, const Sk2f& p0, const Sk2f& p1,
541                                     const Sk2f& p2, const Sk2f& p3, const float chops[],
542                                     int numChops, float localT0, float localT1) {
543     if (numChops) {
544         SkASSERT(numChops > 0);
545         int midChopIdx = numChops/2;
546         float T = chops[midChopIdx];
547         // Chops alternate between literal and approximate mode.
548         AppendCubicMode rightMode = (AppendCubicMode)((bool)mode ^ (midChopIdx & 1) ^ 1);
549 
550         if (T <= localT0) {
551             // T is outside 0..1. Append the right side only.
552             this->appendCubics(rightMode, p0, p1, p2, p3, &chops[midChopIdx + 1],
553                                numChops - midChopIdx - 1, localT0, localT1);
554             return;
555         }
556 
557         if (T >= localT1) {
558             // T is outside 0..1. Append the left side only.
559             this->appendCubics(mode, p0, p1, p2, p3, chops, midChopIdx, localT0, localT1);
560             return;
561         }
562 
563         float localT = (T - localT0) / (localT1 - localT0);
564         Sk2f p01, p02, pT, p11, p12;
565         chop_cubic(p0, p1, p2, p3, localT, &p01, &p02, &pT, &p11, &p12);
566         this->appendCubics(mode, p0, p01, p02, pT, chops, midChopIdx, localT0, T);
567         this->appendCubics(rightMode, pT, p11, p12, p3, &chops[midChopIdx + 1],
568                            numChops - midChopIdx - 1, T, localT1);
569         return;
570     }
571 
572     this->appendCubics(mode, p0, p1, p2, p3);
573 }
574 
appendCubics(AppendCubicMode mode,const Sk2f & p0,const Sk2f & p1,const Sk2f & p2,const Sk2f & p3,int maxSubdivisions)575 void GrCCFillGeometry::appendCubics(AppendCubicMode mode, const Sk2f& p0, const Sk2f& p1,
576                                     const Sk2f& p2, const Sk2f& p3, int maxSubdivisions) {
577     if (SkCubicType::kLoop != fCurrCubicType) {
578         // Serpentines and cusps are always monotonic after chopping around inflection points.
579         SkASSERT(!SkCubicIsDegenerate(fCurrCubicType));
580 
581         if (AppendCubicMode::kApproximate == mode) {
582             // This section passes through an inflection point, so we can get away with a flat line.
583             // This can cause some curves to feel slightly more flat when inspected rigorously back
584             // and forth against another renderer, but for now this seems acceptable given the
585             // simplicity.
586             this->appendLine(p0, p3);
587             return;
588         }
589     } else {
590         Sk2f tan0, tan1;
591         get_cubic_tangents(p0, p1, p2, p3, &tan0, &tan1);
592 
593         if (maxSubdivisions && !is_convex_curve_monotonic(p0, tan0, p3, tan1)) {
594             this->chopAndAppendCubicAtMidTangent(mode, p0, p1, p2, p3, tan0, tan1,
595                                                  maxSubdivisions - 1);
596             return;
597         }
598 
599         if (AppendCubicMode::kApproximate == mode) {
600             Sk2f c;
601             if (!is_cubic_nearly_quadratic(p0, p1, p2, p3, tan0, tan1, &c) && maxSubdivisions) {
602                 this->chopAndAppendCubicAtMidTangent(mode, p0, p1, p2, p3, tan0, tan1,
603                                                      maxSubdivisions - 1);
604                 return;
605             }
606 
607             this->appendMonotonicQuadratic(p0, c, p3);
608             return;
609         }
610     }
611 
612     // Don't send curves to the GPU if we know they are nearly flat (or just very small).
613     // Since the cubic segment is known to be convex at this point, our flatness check is simple.
614     if (are_collinear(p0, (p1 + p2) * .5f, p3)) {
615         this->appendLine(p0, p3);
616         return;
617     }
618 
619     SkASSERT(fPoints.back() == SkPoint::Make(p0[0], p0[1]));
620     SkASSERT((p0 != p3).anyTrue());
621     p1.store(&fPoints.push_back());
622     p2.store(&fPoints.push_back());
623     p3.store(&fPoints.push_back());
624     fVerbs.push_back(Verb::kMonotonicCubicTo);
625     ++fCurrContourTallies.fCubics;
626 }
627 
628 // Given a convex curve segment with the following order-2 tangent function:
629 //
630 //                                                       |C2x  C2y|
631 //     tan = some_scale * |dx/dt  dy/dt| = |t^2  t  1| * |C1x  C1y|
632 //                                                       |C0x  C0y|
633 //
634 // This function finds the T value whose tangent angle is halfway between the tangents at T=0 and
635 // T=1 (tan0 and tan1).
find_midtangent(const Sk2f & tan0,const Sk2f & tan1,const Sk2f & C2,const Sk2f & C1,const Sk2f & C0)636 static inline float find_midtangent(const Sk2f& tan0, const Sk2f& tan1,
637                                     const Sk2f& C2, const Sk2f& C1, const Sk2f& C0) {
638     // Tangents point in the direction of increasing T, so tan0 and -tan1 both point toward the
639     // midtangent. 'n' will therefore bisect tan0 and -tan1, giving us the normal to the midtangent.
640     //
641     //     n dot midtangent = 0
642     //
643     Sk2f n = normalize(tan0) - normalize(tan1);
644 
645     // Find the T value at the midtangent. This is a simple quadratic equation:
646     //
647     //     midtangent dot n = 0
648     //
649     //     (|t^2  t  1| * C) dot n = 0
650     //
651     //     |t^2  t  1| dot C*n = 0
652     //
653     // First find coeffs = C*n.
654     Sk4f C[2];
655     Sk2f::Store4(C, C2, C1, C0, 0);
656     Sk4f coeffs = C[0]*n[0] + C[1]*n[1];
657 
658     // Now solve the quadratic.
659     float a = coeffs[0], b = coeffs[1], c = coeffs[2];
660     float discr = b*b - 4*a*c;
661     if (discr < 0) {
662         return 0; // This will only happen if the curve is a line.
663     }
664 
665     // The roots are q/a and c/q. Pick the one closer to T=.5.
666     float q = -.5f * (b + copysignf(std::sqrt(discr), b));
667     float r = .5f*q*a;
668     return std::abs(q*q - r) < std::abs(a*c - r) ? q/a : c/q;
669 }
670 
chopAndAppendCubicAtMidTangent(AppendCubicMode mode,const Sk2f & p0,const Sk2f & p1,const Sk2f & p2,const Sk2f & p3,const Sk2f & tan0,const Sk2f & tan1,int maxFutureSubdivisions)671 inline void GrCCFillGeometry::chopAndAppendCubicAtMidTangent(AppendCubicMode mode, const Sk2f& p0,
672                                                              const Sk2f& p1, const Sk2f& p2,
673                                                              const Sk2f& p3, const Sk2f& tan0,
674                                                              const Sk2f& tan1,
675                                                              int maxFutureSubdivisions) {
676     float midT = find_midtangent(tan0, tan1, p3 + (p1 - p2)*3 - p0,
677                                              (p0 - p1*2 + p2)*2,
678                                              p1 - p0);
679     // Use positive logic since NaN fails comparisons. (However midT should not be NaN since we cull
680     // near-flat cubics in cubicTo().)
681     if (!(midT > 0 && midT < 1)) {
682         // The cubic is flat. Otherwise there would be a real midtangent inside T=0..1.
683         this->appendLine(p0, p3);
684         return;
685     }
686 
687     Sk2f p01, p02, pT, p11, p12;
688     chop_cubic(p0, p1, p2, p3, midT, &p01, &p02, &pT, &p11, &p12);
689     this->appendCubics(mode, p0, p01, p02, pT, maxFutureSubdivisions);
690     this->appendCubics(mode, pT, p11, p12, p3, maxFutureSubdivisions);
691 }
692 
conicTo(const SkPoint P[3],float w)693 void GrCCFillGeometry::conicTo(const SkPoint P[3], float w) {
694     SkASSERT(fBuildingContour);
695     SkASSERT(P[0] == fPoints.back());
696     Sk2f p0 = Sk2f::Load(P);
697     Sk2f p1 = Sk2f::Load(P+1);
698     Sk2f p2 = Sk2f::Load(P+2);
699 
700     Sk2f tan0 = p1 - p0;
701     Sk2f tan1 = p2 - p1;
702 
703     if (!is_convex_curve_monotonic(p0, tan0, p2, tan1)) {
704         // The derivative of a conic has a cumbersome order-4 denominator. However, this isn't
705         // necessary if we are only interested in a vector in the same *direction* as a given
706         // tangent line. Since the denominator scales dx and dy uniformly, we can throw it out
707         // completely after evaluating the derivative with the standard quotient rule. This leaves
708         // us with a simpler quadratic function that we use to find the midtangent.
709         float midT = find_midtangent(tan0, tan1, (w - 1) * (p2 - p0),
710                                                  (p2 - p0) - 2*w*(p1 - p0),
711                                                  w*(p1 - p0));
712         // Use positive logic since NaN fails comparisons. (However midT should not be NaN since we
713         // cull near-linear conics above. And while w=0 is flat, it's not a line and has valid
714         // midtangents.)
715         if (!(midT > 0 && midT < 1)) {
716             // The conic is flat. Otherwise there would be a real midtangent inside T=0..1.
717             this->appendLine(p0, p2);
718             return;
719         }
720 
721         // Chop the conic at midtangent to produce two monotonic segments.
722         Sk4f p3d0 = Sk4f(p0[0], p0[1], 1, 0);
723         Sk4f p3d1 = Sk4f(p1[0], p1[1], 1, 0) * w;
724         Sk4f p3d2 = Sk4f(p2[0], p2[1], 1, 0);
725         Sk4f midT4 = midT;
726 
727         Sk4f p3d01 = lerp(p3d0, p3d1, midT4);
728         Sk4f p3d12 = lerp(p3d1, p3d2, midT4);
729         Sk4f p3d012 = lerp(p3d01, p3d12, midT4);
730 
731         Sk2f midpoint = Sk2f(p3d012[0], p3d012[1]) / p3d012[2];
732         Sk2f ww = Sk2f(p3d01[2], p3d12[2]) * Sk2f(p3d012[2]).rsqrt();
733 
734         this->appendMonotonicConic(p0, Sk2f(p3d01[0], p3d01[1]) / p3d01[2], midpoint, ww[0]);
735         this->appendMonotonicConic(midpoint, Sk2f(p3d12[0], p3d12[1]) / p3d12[2], p2, ww[1]);
736         return;
737     }
738 
739     this->appendMonotonicConic(p0, p1, p2, w);
740 }
741 
appendMonotonicConic(const Sk2f & p0,const Sk2f & p1,const Sk2f & p2,float w)742 void GrCCFillGeometry::appendMonotonicConic(const Sk2f& p0, const Sk2f& p1, const Sk2f& p2,
743                                             float w) {
744     SkASSERT(w >= 0);
745 
746     Sk2f base = p2 - p0;
747     Sk2f baseAbs = base.abs();
748     float baseWidth = baseAbs[0] + baseAbs[1];
749 
750     // Find the height of the curve. Max height always occurs at T=.5 for conics.
751     Sk2f d = (p1 - p0) * SkNx_shuffle<1,0>(base);
752     float h1 = std::abs(d[1] - d[0]); // Height of p1 above the base.
753     float ht = h1*w, hs = 1 + w; // Height of the conic = ht/hs.
754 
755     // i.e. (ht/hs <= baseWidth * kFlatnessThreshold). Use "<=" in case base == 0.
756     if (ht <= (baseWidth*hs) * kFlatnessThreshold) {
757         // We are flat. (See rationale in are_collinear.)
758         this->appendLine(p0, p2);
759         return;
760     }
761 
762     // i.e. (w > 1 && h1 - ht/hs < baseWidth).
763     if (w > 1 && h1*hs - ht < baseWidth*hs) {
764         // If we get within 1px of p1 when w > 1, we will pick up artifacts from the implicit
765         // function's reflection. Chop at max height (T=.5) and draw a triangle instead.
766         Sk2f p1w = p1*w;
767         Sk2f ab = p0 + p1w;
768         Sk2f bc = p1w + p2;
769         Sk2f highpoint = (ab + bc) / (2*(1 + w));
770         this->appendLine(p0, highpoint);
771         this->appendLine(highpoint, p2);
772         return;
773     }
774 
775     SkASSERT(fPoints.back() == SkPoint::Make(p0[0], p0[1]));
776     SkASSERT((p0 != p2).anyTrue());
777     p1.store(&fPoints.push_back());
778     p2.store(&fPoints.push_back());
779     fConicWeights.push_back(w);
780     fVerbs.push_back(Verb::kMonotonicConicTo);
781     ++fCurrContourTallies.fConics;
782 }
783 
endContour()784 GrCCFillGeometry::PrimitiveTallies GrCCFillGeometry::endContour() {
785     SkASSERT(fBuildingContour);
786     SkASSERT(fVerbs.count() >= fCurrContourTallies.fTriangles);
787 
788     // The fTriangles field currently contains this contour's starting verb index. We can now
789     // use it to calculate the size of the contour's fan.
790     int fanSize = fVerbs.count() - fCurrContourTallies.fTriangles;
791     if (fPoints.back() == fCurrAnchorPoint) {
792         --fanSize;
793         fVerbs.push_back(Verb::kEndClosedContour);
794     } else {
795         fVerbs.push_back(Verb::kEndOpenContour);
796     }
797 
798     fCurrContourTallies.fTriangles = SkTMax(fanSize - 2, 0);
799 
800     SkDEBUGCODE(fBuildingContour = false);
801     return fCurrContourTallies;
802 }
803