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
17static 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
44static 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
69static 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
106ArcLengthIntegrator::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.
157SkScalar 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
172SkCurveMeasure::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
204SkScalar 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
220SkScalar 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
306void 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