[go: nahoru, domu]
Google Git
Sign in
android / platform / external / chromium_org / third_party / skia / src / dd4bd43413220a5532d6a4c35882aced941bb2ee / . / core / SkGeometry.cpp
blob: a12f961948c77b6695ce2e250fa66d94e3eef96c [file] [log] [blame]
epoger@google.com685cfc02011-07-28 14:26:00 +0000[diff] [blame]1
2/*
3 * Copyright 2006 The Android Open Source Project
4 *
5 * Use of this source code is governed by a BSD-style license that can be
6 * found in the LICENSE file.
7 */
8
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]9
10#include "SkGeometry.h"
11#include "Sk64.h"
12#include "SkMatrix.h"
13
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]14bool SkXRayCrossesLine(const SkXRay& pt, const SkPoint pts[2], bool* ambiguous) {
15 if (ambiguous) {
16 *ambiguous = false;
17 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]18 // Determine quick discards.
19 // Consider query line going exactly through point 0 to not
20 // intersect, for symmetry with SkXRayCrossesMonotonicCubic.
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]21 if (pt.fY == pts[0].fY) {
22 if (ambiguous) {
23 *ambiguous = true;
24 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]25 return false;
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]26 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]27 if (pt.fY < pts[0].fY && pt.fY < pts[1].fY)
28 return false;
29 if (pt.fY > pts[0].fY && pt.fY > pts[1].fY)
30 return false;
31 if (pt.fX > pts[0].fX && pt.fX > pts[1].fX)
32 return false;
33 // Determine degenerate cases
34 if (SkScalarNearlyZero(pts[0].fY - pts[1].fY))
35 return false;
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]36 if (SkScalarNearlyZero(pts[0].fX - pts[1].fX)) {
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]37 // We've already determined the query point lies within the
38 // vertical range of the line segment.
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]39 if (pt.fX <= pts[0].fX) {
40 if (ambiguous) {
41 *ambiguous = (pt.fY == pts[1].fY);
42 }
43 return true;
44 }
45 return false;
46 }
47 // Ambiguity check
48 if (pt.fY == pts[1].fY) {
49 if (pt.fX <= pts[1].fX) {
50 if (ambiguous) {
51 *ambiguous = true;
52 }
53 return true;
54 }
55 return false;
56 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]57 // Full line segment evaluation
58 SkScalar delta_y = pts[1].fY - pts[0].fY;
59 SkScalar delta_x = pts[1].fX - pts[0].fX;
60 SkScalar slope = SkScalarDiv(delta_y, delta_x);
61 SkScalar b = pts[0].fY - SkScalarMul(slope, pts[0].fX);
62 // Solve for x coordinate at y = pt.fY
63 SkScalar x = SkScalarDiv(pt.fY - b, slope);
64 return pt.fX <= x;
65}
66
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]67/** If defined, this makes eval_quad and eval_cubic do more setup (sometimes
68 involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul.
69 May also introduce overflow of fixed when we compute our setup.
70*/
71#ifdef SK_SCALAR_IS_FIXED
72 #define DIRECT_EVAL_OF_POLYNOMIALS
73#endif
74
75////////////////////////////////////////////////////////////////////////
76
77#ifdef SK_SCALAR_IS_FIXED
78 static int is_not_monotonic(int a, int b, int c, int d)
79 {
80 return (((a - b) | (b - c) | (c - d)) & ((b - a) | (c - b) | (d - c))) >> 31;
81 }
82
83 static int is_not_monotonic(int a, int b, int c)
84 {
85 return (((a - b) | (b - c)) & ((b - a) | (c - b))) >> 31;
86 }
87#else
88 static int is_not_monotonic(float a, float b, float c)
89 {
90 float ab = a - b;
91 float bc = b - c;
92 if (ab < 0)
93 bc = -bc;
94 return ab == 0 || bc < 0;
95 }
96#endif
97
98////////////////////////////////////////////////////////////////////////
99
100static bool is_unit_interval(SkScalar x)
101{
102 return x > 0 && x < SK_Scalar1;
103}
104
105static int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio)
106{
107 SkASSERT(ratio);
108
109 if (numer < 0)
110 {
111 numer = -numer;
112 denom = -denom;
113 }
114
115 if (denom == 0 || numer == 0 || numer >= denom)
116 return 0;
117
118 SkScalar r = SkScalarDiv(numer, denom);
reed@android.com7c83e1c2010-03-08 17:44:42 +0000[diff] [blame]119 if (SkScalarIsNaN(r)) {
120 return 0;
121 }
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]122 SkASSERT(r >= 0 && r < SK_Scalar1);
123 if (r == 0) // catch underflow if numer <<<< denom
124 return 0;
125 *ratio = r;
126 return 1;
127}
128
129/** From Numerical Recipes in C.
130
131 Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C])
132 x1 = Q / A
133 x2 = C / Q
134*/
135int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2])
136{
137 SkASSERT(roots);
138
139 if (A == 0)
140 return valid_unit_divide(-C, B, roots);
141
142 SkScalar* r = roots;
143
144#ifdef SK_SCALAR_IS_FLOAT
145 float R = B*B - 4*A*C;
reed@android.com7c83e1c2010-03-08 17:44:42 +0000[diff] [blame]146 if (R < 0 || SkScalarIsNaN(R)) { // complex roots
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]147 return 0;
reed@android.com7c83e1c2010-03-08 17:44:42 +0000[diff] [blame]148 }
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]149 R = sk_float_sqrt(R);
150#else
151 Sk64 RR, tmp;
152
153 RR.setMul(B,B);
154 tmp.setMul(A,C);
155 tmp.shiftLeft(2);
156 RR.sub(tmp);
157 if (RR.isNeg())
158 return 0;
159 SkFixed R = RR.getSqrt();
160#endif
161
162 SkScalar Q = (B < 0) ? -(B-R)/2 : -(B+R)/2;
163 r += valid_unit_divide(Q, A, r);
164 r += valid_unit_divide(C, Q, r);
165 if (r - roots == 2)
166 {
167 if (roots[0] > roots[1])
168 SkTSwap<SkScalar>(roots[0], roots[1]);
169 else if (roots[0] == roots[1]) // nearly-equal?
170 r -= 1; // skip the double root
171 }
172 return (int)(r - roots);
173}
174
175#ifdef SK_SCALAR_IS_FIXED
176/** Trim A/B/C down so that they are all <= 32bits
177 and then call SkFindUnitQuadRoots()
178*/
179static int Sk64FindFixedQuadRoots(const Sk64& A, const Sk64& B, const Sk64& C, SkFixed roots[2])
180{
181 int na = A.shiftToMake32();
182 int nb = B.shiftToMake32();
183 int nc = C.shiftToMake32();
184
185 int shift = SkMax32(na, SkMax32(nb, nc));
186 SkASSERT(shift >= 0);
187
188 return SkFindUnitQuadRoots(A.getShiftRight(shift), B.getShiftRight(shift), C.getShiftRight(shift), roots);
189}
190#endif
191
192/////////////////////////////////////////////////////////////////////////////////////
193/////////////////////////////////////////////////////////////////////////////////////
194
195static SkScalar eval_quad(const SkScalar src[], SkScalar t)
196{
197 SkASSERT(src);
198 SkASSERT(t >= 0 && t <= SK_Scalar1);
199
200#ifdef DIRECT_EVAL_OF_POLYNOMIALS
201 SkScalar C = src[0];
202 SkScalar A = src[4] - 2 * src[2] + C;
203 SkScalar B = 2 * (src[2] - C);
204 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
205#else
206 SkScalar ab = SkScalarInterp(src[0], src[2], t);
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]207 SkScalar bc = SkScalarInterp(src[2], src[4], t);
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]208 return SkScalarInterp(ab, bc, t);
209#endif
210}
211
212static SkScalar eval_quad_derivative(const SkScalar src[], SkScalar t)
213{
214 SkScalar A = src[4] - 2 * src[2] + src[0];
215 SkScalar B = src[2] - src[0];
216
217 return 2 * SkScalarMulAdd(A, t, B);
218}
219
220static SkScalar eval_quad_derivative_at_half(const SkScalar src[])
221{
222 SkScalar A = src[4] - 2 * src[2] + src[0];
223 SkScalar B = src[2] - src[0];
224 return A + 2 * B;
225}
226
227void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt, SkVector* tangent)
228{
229 SkASSERT(src);
230 SkASSERT(t >= 0 && t <= SK_Scalar1);
231
232 if (pt)
233 pt->set(eval_quad(&src[0].fX, t), eval_quad(&src[0].fY, t));
234 if (tangent)
235 tangent->set(eval_quad_derivative(&src[0].fX, t),
236 eval_quad_derivative(&src[0].fY, t));
237}
238
239void SkEvalQuadAtHalf(const SkPoint src[3], SkPoint* pt, SkVector* tangent)
240{
241 SkASSERT(src);
242
243 if (pt)
244 {
245 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
246 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
247 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
248 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
249 pt->set(SkScalarAve(x01, x12), SkScalarAve(y01, y12));
250 }
251 if (tangent)
252 tangent->set(eval_quad_derivative_at_half(&src[0].fX),
253 eval_quad_derivative_at_half(&src[0].fY));
254}
255
256static void interp_quad_coords(const SkScalar* src, SkScalar* dst, SkScalar t)
257{
258 SkScalar ab = SkScalarInterp(src[0], src[2], t);
259 SkScalar bc = SkScalarInterp(src[2], src[4], t);
260
261 dst[0] = src[0];
262 dst[2] = ab;
263 dst[4] = SkScalarInterp(ab, bc, t);
264 dst[6] = bc;
265 dst[8] = src[4];
266}
267
268void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t)
269{
270 SkASSERT(t > 0 && t < SK_Scalar1);
271
272 interp_quad_coords(&src[0].fX, &dst[0].fX, t);
273 interp_quad_coords(&src[0].fY, &dst[0].fY, t);
274}
275
276void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5])
277{
278 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
279 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
280 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
281 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
282
283 dst[0] = src[0];
284 dst[1].set(x01, y01);
285 dst[2].set(SkScalarAve(x01, x12), SkScalarAve(y01, y12));
286 dst[3].set(x12, y12);
287 dst[4] = src[2];
288}
289
290/** Quad'(t) = At + B, where
291 A = 2(a - 2b + c)
292 B = 2(b - a)
293 Solve for t, only if it fits between 0 < t < 1
294*/
295int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1])
296{
297 /* At + B == 0
298 t = -B / A
299 */
300#ifdef SK_SCALAR_IS_FIXED
301 return is_not_monotonic(a, b, c) && valid_unit_divide(a - b, a - b - b + c, tValue);
302#else
303 return valid_unit_divide(a - b, a - b - b + c, tValue);
304#endif
305}
306
reed@android.come5dd6cd2009-01-15 14:38:33 +0000[diff] [blame]307static inline void flatten_double_quad_extrema(SkScalar coords[14])
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]308{
309 coords[2] = coords[6] = coords[4];
310}
311
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]312/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is
reed@android.com001bd972009-11-17 18:47:52 +0000[diff] [blame]313 stored in dst[]. Guarantees that the 1/2 quads will be monotonic.
314 */
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]315int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5])
316{
317 SkASSERT(src);
318 SkASSERT(dst);
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]319
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]320#if 0
321 static bool once = true;
322 if (once)
323 {
324 once = false;
325 SkPoint s[3] = { 0, 26398, 0, 26331, 0, 20621428 };
326 SkPoint d[6];
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]327
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]328 int n = SkChopQuadAtYExtrema(s, d);
329 SkDebugf("chop=%d, Y=[%x %x %x %x %x %x]\n", n, d[0].fY, d[1].fY, d[2].fY, d[3].fY, d[4].fY, d[5].fY);
330 }
331#endif
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]332
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]333 SkScalar a = src[0].fY;
334 SkScalar b = src[1].fY;
335 SkScalar c = src[2].fY;
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]336
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]337 if (is_not_monotonic(a, b, c))
338 {
339 SkScalar tValue;
340 if (valid_unit_divide(a - b, a - b - b + c, &tValue))
341 {
342 SkChopQuadAt(src, dst, tValue);
343 flatten_double_quad_extrema(&dst[0].fY);
344 return 1;
345 }
346 // if we get here, we need to force dst to be monotonic, even though
347 // we couldn't compute a unit_divide value (probably underflow).
348 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c;
349 }
350 dst[0].set(src[0].fX, a);
351 dst[1].set(src[1].fX, b);
352 dst[2].set(src[2].fX, c);
353 return 0;
354}
355
reed@android.com001bd972009-11-17 18:47:52 +0000[diff] [blame]356/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is
357 stored in dst[]. Guarantees that the 1/2 quads will be monotonic.
358 */
359int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5])
360{
361 SkASSERT(src);
362 SkASSERT(dst);
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]363
reed@android.com001bd972009-11-17 18:47:52 +0000[diff] [blame]364 SkScalar a = src[0].fX;
365 SkScalar b = src[1].fX;
366 SkScalar c = src[2].fX;
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]367
reed@android.com001bd972009-11-17 18:47:52 +0000[diff] [blame]368 if (is_not_monotonic(a, b, c)) {
369 SkScalar tValue;
370 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) {
371 SkChopQuadAt(src, dst, tValue);
372 flatten_double_quad_extrema(&dst[0].fX);
373 return 1;
374 }
375 // if we get here, we need to force dst to be monotonic, even though
376 // we couldn't compute a unit_divide value (probably underflow).
377 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c;
378 }
379 dst[0].set(a, src[0].fY);
380 dst[1].set(b, src[1].fY);
381 dst[2].set(c, src[2].fY);
382 return 0;
383}
384
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]385// F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2
386// F'(t) = 2 (b - a) + 2 (a - 2b + c) t
387// F''(t) = 2 (a - 2b + c)
388//
389// A = 2 (b - a)
390// B = 2 (a - 2b + c)
391//
392// Maximum curvature for a quadratic means solving
393// Fx' Fx'' + Fy' Fy'' = 0
394//
395// t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2)
396//
397int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5])
398{
399 SkScalar Ax = src[1].fX - src[0].fX;
400 SkScalar Ay = src[1].fY - src[0].fY;
401 SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX;
402 SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY;
403 SkScalar t = 0; // 0 means don't chop
404
405#ifdef SK_SCALAR_IS_FLOAT
406 (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t);
407#else
408 // !!! should I use SkFloat here? seems like it
409 Sk64 numer, denom, tmp;
410
411 numer.setMul(Ax, -Bx);
412 tmp.setMul(Ay, -By);
413 numer.add(tmp);
414
415 if (numer.isPos()) // do nothing if numer <= 0
416 {
417 denom.setMul(Bx, Bx);
418 tmp.setMul(By, By);
419 denom.add(tmp);
420 SkASSERT(!denom.isNeg());
421 if (numer < denom)
422 {
423 t = numer.getFixedDiv(denom);
424 SkASSERT(t >= 0 && t <= SK_Fixed1); // assert that we're numerically stable (ha!)
425 if ((unsigned)t >= SK_Fixed1) // runtime check for numerical stability
426 t = 0; // ignore the chop
427 }
428 }
429#endif
430
431 if (t == 0)
432 {
433 memcpy(dst, src, 3 * sizeof(SkPoint));
434 return 1;
435 }
436 else
437 {
438 SkChopQuadAt(src, dst, t);
439 return 2;
440 }
441}
442
reed@google.com007593e2011-07-27 13:54:36 +0000[diff] [blame]443#ifdef SK_SCALAR_IS_FLOAT
444 #define SK_ScalarTwoThirds (0.666666666f)
445#else
446 #define SK_ScalarTwoThirds ((SkFixed)(43691))
447#endif
448
449void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) {
450 const SkScalar scale = SK_ScalarTwoThirds;
451 dst[0] = src[0];
452 dst[1].set(src[0].fX + SkScalarMul(src[1].fX - src[0].fX, scale),
453 src[0].fY + SkScalarMul(src[1].fY - src[0].fY, scale));
454 dst[2].set(src[2].fX + SkScalarMul(src[1].fX - src[2].fX, scale),
455 src[2].fY + SkScalarMul(src[1].fY - src[2].fY, scale));
456 dst[3] = src[2];
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]457}
458
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]459////////////////////////////////////////////////////////////////////////////////////////
460///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS /////
461////////////////////////////////////////////////////////////////////////////////////////
462
463static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4])
464{
465 coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0];
466 coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]);
467 coeff[2] = 3*(pt[2] - pt[0]);
468 coeff[3] = pt[0];
469}
470
471void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4])
472{
473 SkASSERT(pts);
474
475 if (cx)
476 get_cubic_coeff(&pts[0].fX, cx);
477 if (cy)
478 get_cubic_coeff(&pts[0].fY, cy);
479}
480
481static SkScalar eval_cubic(const SkScalar src[], SkScalar t)
482{
483 SkASSERT(src);
484 SkASSERT(t >= 0 && t <= SK_Scalar1);
485
486 if (t == 0)
487 return src[0];
488
489#ifdef DIRECT_EVAL_OF_POLYNOMIALS
490 SkScalar D = src[0];
491 SkScalar A = src[6] + 3*(src[2] - src[4]) - D;
492 SkScalar B = 3*(src[4] - src[2] - src[2] + D);
493 SkScalar C = 3*(src[2] - D);
494
495 return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D);
496#else
497 SkScalar ab = SkScalarInterp(src[0], src[2], t);
498 SkScalar bc = SkScalarInterp(src[2], src[4], t);
499 SkScalar cd = SkScalarInterp(src[4], src[6], t);
500 SkScalar abc = SkScalarInterp(ab, bc, t);
501 SkScalar bcd = SkScalarInterp(bc, cd, t);
502 return SkScalarInterp(abc, bcd, t);
503#endif
504}
505
506/** return At^2 + Bt + C
507*/
508static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t)
509{
510 SkASSERT(t >= 0 && t <= SK_Scalar1);
511
512 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
513}
514
515static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t)
516{
517 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0];
518 SkScalar B = 2*(src[4] - 2 * src[2] + src[0]);
519 SkScalar C = src[2] - src[0];
520
521 return eval_quadratic(A, B, C, t);
522}
523
524static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t)
525{
526 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0];
527 SkScalar B = src[4] - 2 * src[2] + src[0];
528
529 return SkScalarMulAdd(A, t, B);
530}
531
532void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature)
533{
534 SkASSERT(src);
535 SkASSERT(t >= 0 && t <= SK_Scalar1);
536
537 if (loc)
538 loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t));
539 if (tangent)
540 tangent->set(eval_cubic_derivative(&src[0].fX, t),
541 eval_cubic_derivative(&src[0].fY, t));
542 if (curvature)
543 curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t),
544 eval_cubic_2ndDerivative(&src[0].fY, t));
545}
546
547/** Cubic'(t) = At^2 + Bt + C, where
548 A = 3(-a + 3(b - c) + d)
549 B = 6(a - 2b + c)
550 C = 3(b - a)
551 Solve for t, keeping only those that fit betwee 0 < t < 1
552*/
553int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2])
554{
555#ifdef SK_SCALAR_IS_FIXED
556 if (!is_not_monotonic(a, b, c, d))
557 return 0;
558#endif
559
560 // we divide A,B,C by 3 to simplify
561 SkScalar A = d - a + 3*(b - c);
562 SkScalar B = 2*(a - b - b + c);
563 SkScalar C = b - a;
564
565 return SkFindUnitQuadRoots(A, B, C, tValues);
566}
567
568static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t)
569{
570 SkScalar ab = SkScalarInterp(src[0], src[2], t);
571 SkScalar bc = SkScalarInterp(src[2], src[4], t);
572 SkScalar cd = SkScalarInterp(src[4], src[6], t);
573 SkScalar abc = SkScalarInterp(ab, bc, t);
574 SkScalar bcd = SkScalarInterp(bc, cd, t);
575 SkScalar abcd = SkScalarInterp(abc, bcd, t);
576
577 dst[0] = src[0];
578 dst[2] = ab;
579 dst[4] = abc;
580 dst[6] = abcd;
581 dst[8] = bcd;
582 dst[10] = cd;
583 dst[12] = src[6];
584}
585
586void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t)
587{
588 SkASSERT(t > 0 && t < SK_Scalar1);
589
590 interp_cubic_coords(&src[0].fX, &dst[0].fX, t);
591 interp_cubic_coords(&src[0].fY, &dst[0].fY, t);
592}
593
reed@android.com17bdc092009-08-28 20:06:54 +0000[diff] [blame]594/* http://code.google.com/p/skia/issues/detail?id=32
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]595
reed@android.com17bdc092009-08-28 20:06:54 +0000[diff] [blame]596 This test code would fail when we didn't check the return result of
597 valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is
598 that after the first chop, the parameters to valid_unit_divide are equal
599 (thanks to finite float precision and rounding in the subtracts). Thus
600 even though the 2nd tValue looks < 1.0, after we renormalize it, we end
601 up with 1.0, hence the need to check and just return the last cubic as
602 a degenerate clump of 4 points in the sampe place.
603
604 static void test_cubic() {
605 SkPoint src[4] = {
606 { 556.25000, 523.03003 },
607 { 556.23999, 522.96002 },
608 { 556.21997, 522.89001 },
609 { 556.21997, 522.82001 }
610 };
611 SkPoint dst[10];
612 SkScalar tval[] = { 0.33333334f, 0.99999994f };
613 SkChopCubicAt(src, dst, tval, 2);
614 }
615 */
616
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]617void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots)
618{
619#ifdef SK_DEBUG
620 {
621 for (int i = 0; i < roots - 1; i++)
622 {
623 SkASSERT(is_unit_interval(tValues[i]));
624 SkASSERT(is_unit_interval(tValues[i+1]));
625 SkASSERT(tValues[i] < tValues[i+1]);
626 }
627 }
628#endif
629
630 if (dst)
631 {
632 if (roots == 0) // nothing to chop
633 memcpy(dst, src, 4*sizeof(SkPoint));
634 else
635 {
636 SkScalar t = tValues[0];
637 SkPoint tmp[4];
638
639 for (int i = 0; i < roots; i++)
640 {
641 SkChopCubicAt(src, dst, t);
642 if (i == roots - 1)
643 break;
644
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]645 dst += 3;
reed@android.com17bdc092009-08-28 20:06:54 +0000[diff] [blame]646 // have src point to the remaining cubic (after the chop)
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]647 memcpy(tmp, dst, 4 * sizeof(SkPoint));
648 src = tmp;
reed@android.com17bdc092009-08-28 20:06:54 +0000[diff] [blame]649
650 // watch out in case the renormalized t isn't in range
651 if (!valid_unit_divide(tValues[i+1] - tValues[i],
652 SK_Scalar1 - tValues[i], &t)) {
653 // if we can't, just create a degenerate cubic
654 dst[4] = dst[5] = dst[6] = src[3];
655 break;
656 }
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]657 }
658 }
659 }
660}
661
662void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7])
663{
664 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
665 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
666 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
667 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
668 SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX);
669 SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY);
670
671 SkScalar x012 = SkScalarAve(x01, x12);
672 SkScalar y012 = SkScalarAve(y01, y12);
673 SkScalar x123 = SkScalarAve(x12, x23);
674 SkScalar y123 = SkScalarAve(y12, y23);
675
676 dst[0] = src[0];
677 dst[1].set(x01, y01);
678 dst[2].set(x012, y012);
679 dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123));
680 dst[4].set(x123, y123);
681 dst[5].set(x23, y23);
682 dst[6] = src[3];
683}
684
685static void flatten_double_cubic_extrema(SkScalar coords[14])
686{
687 coords[4] = coords[8] = coords[6];
688}
689
690/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that
691 the resulting beziers are monotonic in Y. This is called by the scan converter.
692 Depending on what is returned, dst[] is treated as follows
693 0 dst[0..3] is the original cubic
694 1 dst[0..3] and dst[3..6] are the two new cubics
695 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics
696 If dst == null, it is ignored and only the count is returned.
697*/
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame]698int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) {
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]699 SkScalar tValues[2];
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame]700 int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY,
701 src[3].fY, tValues);
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]702
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]703 SkChopCubicAt(src, dst, tValues, roots);
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame]704 if (dst && roots > 0) {
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]705 // we do some cleanup to ensure our Y extrema are flat
706 flatten_double_cubic_extrema(&dst[0].fY);
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame]707 if (roots == 2) {
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]708 flatten_double_cubic_extrema(&dst[3].fY);
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame]709 }
710 }
711 return roots;
712}
713
714int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) {
715 SkScalar tValues[2];
716 int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX,
717 src[3].fX, tValues);
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]718
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame]719 SkChopCubicAt(src, dst, tValues, roots);
720 if (dst && roots > 0) {
721 // we do some cleanup to ensure our Y extrema are flat
722 flatten_double_cubic_extrema(&dst[0].fX);
723 if (roots == 2) {
724 flatten_double_cubic_extrema(&dst[3].fX);
725 }
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]726 }
727 return roots;
728}
729
730/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html
731
732 Inflection means that curvature is zero.
733 Curvature is [F' x F''] / [F'^3]
734 So we solve F'x X F''y - F'y X F''y == 0
735 After some canceling of the cubic term, we get
736 A = b - a
737 B = c - 2b + a
738 C = d - 3c + 3b - a
739 (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0
740*/
741int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[])
742{
743 SkScalar Ax = src[1].fX - src[0].fX;
744 SkScalar Ay = src[1].fY - src[0].fY;
745 SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX;
746 SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY;
747 SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX;
748 SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY;
749 int count;
750
751#ifdef SK_SCALAR_IS_FLOAT
752 count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues);
753#else
754 Sk64 A, B, C, tmp;
755
756 A.setMul(Bx, Cy);
757 tmp.setMul(By, Cx);
758 A.sub(tmp);
759
760 B.setMul(Ax, Cy);
761 tmp.setMul(Ay, Cx);
762 B.sub(tmp);
763
764 C.setMul(Ax, By);
765 tmp.setMul(Ay, Bx);
766 C.sub(tmp);
767
768 count = Sk64FindFixedQuadRoots(A, B, C, tValues);
769#endif
770
771 return count;
772}
773
774int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10])
775{
776 SkScalar tValues[2];
777 int count = SkFindCubicInflections(src, tValues);
778
779 if (dst)
780 {
781 if (count == 0)
782 memcpy(dst, src, 4 * sizeof(SkPoint));
783 else
784 SkChopCubicAt(src, dst, tValues, count);
785 }
786 return count + 1;
787}
788
789template <typename T> void bubble_sort(T array[], int count)
790{
791 for (int i = count - 1; i > 0; --i)
792 for (int j = i; j > 0; --j)
793 if (array[j] < array[j-1])
794 {
795 T tmp(array[j]);
796 array[j] = array[j-1];
797 array[j-1] = tmp;
798 }
799}
800
801#include "SkFP.h"
802
803// newton refinement
804#if 0
805static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root)
806{
807 // x1 = x0 - f(t) / f'(t)
808
809 SkFP T = SkScalarToFloat(root);
810 SkFP N, D;
811
812 // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2]
813 D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3);
814 D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2));
815 D = SkFPAdd(D, coeff[2]);
816
817 if (D == 0)
818 return root;
819
820 // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3]
821 N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]);
822 N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1]));
823 N = SkFPAdd(N, SkFPMul(T, coeff[2]));
824 N = SkFPAdd(N, coeff[3]);
825
826 if (N)
827 {
828 SkScalar delta = SkFPToScalar(SkFPDiv(N, D));
829
830 if (delta)
831 root -= delta;
832 }
833 return root;
834}
835#endif
836
reed@google.com7de50932012-02-29 20:59:24 +0000[diff] [blame]837/**
838 * Given an array and count, remove all pair-wise duplicates from the array,
839 * keeping the existing sorting, and return the new count
840 */
841static int collaps_duplicates(float array[], int count) {
reed@google.com7de50932012-02-29 20:59:24 +0000[diff] [blame]842 for (int n = count; n > 1; --n) {
843 if (array[0] == array[1]) {
844 for (int i = 1; i < n; ++i) {
845 array[i - 1] = array[i];
846 }
847 count -= 1;
848 } else {
849 array += 1;
850 }
851 }
852 return count;
853}
854
855#ifdef SK_DEBUG
856
857#define TEST_COLLAPS_ENTRY(array) array, SK_ARRAY_COUNT(array)
858
859static void test_collaps_duplicates() {
860 static bool gOnce;
861 if (gOnce) { return; }
862 gOnce = true;
863 const float src0[] = { 0 };
864 const float src1[] = { 0, 0 };
865 const float src2[] = { 0, 1 };
866 const float src3[] = { 0, 0, 0 };
867 const float src4[] = { 0, 0, 1 };
868 const float src5[] = { 0, 1, 1 };
869 const float src6[] = { 0, 1, 2 };
870 const struct {
871 const float* fData;
872 int fCount;
873 int fCollapsedCount;
874 } data[] = {
875 { TEST_COLLAPS_ENTRY(src0), 1 },
876 { TEST_COLLAPS_ENTRY(src1), 1 },
877 { TEST_COLLAPS_ENTRY(src2), 2 },
878 { TEST_COLLAPS_ENTRY(src3), 1 },
879 { TEST_COLLAPS_ENTRY(src4), 2 },
880 { TEST_COLLAPS_ENTRY(src5), 2 },
881 { TEST_COLLAPS_ENTRY(src6), 3 },
882 };
883 for (size_t i = 0; i < SK_ARRAY_COUNT(data); ++i) {
884 float dst[3];
885 memcpy(dst, data[i].fData, data[i].fCount * sizeof(dst[0]));
886 int count = collaps_duplicates(dst, data[i].fCount);
887 SkASSERT(data[i].fCollapsedCount == count);
888 for (int j = 1; j < count; ++j) {
889 SkASSERT(dst[j-1] < dst[j]);
890 }
891 }
892}
893#endif
894
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]895#if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop
896#pragma warning ( disable : 4702 )
897#endif
898
899/* Solve coeff(t) == 0, returning the number of roots that
900 lie withing 0 < t < 1.
901 coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3]
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]902
reed@google.com7de50932012-02-29 20:59:24 +0000[diff] [blame]903 Eliminates repeated roots (so that all tValues are distinct, and are always
904 in increasing order.
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]905*/
906static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3])
907{
908#ifndef SK_SCALAR_IS_FLOAT
909 return 0; // this is not yet implemented for software float
910#endif
911
912 if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic
913 {
914 return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues);
915 }
916
917 SkFP a, b, c, Q, R;
918
919 {
920 SkASSERT(coeff[0] != 0);
921
922 SkFP inva = SkFPInvert(coeff[0]);
923 a = SkFPMul(coeff[1], inva);
924 b = SkFPMul(coeff[2], inva);
925 c = SkFPMul(coeff[3], inva);
926 }
927 Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9);
928// R = (2*a*a*a - 9*a*b + 27*c) / 54;
929 R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2);
930 R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9));
931 R = SkFPAdd(R, SkFPMulInt(c, 27));
932 R = SkFPDivInt(R, 54);
933
934 SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q);
935 SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3);
936 SkFP adiv3 = SkFPDivInt(a, 3);
937
938 SkScalar* roots = tValues;
939 SkScalar r;
940
941 if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots
942 {
943#ifdef SK_SCALAR_IS_FLOAT
944 float theta = sk_float_acos(R / sk_float_sqrt(Q3));
945 float neg2RootQ = -2 * sk_float_sqrt(Q);
946
947 r = neg2RootQ * sk_float_cos(theta/3) - adiv3;
948 if (is_unit_interval(r))
949 *roots++ = r;
950
951 r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3;
952 if (is_unit_interval(r))
953 *roots++ = r;
954
955 r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3;
956 if (is_unit_interval(r))
957 *roots++ = r;
958
reed@google.com7de50932012-02-29 20:59:24 +0000[diff] [blame]959 SkDEBUGCODE(test_collaps_duplicates();)
960
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]961 // now sort the roots
reed@google.com7de50932012-02-29 20:59:24 +0000[diff] [blame]962 int count = (int)(roots - tValues);
963 SkASSERT((unsigned)count <= 3);
964 bubble_sort(tValues, count);
965 count = collaps_duplicates(tValues, count);
966 roots = tValues + count; // so we compute the proper count below
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]967#endif
968 }
969 else // we have 1 real root
970 {
971 SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3));
972 A = SkFPCubeRoot(A);
973 if (SkFPGT(R, 0))
974 A = SkFPNeg(A);
975
976 if (A != 0)
977 A = SkFPAdd(A, SkFPDiv(Q, A));
978 r = SkFPToScalar(SkFPSub(A, adiv3));
979 if (is_unit_interval(r))
980 *roots++ = r;
981 }
982
983 return (int)(roots - tValues);
984}
985
986/* Looking for F' dot F'' == 0
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]987
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]988 A = b - a
989 B = c - 2b + a
990 C = d - 3c + 3b - a
991
992 F' = 3Ct^2 + 6Bt + 3A
993 F'' = 6Ct + 6B
994
995 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
996*/
997static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4])
998{
999 SkScalar a = src[2] - src[0];
1000 SkScalar b = src[4] - 2 * src[2] + src[0];
1001 SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0];
1002
1003 SkFP A = SkScalarToFP(a);
1004 SkFP B = SkScalarToFP(b);
1005 SkFP C = SkScalarToFP(c);
1006
1007 coeff[0] = SkFPMul(C, C);
1008 coeff[1] = SkFPMulInt(SkFPMul(B, C), 3);
1009 coeff[2] = SkFPMulInt(SkFPMul(B, B), 2);
1010 coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A));
1011 coeff[3] = SkFPMul(A, B);
1012}
1013
1014// EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1
1015//#define kMinTValueForChopping (SK_Scalar1 / 256)
1016#define kMinTValueForChopping 0
1017
1018/* Looking for F' dot F'' == 0
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]1019
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1020 A = b - a
1021 B = c - 2b + a
1022 C = d - 3c + 3b - a
1023
1024 F' = 3Ct^2 + 6Bt + 3A
1025 F'' = 6Ct + 6B
1026
1027 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
1028*/
1029int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3])
1030{
1031 SkFP coeffX[4], coeffY[4];
1032 int i;
1033
1034 formulate_F1DotF2(&src[0].fX, coeffX);
1035 formulate_F1DotF2(&src[0].fY, coeffY);
1036
1037 for (i = 0; i < 4; i++)
1038 coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]);
1039
1040 SkScalar t[3];
1041 int count = solve_cubic_polynomial(coeffX, t);
1042 int maxCount = 0;
1043
1044 // now remove extrema where the curvature is zero (mins)
1045 // !!!! need a test for this !!!!
1046 for (i = 0; i < count; i++)
1047 {
1048 // if (not_min_curvature())
1049 if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping)
1050 tValues[maxCount++] = t[i];
1051 }
1052 return maxCount;
1053}
1054
1055int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3])
1056{
1057 SkScalar t_storage[3];
1058
1059 if (tValues == NULL)
1060 tValues = t_storage;
1061
1062 int count = SkFindCubicMaxCurvature(src, tValues);
1063
1064 if (dst)
1065 {
1066 if (count == 0)
1067 memcpy(dst, src, 4 * sizeof(SkPoint));
1068 else
1069 SkChopCubicAt(src, dst, tValues, count);
1070 }
1071 return count + 1;
1072}
1073
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1074bool SkXRayCrossesMonotonicCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) {
1075 if (ambiguous) {
1076 *ambiguous = false;
1077 }
1078
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1079 // Find the minimum and maximum y of the extrema, which are the
1080 // first and last points since this cubic is monotonic
1081 SkScalar min_y = SkMinScalar(cubic[0].fY, cubic[3].fY);
1082 SkScalar max_y = SkMaxScalar(cubic[0].fY, cubic[3].fY);
1083
1084 if (pt.fY == cubic[0].fY
1085 || pt.fY < min_y
1086 || pt.fY > max_y) {
1087 // The query line definitely does not cross the curve
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1088 if (ambiguous) {
1089 *ambiguous = (pt.fY == cubic[0].fY);
1090 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1091 return false;
1092 }
1093
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1094 bool pt_at_extremum = (pt.fY == cubic[3].fY);
1095
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1096 SkScalar min_x =
1097 SkMinScalar(
1098 SkMinScalar(
1099 SkMinScalar(cubic[0].fX, cubic[1].fX),
1100 cubic[2].fX),
1101 cubic[3].fX);
1102 if (pt.fX < min_x) {
1103 // The query line definitely crosses the curve
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1104 if (ambiguous) {
1105 *ambiguous = pt_at_extremum;
1106 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1107 return true;
1108 }
1109
1110 SkScalar max_x =
1111 SkMaxScalar(
1112 SkMaxScalar(
1113 SkMaxScalar(cubic[0].fX, cubic[1].fX),
1114 cubic[2].fX),
1115 cubic[3].fX);
1116 if (pt.fX > max_x) {
1117 // The query line definitely does not cross the curve
1118 return false;
1119 }
1120
1121 // Do a binary search to find the parameter value which makes y as
1122 // close as possible to the query point. See whether the query
1123 // line's origin is to the left of the associated x coordinate.
1124
1125 // kMaxIter is chosen as the number of mantissa bits for a float,
1126 // since there's no way we are going to get more precision by
1127 // iterating more times than that.
1128 const int kMaxIter = 23;
1129 SkPoint eval;
1130 int iter = 0;
1131 SkScalar upper_t;
1132 SkScalar lower_t;
1133 // Need to invert direction of t parameter if cubic goes up
1134 // instead of down
1135 if (cubic[3].fY > cubic[0].fY) {
1136 upper_t = SK_Scalar1;
1137 lower_t = SkFloatToScalar(0);
1138 } else {
1139 upper_t = SkFloatToScalar(0);
1140 lower_t = SK_Scalar1;
1141 }
1142 do {
1143 SkScalar t = SkScalarAve(upper_t, lower_t);
1144 SkEvalCubicAt(cubic, t, &eval, NULL, NULL);
1145 if (pt.fY > eval.fY) {
1146 lower_t = t;
1147 } else {
1148 upper_t = t;
1149 }
1150 } while (++iter < kMaxIter
1151 && !SkScalarNearlyZero(eval.fY - pt.fY));
1152 if (pt.fX <= eval.fX) {
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1153 if (ambiguous) {
1154 *ambiguous = pt_at_extremum;
1155 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1156 return true;
1157 }
1158 return false;
1159}
1160
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1161int SkNumXRayCrossingsForCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) {
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1162 int num_crossings = 0;
1163 SkPoint monotonic_cubics[10];
1164 int num_monotonic_cubics = SkChopCubicAtYExtrema(cubic, monotonic_cubics);
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1165 if (ambiguous) {
1166 *ambiguous = false;
1167 }
1168 bool locally_ambiguous;
1169 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[0], &locally_ambiguous))
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1170 ++num_crossings;
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1171 if (ambiguous) {
1172 *ambiguous |= locally_ambiguous;
1173 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1174 if (num_monotonic_cubics > 0)
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1175 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[3], &locally_ambiguous))
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1176 ++num_crossings;
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1177 if (ambiguous) {
1178 *ambiguous |= locally_ambiguous;
1179 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1180 if (num_monotonic_cubics > 1)
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1181 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[6], &locally_ambiguous))
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1182 ++num_crossings;
kbr@chromium.orgc1b53332010-07-07 22:20:35 +0000[diff] [blame]1183 if (ambiguous) {
1184 *ambiguous |= locally_ambiguous;
1185 }
reed@android.com5b4541e2010-02-05 20:41:02 +0000[diff] [blame]1186 return num_crossings;
1187}
1188
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1189////////////////////////////////////////////////////////////////////////////////
1190
1191/* Find t value for quadratic [a, b, c] = d.
1192 Return 0 if there is no solution within [0, 1)
1193*/
1194static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d)
1195{
1196 // At^2 + Bt + C = d
1197 SkScalar A = a - 2 * b + c;
1198 SkScalar B = 2 * (b - a);
1199 SkScalar C = a - d;
1200
1201 SkScalar roots[2];
1202 int count = SkFindUnitQuadRoots(A, B, C, roots);
1203
1204 SkASSERT(count <= 1);
1205 return count == 1 ? roots[0] : 0;
1206}
1207
robertphillips@google.comdd4bd432013-07-09 15:03:59 +0000[diff] [blame^]1208/* given a quad-curve and a point (x,y), chop the quad at that point and place
1209 the new off-curve point and endpoint into 'dest'. The new end point is used
1210 (rather than (x,y)) to compensate for numerical inaccuracies.
1211 Should only return false if the computed pos is the start of the curve
1212 (i.e. root == 0)
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1213*/
robertphillips@google.comdd4bd432013-07-09 15:03:59 +0000[diff] [blame^]1214static bool truncate_last_curve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* dest)
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1215{
1216 const SkScalar* base;
1217 SkScalar value;
1218
1219 if (SkScalarAbs(x) < SkScalarAbs(y)) {
1220 base = &quad[0].fX;
1221 value = x;
1222 } else {
1223 base = &quad[0].fY;
1224 value = y;
1225 }
1226
1227 // note: this returns 0 if it thinks value is out of range, meaning the
1228 // root might return something outside of [0, 1)
1229 SkScalar t = quad_solve(base[0], base[2], base[4], value);
1230
1231 if (t > 0)
1232 {
1233 SkPoint tmp[5];
1234 SkChopQuadAt(quad, tmp, t);
robertphillips@google.comdd4bd432013-07-09 15:03:59 +0000[diff] [blame^]1235 dest[0] = tmp[1];
1236 dest[1] = tmp[2];
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1237 return true;
1238 } else {
1239 /* t == 0 means either the value triggered a root outside of [0, 1)
1240 For our purposes, we can ignore the <= 0 roots, but we want to
1241 catch the >= 1 roots (which given our caller, will basically mean
1242 a root of 1, give-or-take numerical instability). If we are in the
1243 >= 1 case, return the existing offCurve point.
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]1244
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1245 The test below checks to see if we are close to the "end" of the
1246 curve (near base[4]). Rather than specifying a tolerance, I just
1247 check to see if value is on to the right/left of the middle point
1248 (depending on the direction/sign of the end points).
1249 */
1250 if ((base[0] < base[4] && value > base[2]) ||
1251 (base[0] > base[4] && value < base[2])) // should root have been 1
1252 {
robertphillips@google.comdd4bd432013-07-09 15:03:59 +0000[diff] [blame^]1253 dest[0] = quad[1];
1254 dest[1].set(x, y);
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1255 return true;
1256 }
1257 }
1258 return false;
1259}
1260
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1261#ifdef SK_SCALAR_IS_FLOAT
reed@google.com41bea9b2013-02-22 14:19:58 +0000[diff] [blame]1262
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1263// Due to floating point issues (i.e., 1.0f - SK_ScalarRoot2Over2 !=
skia.committer@gmail.com6695aed2012-10-19 02:01:19 +0000[diff] [blame]1264// SK_ScalarRoot2Over2 - SK_ScalarTanPIOver8) a cruder root2over2
1265// approximation is required to make the quad circle points convex. The
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1266// root of the problem is that with the root2over2 value in SkScalar.h
1267// the arcs really are ever so slightly concave. Some alternative fixes
1268// to this problem (besides just arbitrarily pushing out the mid-point as
1269// is done here) are:
1270// Adjust all the points (not just the middle) to both better approximate
1271// the curve and remain convex
1272// Switch over to using cubics rather then quads
1273// Use a different method to create the mid-point (e.g., compute
1274// the two side points, average them, then move it out as needed
reed@google.com41bea9b2013-02-22 14:19:58 +0000[diff] [blame]1275#define SK_ScalarRoot2Over2_QuadCircle 0.7072f
robertphillips@google.comc03af442012-10-19 01:26:18 +0000[diff] [blame]1276
1277#else
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1278 #define SK_ScalarRoot2Over2_QuadCircle SK_FixedRoot2Over2
1279#endif
1280
1281
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1282static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = {
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1283 { SK_Scalar1, 0 },
1284 { SK_Scalar1, SK_ScalarTanPIOver8 },
1285 { SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle },
1286 { SK_ScalarTanPIOver8, SK_Scalar1 },
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1287
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1288 { 0, SK_Scalar1 },
1289 { -SK_ScalarTanPIOver8, SK_Scalar1 },
1290 { -SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle },
1291 { -SK_Scalar1, SK_ScalarTanPIOver8 },
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1292
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1293 { -SK_Scalar1, 0 },
1294 { -SK_Scalar1, -SK_ScalarTanPIOver8 },
1295 { -SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle },
1296 { -SK_ScalarTanPIOver8, -SK_Scalar1 },
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1297
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1298 { 0, -SK_Scalar1 },
1299 { SK_ScalarTanPIOver8, -SK_Scalar1 },
1300 { SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle },
1301 { SK_Scalar1, -SK_ScalarTanPIOver8 },
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1302
robertphillips@google.com54820282012-10-18 15:26:12 +0000[diff] [blame]1303 { SK_Scalar1, 0 }
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1304};
1305
1306int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop,
1307 SkRotationDirection dir, const SkMatrix* userMatrix,
1308 SkPoint quadPoints[])
1309{
1310 // rotate by x,y so that uStart is (1.0)
1311 SkScalar x = SkPoint::DotProduct(uStart, uStop);
1312 SkScalar y = SkPoint::CrossProduct(uStart, uStop);
1313
1314 SkScalar absX = SkScalarAbs(x);
1315 SkScalar absY = SkScalarAbs(y);
1316
1317 int pointCount;
1318
1319 // check for (effectively) coincident vectors
1320 // this can happen if our angle is nearly 0 or nearly 180 (y == 0)
1321 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0)
1322 if (absY <= SK_ScalarNearlyZero && x > 0 &&
1323 ((y >= 0 && kCW_SkRotationDirection == dir) ||
1324 (y <= 0 && kCCW_SkRotationDirection == dir))) {
rmistry@google.com935e9f42012-08-23 18:09:54 +0000[diff] [blame]1325
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1326 // just return the start-point
1327 quadPoints[0].set(SK_Scalar1, 0);
1328 pointCount = 1;
1329 } else {
1330 if (dir == kCCW_SkRotationDirection)
1331 y = -y;
1332
1333 // what octant (quadratic curve) is [xy] in?
1334 int oct = 0;
1335 bool sameSign = true;
1336
1337 if (0 == y)
1338 {
1339 oct = 4; // 180
1340 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero);
1341 }
1342 else if (0 == x)
1343 {
1344 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero);
1345 if (y > 0)
1346 oct = 2; // 90
1347 else
1348 oct = 6; // 270
1349 }
1350 else
1351 {
1352 if (y < 0)
1353 oct += 4;
1354 if ((x < 0) != (y < 0))
1355 {
1356 oct += 2;
1357 sameSign = false;
1358 }
1359 if ((absX < absY) == sameSign)
1360 oct += 1;
1361 }
1362
1363 int wholeCount = oct << 1;
1364 memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint));
1365
1366 const SkPoint* arc = &gQuadCirclePts[wholeCount];
robertphillips@google.comdd4bd432013-07-09 15:03:59 +0000[diff] [blame^]1367 if (truncate_last_curve(arc, x, y, &quadPoints[wholeCount + 1]))
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1368 {
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1369 wholeCount += 2;
1370 }
1371 pointCount = wholeCount + 1;
1372 }
1373
1374 // now handle counter-clockwise and the initial unitStart rotation
1375 SkMatrix matrix;
1376 matrix.setSinCos(uStart.fY, uStart.fX);
1377 if (dir == kCCW_SkRotationDirection) {
1378 matrix.preScale(SK_Scalar1, -SK_Scalar1);
1379 }
1380 if (userMatrix) {
1381 matrix.postConcat(*userMatrix);
1382 }
1383 matrix.mapPoints(quadPoints, pointCount);
1384 return pointCount;
1385}
reed@google.com79975882013-04-12 19:11:10 +0000[diff] [blame]1386
1387///////////////////////////////////////////////////////////////////////////////
1388
reed@google.com256f3102013-04-16 21:07:27 +0000[diff] [blame]1389// F = (A (1 - t)^2 + C t^2 + 2 B (1 - t) t w)
1390// ------------------------------------------
1391// ((1 - t)^2 + t^2 + 2 (1 - t) t w)
1392//
1393// = {t^2 (P0 + P2 - 2 P1 w), t (-2 P0 + 2 P1 w), P0}
1394// ------------------------------------------------
1395// {t^2 (2 - 2 w), t (-2 + 2 w), 1}
1396//
reed@google.com256f3102013-04-16 21:07:27 +0000[diff] [blame]1397
1398// Take the parametric specification for the conic (either X or Y) and return
1399// in coeff[] the coefficients for the simple quadratic polynomial
1400// coeff[0] for t^2
1401// coeff[1] for t
1402// coeff[2] for constant term
1403//
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1404static SkScalar conic_eval_pos(const SkScalar src[], SkScalar w, SkScalar t) {
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1405 SkASSERT(src);
1406 SkASSERT(t >= 0 && t <= SK_Scalar1);
skia.committer@gmail.comcd487462013-04-17 07:00:56 +0000[diff] [blame]1407
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1408 SkScalar src2w = SkScalarMul(src[2], w);
1409 SkScalar C = src[0];
1410 SkScalar A = src[4] - 2 * src2w + C;
1411 SkScalar B = 2 * (src2w - C);
1412 SkScalar numer = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
skia.committer@gmail.comcd487462013-04-17 07:00:56 +0000[diff] [blame]1413
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1414 B = 2 * (w - SK_Scalar1);
1415 C = SK_Scalar1;
1416 A = -B;
1417 SkScalar denom = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
skia.committer@gmail.comcd487462013-04-17 07:00:56 +0000[diff] [blame]1418
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1419 return SkScalarDiv(numer, denom);
1420}
1421
1422// F' = 2 (C t (1 + t (-1 + w)) - A (-1 + t) (t (-1 + w) - w) + B (1 - 2 t) w)
1423//
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1424// t^2 : (2 P0 - 2 P2 - 2 P0 w + 2 P2 w)
1425// t^1 : (-2 P0 + 2 P2 + 4 P0 w - 4 P1 w)
1426// t^0 : -2 P0 w + 2 P1 w
1427//
1428// We disregard magnitude, so we can freely ignore the denominator of F', and
1429// divide the numerator by 2
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1430//
reed@google.com256f3102013-04-16 21:07:27 +0000[diff] [blame]1431// coeff[0] for t^2
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1432// coeff[1] for t^1
1433// coeff[2] for t^0
reed@google.com256f3102013-04-16 21:07:27 +0000[diff] [blame]1434//
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1435static void conic_deriv_coeff(const SkScalar src[], SkScalar w, SkScalar coeff[3]) {
1436 const SkScalar P20 = src[4] - src[0];
1437 const SkScalar P10 = src[2] - src[0];
1438 const SkScalar wP10 = w * P10;
1439 coeff[0] = w * P20 - P20;
1440 coeff[1] = P20 - 2 * wP10;
1441 coeff[2] = wP10;
reed@google.com256f3102013-04-16 21:07:27 +0000[diff] [blame]1442}
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1443
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1444static SkScalar conic_eval_tan(const SkScalar coord[], SkScalar w, SkScalar t) {
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1445 SkScalar coeff[3];
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1446 conic_deriv_coeff(coord, w, coeff);
1447 return t * (t * coeff[0] + coeff[1]) + coeff[2];
1448}
1449
1450static bool conic_find_extrema(const SkScalar src[], SkScalar w, SkScalar* t) {
1451 SkScalar coeff[3];
1452 conic_deriv_coeff(src, w, coeff);
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1453
1454 SkScalar tValues[2];
1455 int roots = SkFindUnitQuadRoots(coeff[0], coeff[1], coeff[2], tValues);
1456 SkASSERT(0 == roots || 1 == roots);
skia.committer@gmail.comcd487462013-04-17 07:00:56 +0000[diff] [blame]1457
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1458 if (1 == roots) {
1459 *t = tValues[0];
1460 return true;
1461 }
1462 return false;
1463}
reed@google.com256f3102013-04-16 21:07:27 +0000[diff] [blame]1464
reed@google.com0fef2ef2013-04-12 21:55:26 +0000[diff] [blame]1465struct SkP3D {
1466 SkScalar fX, fY, fZ;
skia.committer@gmail.coma43230c2013-04-13 07:01:15 +0000[diff] [blame]1467
reed@google.com0fef2ef2013-04-12 21:55:26 +0000[diff] [blame]1468 void set(SkScalar x, SkScalar y, SkScalar z) {
1469 fX = x; fY = y; fZ = z;
1470 }
skia.committer@gmail.coma43230c2013-04-13 07:01:15 +0000[diff] [blame]1471
reed@google.com0fef2ef2013-04-12 21:55:26 +0000[diff] [blame]1472 void projectDown(SkPoint* dst) const {
1473 dst->set(fX / fZ, fY / fZ);
1474 }
1475};
1476
1477// we just return the middle 3 points, since the first and last are dups of src
1478//
1479static void p3d_interp(const SkScalar src[3], SkScalar dst[3], SkScalar t) {
1480 SkScalar ab = SkScalarInterp(src[0], src[3], t);
1481 SkScalar bc = SkScalarInterp(src[3], src[6], t);
1482 dst[0] = ab;
1483 dst[3] = SkScalarInterp(ab, bc, t);
1484 dst[6] = bc;
1485}
1486
1487static void ratquad_mapTo3D(const SkPoint src[3], SkScalar w, SkP3D dst[]) {
1488 dst[0].set(src[0].fX * 1, src[0].fY * 1, 1);
1489 dst[1].set(src[1].fX * w, src[1].fY * w, w);
1490 dst[2].set(src[2].fX * 1, src[2].fY * 1, 1);
1491}
1492
reed@google.com4e1502a2013-05-07 20:42:35 +0000[diff] [blame]1493void SkConic::evalAt(SkScalar t, SkPoint* pt, SkVector* tangent) const {
reed@google.com79975882013-04-12 19:11:10 +0000[diff] [blame]1494 SkASSERT(t >= 0 && t <= SK_Scalar1);
skia.committer@gmail.coma43230c2013-04-13 07:01:15 +0000[diff] [blame]1495
reed@google.com79975882013-04-12 19:11:10 +0000[diff] [blame]1496 if (pt) {
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1497 pt->set(conic_eval_pos(&fPts[0].fX, fW, t),
1498 conic_eval_pos(&fPts[0].fY, fW, t));
reed@google.com79975882013-04-12 19:11:10 +0000[diff] [blame]1499 }
reed@google.com4e1502a2013-05-07 20:42:35 +0000[diff] [blame]1500 if (tangent) {
1501 tangent->set(conic_eval_tan(&fPts[0].fX, fW, t),
1502 conic_eval_tan(&fPts[0].fY, fW, t));
1503 }
reed@google.com79975882013-04-12 19:11:10 +0000[diff] [blame]1504}
1505
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1506void SkConic::chopAt(SkScalar t, SkConic dst[2]) const {
reed@google.com0fef2ef2013-04-12 21:55:26 +0000[diff] [blame]1507 SkP3D tmp[3], tmp2[3];
1508
1509 ratquad_mapTo3D(fPts, fW, tmp);
skia.committer@gmail.coma43230c2013-04-13 07:01:15 +0000[diff] [blame]1510
reed@google.com0fef2ef2013-04-12 21:55:26 +0000[diff] [blame]1511 p3d_interp(&tmp[0].fX, &tmp2[0].fX, t);
1512 p3d_interp(&tmp[0].fY, &tmp2[0].fY, t);
1513 p3d_interp(&tmp[0].fZ, &tmp2[0].fZ, t);
skia.committer@gmail.coma43230c2013-04-13 07:01:15 +0000[diff] [blame]1514
reed@google.com0fef2ef2013-04-12 21:55:26 +0000[diff] [blame]1515 dst[0].fPts[0] = fPts[0];
1516 tmp2[0].projectDown(&dst[0].fPts[1]);
1517 tmp2[1].projectDown(&dst[0].fPts[2]); dst[1].fPts[0] = dst[0].fPts[2];
1518 tmp2[2].projectDown(&dst[1].fPts[1]);
1519 dst[1].fPts[2] = fPts[2];
1520
mike@reedtribe.org464a1532013-04-13 10:51:51 +0000[diff] [blame]1521 // to put in "standard form", where w0 and w2 are both 1, we compute the
1522 // new w1 as sqrt(w1*w1/w0*w2)
1523 // or
1524 // w1 /= sqrt(w0*w2)
1525 //
1526 // However, in our case, we know that for dst[0], w0 == 1, and for dst[1], w2 == 1
1527 //
1528 SkScalar root = SkScalarSqrt(tmp2[1].fZ);
1529 dst[0].fW = tmp2[0].fZ / root;
1530 dst[1].fW = tmp2[2].fZ / root;
reed@google.com79975882013-04-12 19:11:10 +0000[diff] [blame]1531}
mike@reedtribe.org91068392013-04-14 02:40:50 +0000[diff] [blame]1532
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1533static SkScalar subdivide_w_value(SkScalar w) {
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1534 return SkScalarSqrt(SK_ScalarHalf + w * SK_ScalarHalf);
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1535}
1536
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1537void SkConic::chop(SkConic dst[2]) const {
mike@reedtribe.org91068392013-04-14 02:40:50 +0000[diff] [blame]1538 SkScalar scale = SkScalarInvert(SK_Scalar1 + fW);
1539 SkScalar p1x = fW * fPts[1].fX;
1540 SkScalar p1y = fW * fPts[1].fY;
1541 SkScalar mx = (fPts[0].fX + 2 * p1x + fPts[2].fX) * scale * SK_ScalarHalf;
1542 SkScalar my = (fPts[0].fY + 2 * p1y + fPts[2].fY) * scale * SK_ScalarHalf;
1543
1544 dst[0].fPts[0] = fPts[0];
1545 dst[0].fPts[1].set((fPts[0].fX + p1x) * scale,
1546 (fPts[0].fY + p1y) * scale);
1547 dst[0].fPts[2].set(mx, my);
1548
1549 dst[1].fPts[0].set(mx, my);
1550 dst[1].fPts[1].set((p1x + fPts[2].fX) * scale,
1551 (p1y + fPts[2].fY) * scale);
1552 dst[1].fPts[2] = fPts[2];
1553
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1554 dst[0].fW = dst[1].fW = subdivide_w_value(fW);
mike@reedtribe.org91068392013-04-14 02:40:50 +0000[diff] [blame]1555}
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1556
mike@reedtribe.orgdc5176e2013-04-27 18:23:16 +0000[diff] [blame]1557/*
1558 * "High order approximation of conic sections by quadratic splines"
1559 * by Michael Floater, 1993
1560 */
mike@reedtribe.org3661c442013-04-30 02:14:58 +0000[diff] [blame]1561#define AS_QUAD_ERROR_SETUP \
1562 SkScalar a = fW - 1; \
1563 SkScalar k = a / (4 * (2 + a)); \
1564 SkScalar x = k * (fPts[0].fX - 2 * fPts[1].fX + fPts[2].fX); \
1565 SkScalar y = k * (fPts[0].fY - 2 * fPts[1].fY + fPts[2].fY);
1566
1567void SkConic::computeAsQuadError(SkVector* err) const {
1568 AS_QUAD_ERROR_SETUP
1569 err->set(x, y);
1570}
1571
1572bool SkConic::asQuadTol(SkScalar tol) const {
1573 AS_QUAD_ERROR_SETUP
1574 return (x * x + y * y) <= tol * tol;
mike@reedtribe.orgdc5176e2013-04-27 18:23:16 +0000[diff] [blame]1575}
1576
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1577int SkConic::computeQuadPOW2(SkScalar tol) const {
mike@reedtribe.org3661c442013-04-30 02:14:58 +0000[diff] [blame]1578 AS_QUAD_ERROR_SETUP
1579 SkScalar error = SkScalarSqrt(x * x + y * y) - tol;
1580
1581 if (error <= 0) {
mike@reedtribe.orgdc5176e2013-04-27 18:23:16 +0000[diff] [blame]1582 return 0;
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1583 }
mike@reedtribe.orgdc5176e2013-04-27 18:23:16 +0000[diff] [blame]1584 uint32_t ierr = (uint32_t)error;
mike@reedtribe.org3661c442013-04-30 02:14:58 +0000[diff] [blame]1585 return (34 - SkCLZ(ierr)) >> 1;
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1586}
1587
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1588static SkPoint* subdivide(const SkConic& src, SkPoint pts[], int level) {
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1589 SkASSERT(level >= 0);
mike@reedtribe.org3661c442013-04-30 02:14:58 +0000[diff] [blame]1590
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1591 if (0 == level) {
1592 memcpy(pts, &src.fPts[1], 2 * sizeof(SkPoint));
1593 return pts + 2;
1594 } else {
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1595 SkConic dst[2];
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1596 src.chop(dst);
1597 --level;
1598 pts = subdivide(dst[0], pts, level);
1599 return subdivide(dst[1], pts, level);
1600 }
1601}
1602
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1603int SkConic::chopIntoQuadsPOW2(SkPoint pts[], int pow2) const {
mike@reedtribe.org3661c442013-04-30 02:14:58 +0000[diff] [blame]1604 SkASSERT(pow2 >= 0);
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1605 *pts = fPts[0];
reed@google.com901641b2013-04-15 15:23:38 +0000[diff] [blame]1606 SkDEBUGCODE(SkPoint* endPts =) subdivide(*this, pts + 1, pow2);
mike@reedtribe.org6c125bd2013-04-15 15:20:52 +0000[diff] [blame]1607 SkASSERT(endPts - pts == (2 * (1 << pow2) + 1));
1608 return 1 << pow2;
1609}
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1610
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1611bool SkConic::findXExtrema(SkScalar* t) const {
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1612 return conic_find_extrema(&fPts[0].fX, fW, t);
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1613}
1614
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1615bool SkConic::findYExtrema(SkScalar* t) const {
mike@reedtribe.org3f0cf542013-05-08 01:55:49 +0000[diff] [blame]1616 return conic_find_extrema(&fPts[0].fY, fW, t);
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1617}
1618
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1619bool SkConic::chopAtXExtrema(SkConic dst[2]) const {
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1620 SkScalar t;
1621 if (this->findXExtrema(&t)) {
1622 this->chopAt(t, dst);
1623 // now clean-up the middle, since we know t was meant to be at
1624 // an X-extrema
1625 SkScalar value = dst[0].fPts[2].fX;
1626 dst[0].fPts[1].fX = value;
1627 dst[1].fPts[0].fX = value;
1628 dst[1].fPts[1].fX = value;
1629 return true;
1630 }
1631 return false;
1632}
1633
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1634bool SkConic::chopAtYExtrema(SkConic dst[2]) const {
mike@reedtribe.org472a49c2013-04-17 01:21:01 +0000[diff] [blame]1635 SkScalar t;
1636 if (this->findYExtrema(&t)) {
1637 this->chopAt(t, dst);
1638 // now clean-up the middle, since we know t was meant to be at
1639 // an Y-extrema
1640 SkScalar value = dst[0].fPts[2].fY;
1641 dst[0].fPts[1].fY = value;
1642 dst[1].fPts[0].fY = value;
1643 dst[1].fPts[1].fY = value;
1644 return true;
1645 }
1646 return false;
1647}
1648
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1649void SkConic::computeTightBounds(SkRect* bounds) const {
mike@reedtribe.orgc7325912013-04-17 02:25:33 +0000[diff] [blame]1650 SkPoint pts[4];
1651 pts[0] = fPts[0];
1652 pts[1] = fPts[2];
1653 int count = 2;
1654
1655 SkScalar t;
1656 if (this->findXExtrema(&t)) {
1657 this->evalAt(t, &pts[count++]);
1658 }
1659 if (this->findYExtrema(&t)) {
1660 this->evalAt(t, &pts[count++]);
1661 }
1662 bounds->set(pts, count);
1663}
1664
mike@reedtribe.org27a0a562013-04-26 00:58:29 +0000[diff] [blame]1665void SkConic::computeFastBounds(SkRect* bounds) const {
mike@reedtribe.orgc7325912013-04-17 02:25:33 +0000[diff] [blame]1666 bounds->set(fPts, 3);
1667}