1 /*
2 * Copyright 2016 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 "SkCurveMeasure.h"
9 #include "SkGeometry.h"
10
11 // for abs
12 #include <cmath>
13
14 #define UNIMPLEMENTED SkDEBUGF(("%s:%d unimplemented\n", __FILE__, __LINE__))
15
16 /// Used inside SkCurveMeasure::getTime's Newton's iteration
evaluate(const SkPoint pts[4],SkSegType segType,SkScalar t)17 static inline SkPoint evaluate(const SkPoint pts[4], SkSegType segType,
18 SkScalar t) {
19 SkPoint pos;
20 switch (segType) {
21 case kQuad_SegType:
22 pos = SkEvalQuadAt(pts, t);
23 break;
24 case kLine_SegType:
25 pos = SkPoint::Make(SkScalarInterp(pts[0].x(), pts[1].x(), t),
26 SkScalarInterp(pts[0].y(), pts[1].y(), t));
27 break;
28 case kCubic_SegType:
29 SkEvalCubicAt(pts, t, &pos, nullptr, nullptr);
30 break;
31 case kConic_SegType: {
32 SkConic conic(pts, pts[3].x());
33 conic.evalAt(t, &pos);
34 }
35 break;
36 default:
37 UNIMPLEMENTED;
38 }
39
40 return pos;
41 }
42
43 /// Used inside SkCurveMeasure::getTime's Newton's iteration
evaluateDerivative(const SkPoint pts[4],SkSegType segType,SkScalar t)44 static inline SkVector evaluateDerivative(const SkPoint pts[4],
45 SkSegType segType, SkScalar t) {
46 SkVector tan;
47 switch (segType) {
48 case kQuad_SegType:
49 tan = SkEvalQuadTangentAt(pts, t);
50 break;
51 case kLine_SegType:
52 tan = pts[1] - pts[0];
53 break;
54 case kCubic_SegType:
55 SkEvalCubicAt(pts, t, nullptr, &tan, nullptr);
56 break;
57 case kConic_SegType: {
58 SkConic conic(pts, pts[3].x());
59 conic.evalAt(t, nullptr, &tan);
60 }
61 break;
62 default:
63 UNIMPLEMENTED;
64 }
65
66 return tan;
67 }
68 /// Used in ArcLengthIntegrator::computeLength
evaluateDerivativeLength(const Sk8f & ts,const float (& xCoeff)[3][8],const float (& yCoeff)[3][8],const SkSegType segType)69 static inline Sk8f evaluateDerivativeLength(const Sk8f& ts,
70 const float (&xCoeff)[3][8],
71 const float (&yCoeff)[3][8],
72 const SkSegType segType) {
73 Sk8f x;
74 Sk8f y;
75
76 Sk8f x0 = Sk8f::Load(&xCoeff[0]),
77 x1 = Sk8f::Load(&xCoeff[1]),
78 x2 = Sk8f::Load(&xCoeff[2]);
79
80 Sk8f y0 = Sk8f::Load(&yCoeff[0]),
81 y1 = Sk8f::Load(&yCoeff[1]),
82 y2 = Sk8f::Load(&yCoeff[2]);
83
84 switch (segType) {
85 case kQuad_SegType:
86 x = x0*ts + x1;
87 y = y0*ts + y1;
88 break;
89 case kCubic_SegType:
90 x = (x0*ts + x1)*ts + x2;
91 y = (y0*ts + y1)*ts + y2;
92 break;
93 case kConic_SegType:
94 UNIMPLEMENTED;
95 break;
96 default:
97 UNIMPLEMENTED;
98 }
99
100 x = x * x;
101 y = y * y;
102
103 return (x + y).sqrt();
104 }
105
ArcLengthIntegrator(const SkPoint * pts,SkSegType segType)106 ArcLengthIntegrator::ArcLengthIntegrator(const SkPoint* pts, SkSegType segType)
107 : fSegType(segType) {
108 switch (fSegType) {
109 case kQuad_SegType: {
110 float Ax = pts[0].x();
111 float Bx = pts[1].x();
112 float Cx = pts[2].x();
113 float Ay = pts[0].y();
114 float By = pts[1].y();
115 float Cy = pts[2].y();
116
117 // precompute coefficients for derivative
118 Sk8f(2*(Ax - 2*Bx + Cx)).store(&xCoeff[0]);
119 Sk8f(2*(Bx - Ax)).store(&xCoeff[1]);
120
121 Sk8f(2*(Ay - 2*By + Cy)).store(&yCoeff[0]);
122 Sk8f(2*(By - Ay)).store(&yCoeff[1]);
123 }
124 break;
125 case kCubic_SegType:
126 {
127 float Ax = pts[0].x();
128 float Bx = pts[1].x();
129 float Cx = pts[2].x();
130 float Dx = pts[3].x();
131 float Ay = pts[0].y();
132 float By = pts[1].y();
133 float Cy = pts[2].y();
134 float Dy = pts[3].y();
135
136 // precompute coefficients for derivative
137 Sk8f(3*(-Ax + 3*(Bx - Cx) + Dx)).store(&xCoeff[0]);
138 Sk8f(6*(Ax - 2*Bx + Cx)).store(&xCoeff[1]);
139 Sk8f(3*(-Ax + Bx)).store(&xCoeff[2]);
140
141 Sk8f(3*(-Ay + 3*(By - Cy) + Dy)).store(&yCoeff[0]);
142 Sk8f(6*(Ay - 2*By + Cy)).store(&yCoeff[1]);
143 Sk8f(3*(-Ay + By)).store(&yCoeff[2]);
144 }
145 break;
146 case kConic_SegType:
147 UNIMPLEMENTED;
148 break;
149 default:
150 UNIMPLEMENTED;
151 }
152 }
153
154 // We use Gaussian quadrature
155 // (https://en.wikipedia.org/wiki/Gaussian_quadrature)
156 // to approximate the arc length integral here, because it is amenable to SIMD.
computeLength(SkScalar t)157 SkScalar ArcLengthIntegrator::computeLength(SkScalar t) {
158 SkScalar length = 0.0f;
159
160 Sk8f lengths = evaluateDerivativeLength(absc*t, xCoeff, yCoeff, fSegType);
161 lengths = weights*lengths;
162 // is it faster or more accurate to sum and then multiply or vice versa?
163 lengths = lengths*(t*0.5f);
164
165 // Why does SkNx index with ints? does negative index mean something?
166 for (int i = 0; i < 8; i++) {
167 length += lengths[i];
168 }
169 return length;
170 }
171
SkCurveMeasure(const SkPoint * pts,SkSegType segType)172 SkCurveMeasure::SkCurveMeasure(const SkPoint* pts, SkSegType segType)
173 : fSegType(segType) {
174 switch (fSegType) {
175 case SkSegType::kQuad_SegType:
176 for (size_t i = 0; i < 3; i++) {
177 fPts[i] = pts[i];
178 }
179 break;
180 case SkSegType::kLine_SegType:
181 fPts[0] = pts[0];
182 fPts[1] = pts[1];
183 fLength = (fPts[1] - fPts[0]).length();
184 break;
185 case SkSegType::kCubic_SegType:
186 for (size_t i = 0; i < 4; i++) {
187 fPts[i] = pts[i];
188 }
189 break;
190 case SkSegType::kConic_SegType:
191 for (size_t i = 0; i < 4; i++) {
192 fPts[i] = pts[i];
193 }
194 break;
195 default:
196 UNIMPLEMENTED;
197 break;
198 }
199 if (kLine_SegType != segType) {
200 fIntegrator = ArcLengthIntegrator(fPts, fSegType);
201 }
202 }
203
getLength()204 SkScalar SkCurveMeasure::getLength() {
205 if (-1.0f == fLength) {
206 fLength = fIntegrator.computeLength(1.0f);
207 }
208 return fLength;
209 }
210
211 // Given an arc length targetLength, we want to determine what t
212 // gives us the corresponding arc length along the curve.
213 // We do this by letting the arc length integral := f(t) and
214 // solving for the root of the equation f(t) - targetLength = 0
215 // using Newton's method and lerp-bisection.
216 // The computationally expensive parts are the integral approximation
217 // at each step, and computing the derivative of the arc length integral,
218 // which is equal to the length of the tangent (so we have to do a sqrt).
219
getTime(SkScalar targetLength)220 SkScalar SkCurveMeasure::getTime(SkScalar targetLength) {
221 if (targetLength <= 0.0f) {
222 return 0.0f;
223 }
224
225 SkScalar currentLength = getLength();
226
227 if (targetLength > currentLength || (SkScalarNearlyEqual(targetLength, currentLength))) {
228 return 1.0f;
229 }
230 if (kLine_SegType == fSegType) {
231 return targetLength / currentLength;
232 }
233
234 // initial estimate of t is percentage of total length
235 SkScalar currentT = targetLength / currentLength;
236 SkScalar prevT = -1.0f;
237 SkScalar newT;
238
239 SkScalar minT = 0.0f;
240 SkScalar maxT = 1.0f;
241
242 int iterations = 0;
243 while (iterations < kNewtonIters + kBisectIters) {
244 currentLength = fIntegrator.computeLength(currentT);
245 SkScalar lengthDiff = currentLength - targetLength;
246
247 // Update root bounds.
248 // If lengthDiff is positive, we have overshot the target, so
249 // we know the current t is an upper bound, and similarly
250 // for the lower bound.
251 if (lengthDiff > 0.0f) {
252 if (currentT < maxT) {
253 maxT = currentT;
254 }
255 } else {
256 if (currentT > minT) {
257 minT = currentT;
258 }
259 }
260
261 // We have a tolerance on both the absolute value of the difference and
262 // on the t value
263 // because we may not have enough precision in the t to get close enough
264 // in the length.
265 if ((std::abs(lengthDiff) < kTolerance) ||
266 (std::abs(prevT - currentT) < kTolerance)) {
267 break;
268 }
269
270 prevT = currentT;
271 if (iterations < kNewtonIters) {
272 // This is just newton's formula.
273 SkScalar dt = evaluateDerivative(fPts, fSegType, currentT).length();
274 newT = currentT - (lengthDiff / dt);
275
276 // If newT is out of bounds, bisect inside newton.
277 if ((newT < 0.0f) || (newT > 1.0f)) {
278 newT = (minT + maxT) * 0.5f;
279 }
280 } else if (iterations < kNewtonIters + kBisectIters) {
281 if (lengthDiff > 0.0f) {
282 maxT = currentT;
283 } else {
284 minT = currentT;
285 }
286 // TODO(hstern) do a lerp here instead of a bisection
287 newT = (minT + maxT) * 0.5f;
288 } else {
289 SkDEBUGF(("%.7f %.7f didn't get close enough after bisection.\n",
290 currentT, currentLength));
291 break;
292 }
293 currentT = newT;
294
295 SkASSERT(minT <= maxT);
296
297 iterations++;
298 }
299
300 // debug. is there an SKDEBUG or something for ifdefs?
301 fIters = iterations;
302
303 return currentT;
304 }
305
getPosTanTime(SkScalar targetLength,SkPoint * pos,SkVector * tan,SkScalar * time)306 void SkCurveMeasure::getPosTanTime(SkScalar targetLength, SkPoint* pos,
307 SkVector* tan, SkScalar* time) {
308 SkScalar t = getTime(targetLength);
309
310 if (time) {
311 *time = t;
312 }
313 if (pos) {
314 *pos = evaluate(fPts, fSegType, t);
315 }
316 if (tan) {
317 *tan = evaluateDerivative(fPts, fSegType, t);
318 }
319 }
320