|
9 | 9 | #include "SkNx.h" |
10 | 10 | #include "../../src/pathops/SkPathOpsCubic.h" |
11 | 11 |
|
12 | | -void SkCubicMap::setPts(SkPoint p1, SkPoint p2) { |
13 | | - Sk2s s1 = Sk2s::Load(&p1) * 3; |
14 | | - Sk2s s2 = Sk2s::Load(&p2) * 3; |
15 | | - |
16 | | - s1 = Sk2s::Min(Sk2s::Max(s1, 0), 3); |
17 | | - s2 = Sk2s::Min(Sk2s::Max(s2, 0), 3); |
18 | | - |
19 | | - (Sk2s(1) + s1 - s2).store(&fCoeff[0]); |
20 | | - (s2 - s1 - s1).store(&fCoeff[1]); |
21 | | - s1.store(&fCoeff[2]); |
| 12 | +static float eval_poly3(float a, float b, float c, float d, float t) { |
| 13 | + return ((a * t + b) * t + c) * t + d; |
| 14 | +} |
22 | 15 |
|
23 | | - this->buildXTable(); |
| 16 | +static float eval_poly2(float a, float b, float c, float t) { |
| 17 | + return (a * t + b) * t + c; |
24 | 18 | } |
25 | 19 |
|
26 | | -SkPoint SkCubicMap::computeFromT(float t) const { |
27 | | - Sk2s a = Sk2s::Load(&fCoeff[0]); |
28 | | - Sk2s b = Sk2s::Load(&fCoeff[1]); |
29 | | - Sk2s c = Sk2s::Load(&fCoeff[2]); |
| 20 | +static float eval_poly1(float a, float b, float t) { |
| 21 | + return a * t + b; |
| 22 | +} |
30 | 23 |
|
31 | | - SkPoint result; |
32 | | - (((a * t + b) * t + c) * t).store(&result); |
33 | | - return result; |
| 24 | +static float guess_nice_cubic_root(float A, float B, float C, float D) { |
| 25 | + return -D; |
34 | 26 | } |
35 | 27 |
|
36 | | -float SkCubicMap::computeYFromX(float x) const { |
37 | | - x = SkTPin<float>(x, 0, 0.99999f) * kTableCount; |
38 | | - float ix = sk_float_floor(x); |
39 | | - int index = (int)ix; |
40 | | - SkASSERT((unsigned)index < SK_ARRAY_COUNT(fXTable)); |
41 | | - return this->computeFromT(fXTable[index].fT0 + fXTable[index].fDT * (x - ix)).fY; |
| 28 | +#ifdef SK_DEBUG |
| 29 | +static bool valid(float r) { |
| 30 | + return r >= 0 && r <= 1; |
42 | 31 | } |
| 32 | +#endif |
43 | 33 |
|
44 | | -float SkCubicMap::hackYFromX(float x) const { |
45 | | - x = SkTPin<float>(x, 0, 0.99999f) * kTableCount; |
46 | | - float ix = sk_float_floor(x); |
47 | | - int index = (int)ix; |
48 | | - SkASSERT((unsigned)index < SK_ARRAY_COUNT(fXTable)); |
49 | | - return fXTable[index].fY0 + fXTable[index].fDY * (x - ix); |
| 34 | +/* |
| 35 | + * TODO: will this be faster if we algebraically compute the polynomials for the numer and denom |
| 36 | + * rather than compute them in parts? |
| 37 | + * |
| 38 | + * TODO: investigate Householder's method, to see if we can get away with even fewer |
| 39 | + * iterations (divides) |
| 40 | + */ |
| 41 | +static float solve_nice_cubic_halley(float A, float B, float C, float D) { |
| 42 | + const int MAX_ITERS = 3; // 3 is accurate to 0.0000112 (anecdotally) |
| 43 | + // 2 is accurate to 0.0188965 (anecdotally) |
| 44 | + const float A3 = 3 * A; |
| 45 | + const float B2 = B + B; |
| 46 | + |
| 47 | + float t = guess_nice_cubic_root(A, B, C, D); |
| 48 | + for (int iters = 0; iters < MAX_ITERS; ++iters) { |
| 49 | + float f = eval_poly3(A, B, C, D, t); // f = At^3 + Bt^2 + Ct + D |
| 50 | + float fp = eval_poly2(A3, B2, C, t); // f' = 3At^2 + 2Bt + C |
| 51 | + float fpp = eval_poly1(A3 + A3, B2, t); // f'' = 6At + 2B |
| 52 | + |
| 53 | + float numer = 2 * fp * f; |
| 54 | + float denom = 2 * fp * fp - f * fpp; |
| 55 | + float delta = numer / denom; |
| 56 | + float new_t = t - delta; |
| 57 | + SkASSERT(valid(new_t)); |
| 58 | + t = new_t; |
| 59 | + } |
| 60 | + SkASSERT(valid(t)); |
| 61 | + return t; |
50 | 62 | } |
51 | 63 |
|
52 | | -static float compute_t_from_x(float A, float B, float C, float x) { |
| 64 | +#ifdef SK_DEBUG |
| 65 | +static float compute_slow(float A, float B, float C, float x) { |
53 | 66 | double roots[3]; |
54 | 67 | SkDEBUGCODE(int count =) SkDCubic::RootsValidT(A, B, C, -x, roots); |
55 | 68 | SkASSERT(count == 1); |
56 | 69 | return (float)roots[0]; |
57 | 70 | } |
58 | 71 |
|
59 | | -void SkCubicMap::buildXTable() { |
60 | | - float prevT = 0; |
| 72 | +static float max_err; |
| 73 | +#endif |
61 | 74 |
|
62 | | - const float dx = 1.0f / kTableCount; |
63 | | - float x = dx; |
| 75 | +static float compute_t_from_x(float A, float B, float C, float x) { |
| 76 | +#ifdef SK_DEBUG |
| 77 | + float answer = compute_slow(A, B, C, x); |
| 78 | +#endif |
| 79 | + float answer2 = solve_nice_cubic_halley(A, B, C, -x); |
| 80 | +#ifdef SK_DEBUG |
| 81 | + float err = sk_float_abs(answer - answer2); |
| 82 | + if (err > max_err) { |
| 83 | + max_err = err; |
| 84 | + SkDebugf("max error %g\n", max_err); |
| 85 | + } |
| 86 | +#endif |
| 87 | + return answer2; |
| 88 | +} |
64 | 89 |
|
65 | | - fXTable[0].fT0 = 0; |
66 | | - fXTable[0].fY0 = 0; |
67 | | - for (int i = 1; i < kTableCount; ++i) { |
68 | | - float t = compute_t_from_x(fCoeff[0].fX, fCoeff[1].fX, fCoeff[2].fX, x); |
69 | | - SkASSERT(t > prevT); |
| 90 | +float SkCubicMap::computeYFromX(float x) const { |
| 91 | + SkASSERT(valid(x)); |
| 92 | + float t = compute_t_from_x(fCoeff[0].fX, fCoeff[1].fX, fCoeff[2].fX, x); |
| 93 | + float a = fCoeff[0].fY; |
| 94 | + float b = fCoeff[1].fY; |
| 95 | + float c = fCoeff[2].fY; |
| 96 | + return ((a * t + b) * t + c) * t; |
| 97 | +} |
70 | 98 |
|
71 | | - fXTable[i - 1].fDT = t - prevT; |
72 | | - fXTable[i].fT0 = t; |
| 99 | +void SkCubicMap::setPts(SkPoint p1, SkPoint p2) { |
| 100 | + Sk2s s1 = Sk2s::Load(&p1) * 3; |
| 101 | + Sk2s s2 = Sk2s::Load(&p2) * 3; |
73 | 102 |
|
74 | | - SkPoint p = this->computeFromT(t); |
75 | | - fXTable[i - 1].fDY = p.fY - fXTable[i - 1].fY0; |
76 | | - fXTable[i].fY0 = p.fY; |
| 103 | + s1 = Sk2s::Min(Sk2s::Max(s1, 0), 3); |
| 104 | + s2 = Sk2s::Min(Sk2s::Max(s2, 0), 3); |
77 | 105 |
|
78 | | - prevT = t; |
79 | | - x += dx; |
80 | | - } |
81 | | - fXTable[kTableCount - 1].fDT = 1 - prevT; |
82 | | - fXTable[kTableCount - 1].fDY = 1 - fXTable[kTableCount - 1].fY0; |
| 106 | + (Sk2s(1) + s1 - s2).store(&fCoeff[0]); |
| 107 | + (s2 - s1 - s1).store(&fCoeff[1]); |
| 108 | + s1.store(&fCoeff[2]); |
| 109 | +} |
| 110 | + |
| 111 | +SkPoint SkCubicMap::computeFromT(float t) const { |
| 112 | + Sk2s a = Sk2s::Load(&fCoeff[0]); |
| 113 | + Sk2s b = Sk2s::Load(&fCoeff[1]); |
| 114 | + Sk2s c = Sk2s::Load(&fCoeff[2]); |
| 115 | + |
| 116 | + SkPoint result; |
| 117 | + (((a * t + b) * t + c) * t).store(&result); |
| 118 | + return result; |
83 | 119 | } |
0 commit comments