[go: nahoru, domu]
Google Git
Sign in
android / platform / external / chromium_org / third_party / skia / src / 001bd97fe25bbaefe62915ac6f803b765230282b / . / core / SkGeometry.cpp
blob: 483c08e35be4b2006fad44634d97b8b831f4f672 [file] [log] [blame]
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]1/* libs/graphics/sgl/SkGeometry.cpp
2**
3** Copyright 2006, The Android Open Source Project
4**
5** Licensed under the Apache License, Version 2.0 (the "License");
6** you may not use this file except in compliance with the License.
7** You may obtain a copy of the License at
8**
9** http://www.apache.org/licenses/LICENSE-2.0
10**
11** Unless required by applicable law or agreed to in writing, software
12** distributed under the License is distributed on an "AS IS" BASIS,
13** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14** See the License for the specific language governing permissions and
15** limitations under the License.
16*/
17
18#include "SkGeometry.h"
19#include "Sk64.h"
20#include "SkMatrix.h"
21
22/** If defined, this makes eval_quad and eval_cubic do more setup (sometimes
23 involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul.
24 May also introduce overflow of fixed when we compute our setup.
25*/
26#ifdef SK_SCALAR_IS_FIXED
27 #define DIRECT_EVAL_OF_POLYNOMIALS
28#endif
29
30////////////////////////////////////////////////////////////////////////
31
32#ifdef SK_SCALAR_IS_FIXED
33 static int is_not_monotonic(int a, int b, int c, int d)
34 {
35 return (((a - b) | (b - c) | (c - d)) & ((b - a) | (c - b) | (d - c))) >> 31;
36 }
37
38 static int is_not_monotonic(int a, int b, int c)
39 {
40 return (((a - b) | (b - c)) & ((b - a) | (c - b))) >> 31;
41 }
42#else
43 static int is_not_monotonic(float a, float b, float c)
44 {
45 float ab = a - b;
46 float bc = b - c;
47 if (ab < 0)
48 bc = -bc;
49 return ab == 0 || bc < 0;
50 }
51#endif
52
53////////////////////////////////////////////////////////////////////////
54
55static bool is_unit_interval(SkScalar x)
56{
57 return x > 0 && x < SK_Scalar1;
58}
59
60static int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio)
61{
62 SkASSERT(ratio);
63
64 if (numer < 0)
65 {
66 numer = -numer;
67 denom = -denom;
68 }
69
70 if (denom == 0 || numer == 0 || numer >= denom)
71 return 0;
72
73 SkScalar r = SkScalarDiv(numer, denom);
74 SkASSERT(r >= 0 && r < SK_Scalar1);
75 if (r == 0) // catch underflow if numer <<<< denom
76 return 0;
77 *ratio = r;
78 return 1;
79}
80
81/** From Numerical Recipes in C.
82
83 Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C])
84 x1 = Q / A
85 x2 = C / Q
86*/
87int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2])
88{
89 SkASSERT(roots);
90
91 if (A == 0)
92 return valid_unit_divide(-C, B, roots);
93
94 SkScalar* r = roots;
95
96#ifdef SK_SCALAR_IS_FLOAT
97 float R = B*B - 4*A*C;
98 if (R < 0) // complex roots
99 return 0;
100 R = sk_float_sqrt(R);
101#else
102 Sk64 RR, tmp;
103
104 RR.setMul(B,B);
105 tmp.setMul(A,C);
106 tmp.shiftLeft(2);
107 RR.sub(tmp);
108 if (RR.isNeg())
109 return 0;
110 SkFixed R = RR.getSqrt();
111#endif
112
113 SkScalar Q = (B < 0) ? -(B-R)/2 : -(B+R)/2;
114 r += valid_unit_divide(Q, A, r);
115 r += valid_unit_divide(C, Q, r);
116 if (r - roots == 2)
117 {
118 if (roots[0] > roots[1])
119 SkTSwap<SkScalar>(roots[0], roots[1]);
120 else if (roots[0] == roots[1]) // nearly-equal?
121 r -= 1; // skip the double root
122 }
123 return (int)(r - roots);
124}
125
126#ifdef SK_SCALAR_IS_FIXED
127/** Trim A/B/C down so that they are all <= 32bits
128 and then call SkFindUnitQuadRoots()
129*/
130static int Sk64FindFixedQuadRoots(const Sk64& A, const Sk64& B, const Sk64& C, SkFixed roots[2])
131{
132 int na = A.shiftToMake32();
133 int nb = B.shiftToMake32();
134 int nc = C.shiftToMake32();
135
136 int shift = SkMax32(na, SkMax32(nb, nc));
137 SkASSERT(shift >= 0);
138
139 return SkFindUnitQuadRoots(A.getShiftRight(shift), B.getShiftRight(shift), C.getShiftRight(shift), roots);
140}
141#endif
142
143/////////////////////////////////////////////////////////////////////////////////////
144/////////////////////////////////////////////////////////////////////////////////////
145
146static SkScalar eval_quad(const SkScalar src[], SkScalar t)
147{
148 SkASSERT(src);
149 SkASSERT(t >= 0 && t <= SK_Scalar1);
150
151#ifdef DIRECT_EVAL_OF_POLYNOMIALS
152 SkScalar C = src[0];
153 SkScalar A = src[4] - 2 * src[2] + C;
154 SkScalar B = 2 * (src[2] - C);
155 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
156#else
157 SkScalar ab = SkScalarInterp(src[0], src[2], t);
158 SkScalar bc = SkScalarInterp(src[2], src[4], t);
159 return SkScalarInterp(ab, bc, t);
160#endif
161}
162
163static SkScalar eval_quad_derivative(const SkScalar src[], SkScalar t)
164{
165 SkScalar A = src[4] - 2 * src[2] + src[0];
166 SkScalar B = src[2] - src[0];
167
168 return 2 * SkScalarMulAdd(A, t, B);
169}
170
171static SkScalar eval_quad_derivative_at_half(const SkScalar src[])
172{
173 SkScalar A = src[4] - 2 * src[2] + src[0];
174 SkScalar B = src[2] - src[0];
175 return A + 2 * B;
176}
177
178void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt, SkVector* tangent)
179{
180 SkASSERT(src);
181 SkASSERT(t >= 0 && t <= SK_Scalar1);
182
183 if (pt)
184 pt->set(eval_quad(&src[0].fX, t), eval_quad(&src[0].fY, t));
185 if (tangent)
186 tangent->set(eval_quad_derivative(&src[0].fX, t),
187 eval_quad_derivative(&src[0].fY, t));
188}
189
190void SkEvalQuadAtHalf(const SkPoint src[3], SkPoint* pt, SkVector* tangent)
191{
192 SkASSERT(src);
193
194 if (pt)
195 {
196 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
197 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
198 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
199 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
200 pt->set(SkScalarAve(x01, x12), SkScalarAve(y01, y12));
201 }
202 if (tangent)
203 tangent->set(eval_quad_derivative_at_half(&src[0].fX),
204 eval_quad_derivative_at_half(&src[0].fY));
205}
206
207static void interp_quad_coords(const SkScalar* src, SkScalar* dst, SkScalar t)
208{
209 SkScalar ab = SkScalarInterp(src[0], src[2], t);
210 SkScalar bc = SkScalarInterp(src[2], src[4], t);
211
212 dst[0] = src[0];
213 dst[2] = ab;
214 dst[4] = SkScalarInterp(ab, bc, t);
215 dst[6] = bc;
216 dst[8] = src[4];
217}
218
219void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t)
220{
221 SkASSERT(t > 0 && t < SK_Scalar1);
222
223 interp_quad_coords(&src[0].fX, &dst[0].fX, t);
224 interp_quad_coords(&src[0].fY, &dst[0].fY, t);
225}
226
227void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5])
228{
229 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
230 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
231 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
232 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
233
234 dst[0] = src[0];
235 dst[1].set(x01, y01);
236 dst[2].set(SkScalarAve(x01, x12), SkScalarAve(y01, y12));
237 dst[3].set(x12, y12);
238 dst[4] = src[2];
239}
240
241/** Quad'(t) = At + B, where
242 A = 2(a - 2b + c)
243 B = 2(b - a)
244 Solve for t, only if it fits between 0 < t < 1
245*/
246int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1])
247{
248 /* At + B == 0
249 t = -B / A
250 */
251#ifdef SK_SCALAR_IS_FIXED
252 return is_not_monotonic(a, b, c) && valid_unit_divide(a - b, a - b - b + c, tValue);
253#else
254 return valid_unit_divide(a - b, a - b - b + c, tValue);
255#endif
256}
257
reed@android.come5dd6cd2009-01-15 14:38:33 +0000[diff] [blame]258static inline void flatten_double_quad_extrema(SkScalar coords[14])
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]259{
260 coords[2] = coords[6] = coords[4];
261}
262
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]263/* 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^]264 stored in dst[]. Guarantees that the 1/2 quads will be monotonic.
265 */
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]266int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5])
267{
268 SkASSERT(src);
269 SkASSERT(dst);
reed@android.com001bd972009-11-17 18:47:52 +0000[diff] [blame^]270
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]271#if 0
272 static bool once = true;
273 if (once)
274 {
275 once = false;
276 SkPoint s[3] = { 0, 26398, 0, 26331, 0, 20621428 };
277 SkPoint d[6];
278
279 int n = SkChopQuadAtYExtrema(s, d);
280 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);
281 }
282#endif
reed@android.com001bd972009-11-17 18:47:52 +0000[diff] [blame^]283
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]284 SkScalar a = src[0].fY;
285 SkScalar b = src[1].fY;
286 SkScalar c = src[2].fY;
reed@android.com001bd972009-11-17 18:47:52 +0000[diff] [blame^]287
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]288 if (is_not_monotonic(a, b, c))
289 {
290 SkScalar tValue;
291 if (valid_unit_divide(a - b, a - b - b + c, &tValue))
292 {
293 SkChopQuadAt(src, dst, tValue);
294 flatten_double_quad_extrema(&dst[0].fY);
295 return 1;
296 }
297 // if we get here, we need to force dst to be monotonic, even though
298 // we couldn't compute a unit_divide value (probably underflow).
299 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c;
300 }
301 dst[0].set(src[0].fX, a);
302 dst[1].set(src[1].fX, b);
303 dst[2].set(src[2].fX, c);
304 return 0;
305}
306
reed@android.com001bd972009-11-17 18:47:52 +0000[diff] [blame^]307/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is
308 stored in dst[]. Guarantees that the 1/2 quads will be monotonic.
309 */
310int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5])
311{
312 SkASSERT(src);
313 SkASSERT(dst);
314
315 SkScalar a = src[0].fX;
316 SkScalar b = src[1].fX;
317 SkScalar c = src[2].fX;
318
319 if (is_not_monotonic(a, b, c)) {
320 SkScalar tValue;
321 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) {
322 SkChopQuadAt(src, dst, tValue);
323 flatten_double_quad_extrema(&dst[0].fX);
324 return 1;
325 }
326 // if we get here, we need to force dst to be monotonic, even though
327 // we couldn't compute a unit_divide value (probably underflow).
328 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c;
329 }
330 dst[0].set(a, src[0].fY);
331 dst[1].set(b, src[1].fY);
332 dst[2].set(c, src[2].fY);
333 return 0;
334}
335
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]336// F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2
337// F'(t) = 2 (b - a) + 2 (a - 2b + c) t
338// F''(t) = 2 (a - 2b + c)
339//
340// A = 2 (b - a)
341// B = 2 (a - 2b + c)
342//
343// Maximum curvature for a quadratic means solving
344// Fx' Fx'' + Fy' Fy'' = 0
345//
346// t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2)
347//
348int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5])
349{
350 SkScalar Ax = src[1].fX - src[0].fX;
351 SkScalar Ay = src[1].fY - src[0].fY;
352 SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX;
353 SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY;
354 SkScalar t = 0; // 0 means don't chop
355
356#ifdef SK_SCALAR_IS_FLOAT
357 (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t);
358#else
359 // !!! should I use SkFloat here? seems like it
360 Sk64 numer, denom, tmp;
361
362 numer.setMul(Ax, -Bx);
363 tmp.setMul(Ay, -By);
364 numer.add(tmp);
365
366 if (numer.isPos()) // do nothing if numer <= 0
367 {
368 denom.setMul(Bx, Bx);
369 tmp.setMul(By, By);
370 denom.add(tmp);
371 SkASSERT(!denom.isNeg());
372 if (numer < denom)
373 {
374 t = numer.getFixedDiv(denom);
375 SkASSERT(t >= 0 && t <= SK_Fixed1); // assert that we're numerically stable (ha!)
376 if ((unsigned)t >= SK_Fixed1) // runtime check for numerical stability
377 t = 0; // ignore the chop
378 }
379 }
380#endif
381
382 if (t == 0)
383 {
384 memcpy(dst, src, 3 * sizeof(SkPoint));
385 return 1;
386 }
387 else
388 {
389 SkChopQuadAt(src, dst, t);
390 return 2;
391 }
392}
393
394////////////////////////////////////////////////////////////////////////////////////////
395///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS /////
396////////////////////////////////////////////////////////////////////////////////////////
397
398static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4])
399{
400 coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0];
401 coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]);
402 coeff[2] = 3*(pt[2] - pt[0]);
403 coeff[3] = pt[0];
404}
405
406void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4])
407{
408 SkASSERT(pts);
409
410 if (cx)
411 get_cubic_coeff(&pts[0].fX, cx);
412 if (cy)
413 get_cubic_coeff(&pts[0].fY, cy);
414}
415
416static SkScalar eval_cubic(const SkScalar src[], SkScalar t)
417{
418 SkASSERT(src);
419 SkASSERT(t >= 0 && t <= SK_Scalar1);
420
421 if (t == 0)
422 return src[0];
423
424#ifdef DIRECT_EVAL_OF_POLYNOMIALS
425 SkScalar D = src[0];
426 SkScalar A = src[6] + 3*(src[2] - src[4]) - D;
427 SkScalar B = 3*(src[4] - src[2] - src[2] + D);
428 SkScalar C = 3*(src[2] - D);
429
430 return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D);
431#else
432 SkScalar ab = SkScalarInterp(src[0], src[2], t);
433 SkScalar bc = SkScalarInterp(src[2], src[4], t);
434 SkScalar cd = SkScalarInterp(src[4], src[6], t);
435 SkScalar abc = SkScalarInterp(ab, bc, t);
436 SkScalar bcd = SkScalarInterp(bc, cd, t);
437 return SkScalarInterp(abc, bcd, t);
438#endif
439}
440
441/** return At^2 + Bt + C
442*/
443static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t)
444{
445 SkASSERT(t >= 0 && t <= SK_Scalar1);
446
447 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
448}
449
450static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t)
451{
452 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0];
453 SkScalar B = 2*(src[4] - 2 * src[2] + src[0]);
454 SkScalar C = src[2] - src[0];
455
456 return eval_quadratic(A, B, C, t);
457}
458
459static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t)
460{
461 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0];
462 SkScalar B = src[4] - 2 * src[2] + src[0];
463
464 return SkScalarMulAdd(A, t, B);
465}
466
467void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature)
468{
469 SkASSERT(src);
470 SkASSERT(t >= 0 && t <= SK_Scalar1);
471
472 if (loc)
473 loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t));
474 if (tangent)
475 tangent->set(eval_cubic_derivative(&src[0].fX, t),
476 eval_cubic_derivative(&src[0].fY, t));
477 if (curvature)
478 curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t),
479 eval_cubic_2ndDerivative(&src[0].fY, t));
480}
481
482/** Cubic'(t) = At^2 + Bt + C, where
483 A = 3(-a + 3(b - c) + d)
484 B = 6(a - 2b + c)
485 C = 3(b - a)
486 Solve for t, keeping only those that fit betwee 0 < t < 1
487*/
488int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2])
489{
490#ifdef SK_SCALAR_IS_FIXED
491 if (!is_not_monotonic(a, b, c, d))
492 return 0;
493#endif
494
495 // we divide A,B,C by 3 to simplify
496 SkScalar A = d - a + 3*(b - c);
497 SkScalar B = 2*(a - b - b + c);
498 SkScalar C = b - a;
499
500 return SkFindUnitQuadRoots(A, B, C, tValues);
501}
502
503static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t)
504{
505 SkScalar ab = SkScalarInterp(src[0], src[2], t);
506 SkScalar bc = SkScalarInterp(src[2], src[4], t);
507 SkScalar cd = SkScalarInterp(src[4], src[6], t);
508 SkScalar abc = SkScalarInterp(ab, bc, t);
509 SkScalar bcd = SkScalarInterp(bc, cd, t);
510 SkScalar abcd = SkScalarInterp(abc, bcd, t);
511
512 dst[0] = src[0];
513 dst[2] = ab;
514 dst[4] = abc;
515 dst[6] = abcd;
516 dst[8] = bcd;
517 dst[10] = cd;
518 dst[12] = src[6];
519}
520
521void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t)
522{
523 SkASSERT(t > 0 && t < SK_Scalar1);
524
525 interp_cubic_coords(&src[0].fX, &dst[0].fX, t);
526 interp_cubic_coords(&src[0].fY, &dst[0].fY, t);
527}
528
reed@android.com17bdc092009-08-28 20:06:54 +0000[diff] [blame]529/* http://code.google.com/p/skia/issues/detail?id=32
530
531 This test code would fail when we didn't check the return result of
532 valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is
533 that after the first chop, the parameters to valid_unit_divide are equal
534 (thanks to finite float precision and rounding in the subtracts). Thus
535 even though the 2nd tValue looks < 1.0, after we renormalize it, we end
536 up with 1.0, hence the need to check and just return the last cubic as
537 a degenerate clump of 4 points in the sampe place.
538
539 static void test_cubic() {
540 SkPoint src[4] = {
541 { 556.25000, 523.03003 },
542 { 556.23999, 522.96002 },
543 { 556.21997, 522.89001 },
544 { 556.21997, 522.82001 }
545 };
546 SkPoint dst[10];
547 SkScalar tval[] = { 0.33333334f, 0.99999994f };
548 SkChopCubicAt(src, dst, tval, 2);
549 }
550 */
551
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]552void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots)
553{
554#ifdef SK_DEBUG
555 {
556 for (int i = 0; i < roots - 1; i++)
557 {
558 SkASSERT(is_unit_interval(tValues[i]));
559 SkASSERT(is_unit_interval(tValues[i+1]));
560 SkASSERT(tValues[i] < tValues[i+1]);
561 }
562 }
563#endif
564
565 if (dst)
566 {
567 if (roots == 0) // nothing to chop
568 memcpy(dst, src, 4*sizeof(SkPoint));
569 else
570 {
571 SkScalar t = tValues[0];
572 SkPoint tmp[4];
573
574 for (int i = 0; i < roots; i++)
575 {
576 SkChopCubicAt(src, dst, t);
577 if (i == roots - 1)
578 break;
579
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]580 dst += 3;
reed@android.com17bdc092009-08-28 20:06:54 +0000[diff] [blame]581 // have src point to the remaining cubic (after the chop)
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]582 memcpy(tmp, dst, 4 * sizeof(SkPoint));
583 src = tmp;
reed@android.com17bdc092009-08-28 20:06:54 +0000[diff] [blame]584
585 // watch out in case the renormalized t isn't in range
586 if (!valid_unit_divide(tValues[i+1] - tValues[i],
587 SK_Scalar1 - tValues[i], &t)) {
588 // if we can't, just create a degenerate cubic
589 dst[4] = dst[5] = dst[6] = src[3];
590 break;
591 }
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]592 }
593 }
594 }
595}
596
597void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7])
598{
599 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
600 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
601 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
602 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
603 SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX);
604 SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY);
605
606 SkScalar x012 = SkScalarAve(x01, x12);
607 SkScalar y012 = SkScalarAve(y01, y12);
608 SkScalar x123 = SkScalarAve(x12, x23);
609 SkScalar y123 = SkScalarAve(y12, y23);
610
611 dst[0] = src[0];
612 dst[1].set(x01, y01);
613 dst[2].set(x012, y012);
614 dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123));
615 dst[4].set(x123, y123);
616 dst[5].set(x23, y23);
617 dst[6] = src[3];
618}
619
620static void flatten_double_cubic_extrema(SkScalar coords[14])
621{
622 coords[4] = coords[8] = coords[6];
623}
624
625/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that
626 the resulting beziers are monotonic in Y. This is called by the scan converter.
627 Depending on what is returned, dst[] is treated as follows
628 0 dst[0..3] is the original cubic
629 1 dst[0..3] and dst[3..6] are the two new cubics
630 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics
631 If dst == null, it is ignored and only the count is returned.
632*/
633int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10])
634{
635 SkScalar tValues[2];
636 int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY, src[3].fY, tValues);
637
638 SkChopCubicAt(src, dst, tValues, roots);
639 if (dst && roots > 0)
640 {
641 // we do some cleanup to ensure our Y extrema are flat
642 flatten_double_cubic_extrema(&dst[0].fY);
643 if (roots == 2)
644 flatten_double_cubic_extrema(&dst[3].fY);
645 }
646 return roots;
647}
648
649/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html
650
651 Inflection means that curvature is zero.
652 Curvature is [F' x F''] / [F'^3]
653 So we solve F'x X F''y - F'y X F''y == 0
654 After some canceling of the cubic term, we get
655 A = b - a
656 B = c - 2b + a
657 C = d - 3c + 3b - a
658 (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0
659*/
660int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[])
661{
662 SkScalar Ax = src[1].fX - src[0].fX;
663 SkScalar Ay = src[1].fY - src[0].fY;
664 SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX;
665 SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY;
666 SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX;
667 SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY;
668 int count;
669
670#ifdef SK_SCALAR_IS_FLOAT
671 count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues);
672#else
673 Sk64 A, B, C, tmp;
674
675 A.setMul(Bx, Cy);
676 tmp.setMul(By, Cx);
677 A.sub(tmp);
678
679 B.setMul(Ax, Cy);
680 tmp.setMul(Ay, Cx);
681 B.sub(tmp);
682
683 C.setMul(Ax, By);
684 tmp.setMul(Ay, Bx);
685 C.sub(tmp);
686
687 count = Sk64FindFixedQuadRoots(A, B, C, tValues);
688#endif
689
690 return count;
691}
692
693int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10])
694{
695 SkScalar tValues[2];
696 int count = SkFindCubicInflections(src, tValues);
697
698 if (dst)
699 {
700 if (count == 0)
701 memcpy(dst, src, 4 * sizeof(SkPoint));
702 else
703 SkChopCubicAt(src, dst, tValues, count);
704 }
705 return count + 1;
706}
707
708template <typename T> void bubble_sort(T array[], int count)
709{
710 for (int i = count - 1; i > 0; --i)
711 for (int j = i; j > 0; --j)
712 if (array[j] < array[j-1])
713 {
714 T tmp(array[j]);
715 array[j] = array[j-1];
716 array[j-1] = tmp;
717 }
718}
719
720#include "SkFP.h"
721
722// newton refinement
723#if 0
724static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root)
725{
726 // x1 = x0 - f(t) / f'(t)
727
728 SkFP T = SkScalarToFloat(root);
729 SkFP N, D;
730
731 // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2]
732 D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3);
733 D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2));
734 D = SkFPAdd(D, coeff[2]);
735
736 if (D == 0)
737 return root;
738
739 // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3]
740 N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]);
741 N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1]));
742 N = SkFPAdd(N, SkFPMul(T, coeff[2]));
743 N = SkFPAdd(N, coeff[3]);
744
745 if (N)
746 {
747 SkScalar delta = SkFPToScalar(SkFPDiv(N, D));
748
749 if (delta)
750 root -= delta;
751 }
752 return root;
753}
754#endif
755
756#if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop
757#pragma warning ( disable : 4702 )
758#endif
759
760/* Solve coeff(t) == 0, returning the number of roots that
761 lie withing 0 < t < 1.
762 coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3]
763*/
764static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3])
765{
766#ifndef SK_SCALAR_IS_FLOAT
767 return 0; // this is not yet implemented for software float
768#endif
769
770 if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic
771 {
772 return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues);
773 }
774
775 SkFP a, b, c, Q, R;
776
777 {
778 SkASSERT(coeff[0] != 0);
779
780 SkFP inva = SkFPInvert(coeff[0]);
781 a = SkFPMul(coeff[1], inva);
782 b = SkFPMul(coeff[2], inva);
783 c = SkFPMul(coeff[3], inva);
784 }
785 Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9);
786// R = (2*a*a*a - 9*a*b + 27*c) / 54;
787 R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2);
788 R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9));
789 R = SkFPAdd(R, SkFPMulInt(c, 27));
790 R = SkFPDivInt(R, 54);
791
792 SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q);
793 SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3);
794 SkFP adiv3 = SkFPDivInt(a, 3);
795
796 SkScalar* roots = tValues;
797 SkScalar r;
798
799 if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots
800 {
801#ifdef SK_SCALAR_IS_FLOAT
802 float theta = sk_float_acos(R / sk_float_sqrt(Q3));
803 float neg2RootQ = -2 * sk_float_sqrt(Q);
804
805 r = neg2RootQ * sk_float_cos(theta/3) - adiv3;
806 if (is_unit_interval(r))
807 *roots++ = r;
808
809 r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3;
810 if (is_unit_interval(r))
811 *roots++ = r;
812
813 r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3;
814 if (is_unit_interval(r))
815 *roots++ = r;
816
817 // now sort the roots
818 bubble_sort(tValues, (int)(roots - tValues));
819#endif
820 }
821 else // we have 1 real root
822 {
823 SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3));
824 A = SkFPCubeRoot(A);
825 if (SkFPGT(R, 0))
826 A = SkFPNeg(A);
827
828 if (A != 0)
829 A = SkFPAdd(A, SkFPDiv(Q, A));
830 r = SkFPToScalar(SkFPSub(A, adiv3));
831 if (is_unit_interval(r))
832 *roots++ = r;
833 }
834
835 return (int)(roots - tValues);
836}
837
838/* Looking for F' dot F'' == 0
839
840 A = b - a
841 B = c - 2b + a
842 C = d - 3c + 3b - a
843
844 F' = 3Ct^2 + 6Bt + 3A
845 F'' = 6Ct + 6B
846
847 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
848*/
849static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4])
850{
851 SkScalar a = src[2] - src[0];
852 SkScalar b = src[4] - 2 * src[2] + src[0];
853 SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0];
854
855 SkFP A = SkScalarToFP(a);
856 SkFP B = SkScalarToFP(b);
857 SkFP C = SkScalarToFP(c);
858
859 coeff[0] = SkFPMul(C, C);
860 coeff[1] = SkFPMulInt(SkFPMul(B, C), 3);
861 coeff[2] = SkFPMulInt(SkFPMul(B, B), 2);
862 coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A));
863 coeff[3] = SkFPMul(A, B);
864}
865
866// EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1
867//#define kMinTValueForChopping (SK_Scalar1 / 256)
868#define kMinTValueForChopping 0
869
870/* Looking for F' dot F'' == 0
871
872 A = b - a
873 B = c - 2b + a
874 C = d - 3c + 3b - a
875
876 F' = 3Ct^2 + 6Bt + 3A
877 F'' = 6Ct + 6B
878
879 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
880*/
881int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3])
882{
883 SkFP coeffX[4], coeffY[4];
884 int i;
885
886 formulate_F1DotF2(&src[0].fX, coeffX);
887 formulate_F1DotF2(&src[0].fY, coeffY);
888
889 for (i = 0; i < 4; i++)
890 coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]);
891
892 SkScalar t[3];
893 int count = solve_cubic_polynomial(coeffX, t);
894 int maxCount = 0;
895
896 // now remove extrema where the curvature is zero (mins)
897 // !!!! need a test for this !!!!
898 for (i = 0; i < count; i++)
899 {
900 // if (not_min_curvature())
901 if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping)
902 tValues[maxCount++] = t[i];
903 }
904 return maxCount;
905}
906
907int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3])
908{
909 SkScalar t_storage[3];
910
911 if (tValues == NULL)
912 tValues = t_storage;
913
914 int count = SkFindCubicMaxCurvature(src, tValues);
915
916 if (dst)
917 {
918 if (count == 0)
919 memcpy(dst, src, 4 * sizeof(SkPoint));
920 else
921 SkChopCubicAt(src, dst, tValues, count);
922 }
923 return count + 1;
924}
925
926////////////////////////////////////////////////////////////////////////////////
927
928/* Find t value for quadratic [a, b, c] = d.
929 Return 0 if there is no solution within [0, 1)
930*/
931static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d)
932{
933 // At^2 + Bt + C = d
934 SkScalar A = a - 2 * b + c;
935 SkScalar B = 2 * (b - a);
936 SkScalar C = a - d;
937
938 SkScalar roots[2];
939 int count = SkFindUnitQuadRoots(A, B, C, roots);
940
941 SkASSERT(count <= 1);
942 return count == 1 ? roots[0] : 0;
943}
944
945/* given a quad-curve and a point (x,y), chop the quad at that point and return
946 the new quad's offCurve point. Should only return false if the computed pos
947 is the start of the curve (i.e. root == 0)
948*/
949static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve)
950{
951 const SkScalar* base;
952 SkScalar value;
953
954 if (SkScalarAbs(x) < SkScalarAbs(y)) {
955 base = &quad[0].fX;
956 value = x;
957 } else {
958 base = &quad[0].fY;
959 value = y;
960 }
961
962 // note: this returns 0 if it thinks value is out of range, meaning the
963 // root might return something outside of [0, 1)
964 SkScalar t = quad_solve(base[0], base[2], base[4], value);
965
966 if (t > 0)
967 {
968 SkPoint tmp[5];
969 SkChopQuadAt(quad, tmp, t);
970 *offCurve = tmp[1];
971 return true;
972 } else {
973 /* t == 0 means either the value triggered a root outside of [0, 1)
974 For our purposes, we can ignore the <= 0 roots, but we want to
975 catch the >= 1 roots (which given our caller, will basically mean
976 a root of 1, give-or-take numerical instability). If we are in the
977 >= 1 case, return the existing offCurve point.
978
979 The test below checks to see if we are close to the "end" of the
980 curve (near base[4]). Rather than specifying a tolerance, I just
981 check to see if value is on to the right/left of the middle point
982 (depending on the direction/sign of the end points).
983 */
984 if ((base[0] < base[4] && value > base[2]) ||
985 (base[0] > base[4] && value < base[2])) // should root have been 1
986 {
987 *offCurve = quad[1];
988 return true;
989 }
990 }
991 return false;
992}
993
994static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = {
995 { SK_Scalar1, 0 },
996 { SK_Scalar1, SK_ScalarTanPIOver8 },
997 { SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 },
998 { SK_ScalarTanPIOver8, SK_Scalar1 },
999
1000 { 0, SK_Scalar1 },
1001 { -SK_ScalarTanPIOver8, SK_Scalar1 },
1002 { -SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 },
1003 { -SK_Scalar1, SK_ScalarTanPIOver8 },
1004
1005 { -SK_Scalar1, 0 },
1006 { -SK_Scalar1, -SK_ScalarTanPIOver8 },
1007 { -SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 },
1008 { -SK_ScalarTanPIOver8, -SK_Scalar1 },
1009
1010 { 0, -SK_Scalar1 },
1011 { SK_ScalarTanPIOver8, -SK_Scalar1 },
1012 { SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 },
1013 { SK_Scalar1, -SK_ScalarTanPIOver8 },
1014
1015 { SK_Scalar1, 0 }
1016};
1017
1018int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop,
1019 SkRotationDirection dir, const SkMatrix* userMatrix,
1020 SkPoint quadPoints[])
1021{
1022 // rotate by x,y so that uStart is (1.0)
1023 SkScalar x = SkPoint::DotProduct(uStart, uStop);
1024 SkScalar y = SkPoint::CrossProduct(uStart, uStop);
1025
1026 SkScalar absX = SkScalarAbs(x);
1027 SkScalar absY = SkScalarAbs(y);
1028
1029 int pointCount;
1030
1031 // check for (effectively) coincident vectors
1032 // this can happen if our angle is nearly 0 or nearly 180 (y == 0)
1033 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0)
1034 if (absY <= SK_ScalarNearlyZero && x > 0 &&
1035 ((y >= 0 && kCW_SkRotationDirection == dir) ||
1036 (y <= 0 && kCCW_SkRotationDirection == dir))) {
1037
1038 // just return the start-point
1039 quadPoints[0].set(SK_Scalar1, 0);
1040 pointCount = 1;
1041 } else {
1042 if (dir == kCCW_SkRotationDirection)
1043 y = -y;
1044
1045 // what octant (quadratic curve) is [xy] in?
1046 int oct = 0;
1047 bool sameSign = true;
1048
1049 if (0 == y)
1050 {
1051 oct = 4; // 180
1052 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero);
1053 }
1054 else if (0 == x)
1055 {
1056 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero);
1057 if (y > 0)
1058 oct = 2; // 90
1059 else
1060 oct = 6; // 270
1061 }
1062 else
1063 {
1064 if (y < 0)
1065 oct += 4;
1066 if ((x < 0) != (y < 0))
1067 {
1068 oct += 2;
1069 sameSign = false;
1070 }
1071 if ((absX < absY) == sameSign)
1072 oct += 1;
1073 }
1074
1075 int wholeCount = oct << 1;
1076 memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint));
1077
1078 const SkPoint* arc = &gQuadCirclePts[wholeCount];
1079 if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1]))
1080 {
1081 quadPoints[wholeCount + 2].set(x, y);
1082 wholeCount += 2;
1083 }
1084 pointCount = wholeCount + 1;
1085 }
1086
1087 // now handle counter-clockwise and the initial unitStart rotation
1088 SkMatrix matrix;
1089 matrix.setSinCos(uStart.fY, uStart.fX);
1090 if (dir == kCCW_SkRotationDirection) {
1091 matrix.preScale(SK_Scalar1, -SK_Scalar1);
1092 }
1093 if (userMatrix) {
1094 matrix.postConcat(*userMatrix);
1095 }
1096 matrix.mapPoints(quadPoints, pointCount);
1097 return pointCount;
1098}
1099