[go: nahoru, domu]
Google Git
Sign in
android / platform / external / chromium_org / third_party / skia / src / 68779c3e307771162ed2a55d2bfc9c63c807fc09 / . / core / SkGeometry.cpp
blob: 4704236fac76420fc4b951cc55fcfd9afc290f95 [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*/
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame^]633int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) {
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]634 SkScalar tValues[2];
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame^]635 int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY,
636 src[3].fY, tValues);
637
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]638 SkChopCubicAt(src, dst, tValues, roots);
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame^]639 if (dst && roots > 0) {
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]640 // we do some cleanup to ensure our Y extrema are flat
641 flatten_double_cubic_extrema(&dst[0].fY);
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame^]642 if (roots == 2) {
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]643 flatten_double_cubic_extrema(&dst[3].fY);
reed@android.com68779c32009-11-18 13:47:40 +0000[diff] [blame^]644 }
645 }
646 return roots;
647}
648
649int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) {
650 SkScalar tValues[2];
651 int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX,
652 src[3].fX, tValues);
653
654 SkChopCubicAt(src, dst, tValues, roots);
655 if (dst && roots > 0) {
656 // we do some cleanup to ensure our Y extrema are flat
657 flatten_double_cubic_extrema(&dst[0].fX);
658 if (roots == 2) {
659 flatten_double_cubic_extrema(&dst[3].fX);
660 }
reed@android.combcd4d5a2008-12-17 15:59:43 +0000[diff] [blame]661 }
662 return roots;
663}
664
665/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html
666
667 Inflection means that curvature is zero.
668 Curvature is [F' x F''] / [F'^3]
669 So we solve F'x X F''y - F'y X F''y == 0
670 After some canceling of the cubic term, we get
671 A = b - a
672 B = c - 2b + a
673 C = d - 3c + 3b - a
674 (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0
675*/
676int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[])
677{
678 SkScalar Ax = src[1].fX - src[0].fX;
679 SkScalar Ay = src[1].fY - src[0].fY;
680 SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX;
681 SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY;
682 SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX;
683 SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY;
684 int count;
685
686#ifdef SK_SCALAR_IS_FLOAT
687 count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues);
688#else
689 Sk64 A, B, C, tmp;
690
691 A.setMul(Bx, Cy);
692 tmp.setMul(By, Cx);
693 A.sub(tmp);
694
695 B.setMul(Ax, Cy);
696 tmp.setMul(Ay, Cx);
697 B.sub(tmp);
698
699 C.setMul(Ax, By);
700 tmp.setMul(Ay, Bx);
701 C.sub(tmp);
702
703 count = Sk64FindFixedQuadRoots(A, B, C, tValues);
704#endif
705
706 return count;
707}
708
709int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10])
710{
711 SkScalar tValues[2];
712 int count = SkFindCubicInflections(src, tValues);
713
714 if (dst)
715 {
716 if (count == 0)
717 memcpy(dst, src, 4 * sizeof(SkPoint));
718 else
719 SkChopCubicAt(src, dst, tValues, count);
720 }
721 return count + 1;
722}
723
724template <typename T> void bubble_sort(T array[], int count)
725{
726 for (int i = count - 1; i > 0; --i)
727 for (int j = i; j > 0; --j)
728 if (array[j] < array[j-1])
729 {
730 T tmp(array[j]);
731 array[j] = array[j-1];
732 array[j-1] = tmp;
733 }
734}
735
736#include "SkFP.h"
737
738// newton refinement
739#if 0
740static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root)
741{
742 // x1 = x0 - f(t) / f'(t)
743
744 SkFP T = SkScalarToFloat(root);
745 SkFP N, D;
746
747 // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2]
748 D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3);
749 D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2));
750 D = SkFPAdd(D, coeff[2]);
751
752 if (D == 0)
753 return root;
754
755 // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3]
756 N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]);
757 N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1]));
758 N = SkFPAdd(N, SkFPMul(T, coeff[2]));
759 N = SkFPAdd(N, coeff[3]);
760
761 if (N)
762 {
763 SkScalar delta = SkFPToScalar(SkFPDiv(N, D));
764
765 if (delta)
766 root -= delta;
767 }
768 return root;
769}
770#endif
771
772#if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop
773#pragma warning ( disable : 4702 )
774#endif
775
776/* Solve coeff(t) == 0, returning the number of roots that
777 lie withing 0 < t < 1.
778 coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3]
779*/
780static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3])
781{
782#ifndef SK_SCALAR_IS_FLOAT
783 return 0; // this is not yet implemented for software float
784#endif
785
786 if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic
787 {
788 return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues);
789 }
790
791 SkFP a, b, c, Q, R;
792
793 {
794 SkASSERT(coeff[0] != 0);
795
796 SkFP inva = SkFPInvert(coeff[0]);
797 a = SkFPMul(coeff[1], inva);
798 b = SkFPMul(coeff[2], inva);
799 c = SkFPMul(coeff[3], inva);
800 }
801 Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9);
802// R = (2*a*a*a - 9*a*b + 27*c) / 54;
803 R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2);
804 R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9));
805 R = SkFPAdd(R, SkFPMulInt(c, 27));
806 R = SkFPDivInt(R, 54);
807
808 SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q);
809 SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3);
810 SkFP adiv3 = SkFPDivInt(a, 3);
811
812 SkScalar* roots = tValues;
813 SkScalar r;
814
815 if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots
816 {
817#ifdef SK_SCALAR_IS_FLOAT
818 float theta = sk_float_acos(R / sk_float_sqrt(Q3));
819 float neg2RootQ = -2 * sk_float_sqrt(Q);
820
821 r = neg2RootQ * sk_float_cos(theta/3) - adiv3;
822 if (is_unit_interval(r))
823 *roots++ = r;
824
825 r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3;
826 if (is_unit_interval(r))
827 *roots++ = r;
828
829 r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3;
830 if (is_unit_interval(r))
831 *roots++ = r;
832
833 // now sort the roots
834 bubble_sort(tValues, (int)(roots - tValues));
835#endif
836 }
837 else // we have 1 real root
838 {
839 SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3));
840 A = SkFPCubeRoot(A);
841 if (SkFPGT(R, 0))
842 A = SkFPNeg(A);
843
844 if (A != 0)
845 A = SkFPAdd(A, SkFPDiv(Q, A));
846 r = SkFPToScalar(SkFPSub(A, adiv3));
847 if (is_unit_interval(r))
848 *roots++ = r;
849 }
850
851 return (int)(roots - tValues);
852}
853
854/* Looking for F' dot F'' == 0
855
856 A = b - a
857 B = c - 2b + a
858 C = d - 3c + 3b - a
859
860 F' = 3Ct^2 + 6Bt + 3A
861 F'' = 6Ct + 6B
862
863 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
864*/
865static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4])
866{
867 SkScalar a = src[2] - src[0];
868 SkScalar b = src[4] - 2 * src[2] + src[0];
869 SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0];
870
871 SkFP A = SkScalarToFP(a);
872 SkFP B = SkScalarToFP(b);
873 SkFP C = SkScalarToFP(c);
874
875 coeff[0] = SkFPMul(C, C);
876 coeff[1] = SkFPMulInt(SkFPMul(B, C), 3);
877 coeff[2] = SkFPMulInt(SkFPMul(B, B), 2);
878 coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A));
879 coeff[3] = SkFPMul(A, B);
880}
881
882// EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1
883//#define kMinTValueForChopping (SK_Scalar1 / 256)
884#define kMinTValueForChopping 0
885
886/* Looking for F' dot F'' == 0
887
888 A = b - a
889 B = c - 2b + a
890 C = d - 3c + 3b - a
891
892 F' = 3Ct^2 + 6Bt + 3A
893 F'' = 6Ct + 6B
894
895 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
896*/
897int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3])
898{
899 SkFP coeffX[4], coeffY[4];
900 int i;
901
902 formulate_F1DotF2(&src[0].fX, coeffX);
903 formulate_F1DotF2(&src[0].fY, coeffY);
904
905 for (i = 0; i < 4; i++)
906 coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]);
907
908 SkScalar t[3];
909 int count = solve_cubic_polynomial(coeffX, t);
910 int maxCount = 0;
911
912 // now remove extrema where the curvature is zero (mins)
913 // !!!! need a test for this !!!!
914 for (i = 0; i < count; i++)
915 {
916 // if (not_min_curvature())
917 if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping)
918 tValues[maxCount++] = t[i];
919 }
920 return maxCount;
921}
922
923int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3])
924{
925 SkScalar t_storage[3];
926
927 if (tValues == NULL)
928 tValues = t_storage;
929
930 int count = SkFindCubicMaxCurvature(src, tValues);
931
932 if (dst)
933 {
934 if (count == 0)
935 memcpy(dst, src, 4 * sizeof(SkPoint));
936 else
937 SkChopCubicAt(src, dst, tValues, count);
938 }
939 return count + 1;
940}
941
942////////////////////////////////////////////////////////////////////////////////
943
944/* Find t value for quadratic [a, b, c] = d.
945 Return 0 if there is no solution within [0, 1)
946*/
947static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d)
948{
949 // At^2 + Bt + C = d
950 SkScalar A = a - 2 * b + c;
951 SkScalar B = 2 * (b - a);
952 SkScalar C = a - d;
953
954 SkScalar roots[2];
955 int count = SkFindUnitQuadRoots(A, B, C, roots);
956
957 SkASSERT(count <= 1);
958 return count == 1 ? roots[0] : 0;
959}
960
961/* given a quad-curve and a point (x,y), chop the quad at that point and return
962 the new quad's offCurve point. Should only return false if the computed pos
963 is the start of the curve (i.e. root == 0)
964*/
965static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve)
966{
967 const SkScalar* base;
968 SkScalar value;
969
970 if (SkScalarAbs(x) < SkScalarAbs(y)) {
971 base = &quad[0].fX;
972 value = x;
973 } else {
974 base = &quad[0].fY;
975 value = y;
976 }
977
978 // note: this returns 0 if it thinks value is out of range, meaning the
979 // root might return something outside of [0, 1)
980 SkScalar t = quad_solve(base[0], base[2], base[4], value);
981
982 if (t > 0)
983 {
984 SkPoint tmp[5];
985 SkChopQuadAt(quad, tmp, t);
986 *offCurve = tmp[1];
987 return true;
988 } else {
989 /* t == 0 means either the value triggered a root outside of [0, 1)
990 For our purposes, we can ignore the <= 0 roots, but we want to
991 catch the >= 1 roots (which given our caller, will basically mean
992 a root of 1, give-or-take numerical instability). If we are in the
993 >= 1 case, return the existing offCurve point.
994
995 The test below checks to see if we are close to the "end" of the
996 curve (near base[4]). Rather than specifying a tolerance, I just
997 check to see if value is on to the right/left of the middle point
998 (depending on the direction/sign of the end points).
999 */
1000 if ((base[0] < base[4] && value > base[2]) ||
1001 (base[0] > base[4] && value < base[2])) // should root have been 1
1002 {
1003 *offCurve = quad[1];
1004 return true;
1005 }
1006 }
1007 return false;
1008}
1009
1010static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = {
1011 { SK_Scalar1, 0 },
1012 { SK_Scalar1, SK_ScalarTanPIOver8 },
1013 { SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 },
1014 { SK_ScalarTanPIOver8, SK_Scalar1 },
1015
1016 { 0, SK_Scalar1 },
1017 { -SK_ScalarTanPIOver8, SK_Scalar1 },
1018 { -SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 },
1019 { -SK_Scalar1, SK_ScalarTanPIOver8 },
1020
1021 { -SK_Scalar1, 0 },
1022 { -SK_Scalar1, -SK_ScalarTanPIOver8 },
1023 { -SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 },
1024 { -SK_ScalarTanPIOver8, -SK_Scalar1 },
1025
1026 { 0, -SK_Scalar1 },
1027 { SK_ScalarTanPIOver8, -SK_Scalar1 },
1028 { SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 },
1029 { SK_Scalar1, -SK_ScalarTanPIOver8 },
1030
1031 { SK_Scalar1, 0 }
1032};
1033
1034int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop,
1035 SkRotationDirection dir, const SkMatrix* userMatrix,
1036 SkPoint quadPoints[])
1037{
1038 // rotate by x,y so that uStart is (1.0)
1039 SkScalar x = SkPoint::DotProduct(uStart, uStop);
1040 SkScalar y = SkPoint::CrossProduct(uStart, uStop);
1041
1042 SkScalar absX = SkScalarAbs(x);
1043 SkScalar absY = SkScalarAbs(y);
1044
1045 int pointCount;
1046
1047 // check for (effectively) coincident vectors
1048 // this can happen if our angle is nearly 0 or nearly 180 (y == 0)
1049 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0)
1050 if (absY <= SK_ScalarNearlyZero && x > 0 &&
1051 ((y >= 0 && kCW_SkRotationDirection == dir) ||
1052 (y <= 0 && kCCW_SkRotationDirection == dir))) {
1053
1054 // just return the start-point
1055 quadPoints[0].set(SK_Scalar1, 0);
1056 pointCount = 1;
1057 } else {
1058 if (dir == kCCW_SkRotationDirection)
1059 y = -y;
1060
1061 // what octant (quadratic curve) is [xy] in?
1062 int oct = 0;
1063 bool sameSign = true;
1064
1065 if (0 == y)
1066 {
1067 oct = 4; // 180
1068 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero);
1069 }
1070 else if (0 == x)
1071 {
1072 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero);
1073 if (y > 0)
1074 oct = 2; // 90
1075 else
1076 oct = 6; // 270
1077 }
1078 else
1079 {
1080 if (y < 0)
1081 oct += 4;
1082 if ((x < 0) != (y < 0))
1083 {
1084 oct += 2;
1085 sameSign = false;
1086 }
1087 if ((absX < absY) == sameSign)
1088 oct += 1;
1089 }
1090
1091 int wholeCount = oct << 1;
1092 memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint));
1093
1094 const SkPoint* arc = &gQuadCirclePts[wholeCount];
1095 if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1]))
1096 {
1097 quadPoints[wholeCount + 2].set(x, y);
1098 wholeCount += 2;
1099 }
1100 pointCount = wholeCount + 1;
1101 }
1102
1103 // now handle counter-clockwise and the initial unitStart rotation
1104 SkMatrix matrix;
1105 matrix.setSinCos(uStart.fY, uStart.fX);
1106 if (dir == kCCW_SkRotationDirection) {
1107 matrix.preScale(SK_Scalar1, -SK_Scalar1);
1108 }
1109 if (userMatrix) {
1110 matrix.postConcat(*userMatrix);
1111 }
1112 matrix.mapPoints(quadPoints, pointCount);
1113 return pointCount;
1114}
1115