1 /*
2 * Copyright 2018 Google Inc.
3 *
4 * Use of this source code is governed by a BSD-style license that can be
5 * found in the LICENSE file.
6 */
7
8 #include "SkCubicMap.h"
9 #include "SkNx.h"
10
11 //#define CUBICMAP_TRACK_MAX_ERROR
12
13 #ifdef CUBICMAP_TRACK_MAX_ERROR
14 #include "../../src/pathops/SkPathOpsCubic.h"
15 #endif
16
eval_poly3(float a,float b,float c,float d,float t)17 static float eval_poly3(float a, float b, float c, float d, float t) {
18 return ((a * t + b) * t + c) * t + d;
19 }
20
eval_poly2(float a,float b,float c,float t)21 static float eval_poly2(float a, float b, float c, float t) {
22 return (a * t + b) * t + c;
23 }
24
eval_poly1(float a,float b,float t)25 static float eval_poly1(float a, float b, float t) {
26 return a * t + b;
27 }
28
guess_nice_cubic_root(float A,float B,float C,float D)29 static float guess_nice_cubic_root(float A, float B, float C, float D) {
30 return -D;
31 }
32
33 #ifdef SK_DEBUG
valid(float r)34 static bool valid(float r) {
35 return r >= 0 && r <= 1;
36 }
37 #endif
38
nearly_zero(SkScalar x)39 static inline bool nearly_zero(SkScalar x) {
40 SkASSERT(x >= 0);
41 return x <= 0.0000000001f;
42 }
43
delta_nearly_zero(float delta)44 static inline bool delta_nearly_zero(float delta) {
45 return sk_float_abs(delta) <= 0.0001f;
46 }
47
48 #ifdef CUBICMAP_TRACK_MAX_ERROR
49 static int max_iters;
50 #endif
51
52 /*
53 * TODO: will this be faster if we algebraically compute the polynomials for the numer and denom
54 * rather than compute them in parts?
55 */
solve_nice_cubic_halley(float A,float B,float C,float D)56 static float solve_nice_cubic_halley(float A, float B, float C, float D) {
57 const int MAX_ITERS = 8;
58 const float A3 = 3 * A;
59 const float B2 = B + B;
60
61 float t = guess_nice_cubic_root(A, B, C, D);
62 int iters = 0;
63 for (; iters < MAX_ITERS; ++iters) {
64 float f = eval_poly3(A, B, C, D, t); // f = At^3 + Bt^2 + Ct + D
65 float fp = eval_poly2(A3, B2, C, t); // f' = 3At^2 + 2Bt + C
66 float fpp = eval_poly1(A3 + A3, B2, t); // f'' = 6At + 2B
67
68 float numer = 2 * fp * f;
69 if (numer == 0) {
70 break;
71 }
72 float denom = 2 * fp * fp - f * fpp;
73 float delta = numer / denom;
74 // SkDebugf("[%d] delta %g t %g\n", iters, delta, t);
75 if (delta_nearly_zero(delta)) {
76 break;
77 }
78 float new_t = t - delta;
79 SkASSERT(valid(new_t));
80 t = new_t;
81 }
82 SkASSERT(valid(t));
83 #ifdef CUBICMAP_TRACK_MAX_ERROR
84 if (iters > max_iters) {
85 max_iters = iters;
86 SkDebugf("max_iters %d\n", max_iters);
87 }
88 #endif
89 return t;
90 }
91
92 // At the moment, this technique does not appear to be better (i.e. faster at same precision)
93 // but the code is left here (at least for a while) to document the attempt.
solve_nice_cubic_householder(float A,float B,float C,float D)94 static float solve_nice_cubic_householder(float A, float B, float C, float D) {
95 const int MAX_ITERS = 8;
96 const float A3 = 3 * A;
97 const float B2 = B + B;
98
99 float t = guess_nice_cubic_root(A, B, C, D);
100 int iters = 0;
101 for (; iters < MAX_ITERS; ++iters) {
102 float f = eval_poly3(A, B, C, D, t); // f = At^3 + Bt^2 + Ct + D
103 float fp = eval_poly2(A3, B2, C, t); // f' = 3At^2 + 2Bt + C
104 float fpp = eval_poly1(A3 + A3, B2, t); // f'' = 6At + 2B
105 float fppp = A3 + A3; // f''' = 6A
106
107 float f2 = f * f;
108 float fp2 = fp * fp;
109
110 // float numer = 6 * f * fp * fp - 3 * f * f * fpp;
111 // float denom = 6 * fp * fp * fp - 6 * f * fp * fpp + f * f * fppp;
112
113 float numer = 6 * f * fp2 - 3 * f2 * fpp;
114 if (numer == 0) {
115 break;
116 }
117 float denom = 6 * (fp2 * fp - f * fp * fpp) + f2 * fppp;
118 float delta = numer / denom;
119 // SkDebugf("[%d] delta %g t %g\n", iters, delta, t);
120 if (delta_nearly_zero(delta)) {
121 break;
122 }
123 float new_t = t - delta;
124 SkASSERT(valid(new_t));
125 t = new_t;
126 }
127 SkASSERT(valid(t));
128 #ifdef CUBICMAP_TRACK_MAX_ERROR
129 if (iters > max_iters) {
130 max_iters = iters;
131 SkDebugf("max_iters %d\n", max_iters);
132 }
133 #endif
134 return t;
135 }
136
137 #ifdef CUBICMAP_TRACK_MAX_ERROR
compute_slow(float A,float B,float C,float x)138 static float compute_slow(float A, float B, float C, float x) {
139 double roots[3];
140 SkDEBUGCODE(int count =) SkDCubic::RootsValidT(A, B, C, -x, roots);
141 SkASSERT(count == 1);
142 return (float)roots[0];
143 }
144
145 static float max_err;
146 #endif
147
compute_t_from_x(float A,float B,float C,float x)148 static float compute_t_from_x(float A, float B, float C, float x) {
149 #ifdef CUBICMAP_TRACK_MAX_ERROR
150 float answer = compute_slow(A, B, C, x);
151 #endif
152 float answer2 = true ?
153 solve_nice_cubic_halley(A, B, C, -x) :
154 solve_nice_cubic_householder(A, B, C, -x);
155 #ifdef CUBICMAP_TRACK_MAX_ERROR
156 float err = sk_float_abs(answer - answer2);
157 if (err > max_err) {
158 max_err = err;
159 SkDebugf("max error %g\n", max_err);
160 }
161 #endif
162 return answer2;
163 }
164
computeYFromX(float x) const165 float SkCubicMap::computeYFromX(float x) const {
166 x = SkScalarPin(x, 0, 1);
167
168 if (nearly_zero(x) || nearly_zero(1 - x)) {
169 return x;
170 }
171 if (fType == kLine_Type) {
172 return x;
173 }
174 float t;
175 if (fType == kCubeRoot_Type) {
176 t = sk_float_pow(x / fCoeff[0].fX, 1.0f / 3);
177 } else {
178 t = compute_t_from_x(fCoeff[0].fX, fCoeff[1].fX, fCoeff[2].fX, x);
179 }
180 float a = fCoeff[0].fY;
181 float b = fCoeff[1].fY;
182 float c = fCoeff[2].fY;
183 float y = ((a * t + b) * t + c) * t;
184 SkASSERT(y >= 0);
185 return std::min(y, 1.0f);
186 }
187
coeff_nearly_zero(float delta)188 static inline bool coeff_nearly_zero(float delta) {
189 return sk_float_abs(delta) <= 0.0000001f;
190 }
191
setPts(SkPoint p1,SkPoint p2)192 void SkCubicMap::setPts(SkPoint p1, SkPoint p2) {
193 Sk2s s1 = Sk2s::Load(&p1) * 3;
194 Sk2s s2 = Sk2s::Load(&p2) * 3;
195
196 s1 = Sk2s::Min(Sk2s::Max(s1, 0), 3);
197 s2 = Sk2s::Min(Sk2s::Max(s2, 0), 3);
198
199 (Sk2s(1) + s1 - s2).store(&fCoeff[0]);
200 (s2 - s1 - s1).store(&fCoeff[1]);
201 s1.store(&fCoeff[2]);
202
203 fType = kSolver_Type;
204 if (SkScalarNearlyEqual(p1.fX, p1.fY) && SkScalarNearlyEqual(p2.fX, p2.fY)) {
205 fType = kLine_Type;
206 } else if (coeff_nearly_zero(fCoeff[1].fX) && coeff_nearly_zero(fCoeff[2].fX)) {
207 fType = kCubeRoot_Type;
208 }
209 }
210
computeFromT(float t) const211 SkPoint SkCubicMap::computeFromT(float t) const {
212 Sk2s a = Sk2s::Load(&fCoeff[0]);
213 Sk2s b = Sk2s::Load(&fCoeff[1]);
214 Sk2s c = Sk2s::Load(&fCoeff[2]);
215
216 SkPoint result;
217 (((a * t + b) * t + c) * t).store(&result);
218 return result;
219 }
220