1/*
2 * Copyright 2012 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#include "CubicUtilities.h"
8#include "Extrema.h"
9#include "LineUtilities.h"
10#include "QuadraticUtilities.h"
11
12const int gPrecisionUnit = 256; // FIXME: arbitrary -- should try different values in test framework
13
14// FIXME: cache keep the bounds and/or precision with the caller?
15double calcPrecision(const Cubic& cubic) {
16    _Rect dRect;
17    dRect.setBounds(cubic); // OPTIMIZATION: just use setRawBounds ?
18    double width = dRect.right - dRect.left;
19    double height = dRect.bottom - dRect.top;
20    return (width > height ? width : height) / gPrecisionUnit;
21}
22
23#if SK_DEBUG
24double calcPrecision(const Cubic& cubic, double t, double scale) {
25    Cubic part;
26    sub_divide(cubic, SkTMax(0., t - scale), SkTMin(1., t + scale), part);
27    return calcPrecision(part);
28}
29#endif
30
31bool clockwise(const Cubic& c) {
32    double sum = (c[0].x - c[3].x) * (c[0].y + c[3].y);
33    for (int idx = 0; idx < 3; ++idx){
34        sum += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y);
35    }
36    return sum <= 0;
37}
38
39void coefficients(const double* cubic, double& A, double& B, double& C, double& D) {
40    A = cubic[6]; // d
41    B = cubic[4] * 3; // 3*c
42    C = cubic[2] * 3; // 3*b
43    D = cubic[0]; // a
44    A -= D - C + B;     // A =   -a + 3*b - 3*c + d
45    B += 3 * D - 2 * C; // B =  3*a - 6*b + 3*c
46    C -= 3 * D;         // C = -3*a + 3*b
47}
48
49bool controls_contained_by_ends(const Cubic& c) {
50    _Vector startTan = c[1] - c[0];
51    if (startTan.x == 0 && startTan.y == 0) {
52        startTan = c[2] - c[0];
53    }
54    _Vector endTan = c[2] - c[3];
55    if (endTan.x == 0 && endTan.y == 0) {
56        endTan = c[1] - c[3];
57    }
58    if (startTan.dot(endTan) >= 0) {
59        return false;
60    }
61    _Line startEdge = {c[0], c[0]};
62    startEdge[1].x -= startTan.y;
63    startEdge[1].y += startTan.x;
64    _Line endEdge = {c[3], c[3]};
65    endEdge[1].x -= endTan.y;
66    endEdge[1].y += endTan.x;
67    double leftStart1 = is_left(startEdge, c[1]);
68    if (leftStart1 * is_left(startEdge, c[2]) < 0) {
69        return false;
70    }
71    double leftEnd1 = is_left(endEdge, c[1]);
72    if (leftEnd1 * is_left(endEdge, c[2]) < 0) {
73        return false;
74    }
75    return leftStart1 * leftEnd1 >= 0;
76}
77
78bool ends_are_extrema_in_x_or_y(const Cubic& c) {
79    return (between(c[0].x, c[1].x, c[3].x) && between(c[0].x, c[2].x, c[3].x))
80            || (between(c[0].y, c[1].y, c[3].y) && between(c[0].y, c[2].y, c[3].y));
81}
82
83bool monotonic_in_y(const Cubic& c) {
84    return between(c[0].y, c[1].y, c[3].y) && between(c[0].y, c[2].y, c[3].y);
85}
86
87bool serpentine(const Cubic& c) {
88    if (!controls_contained_by_ends(c)) {
89        return false;
90    }
91    double wiggle = (c[0].x - c[2].x) * (c[0].y + c[2].y);
92    for (int idx = 0; idx < 2; ++idx){
93        wiggle += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y);
94    }
95    double waggle = (c[1].x - c[3].x) * (c[1].y + c[3].y);
96    for (int idx = 1; idx < 3; ++idx){
97        waggle += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y);
98    }
99    return wiggle * waggle < 0;
100}
101
102// cubic roots
103
104const double PI = 4 * atan(1);
105
106// from SkGeometry.cpp (and Numeric Solutions, 5.6)
107int cubicRootsValidT(double A, double B, double C, double D, double t[3]) {
108#if 0
109    if (approximately_zero(A)) {  // we're just a quadratic
110        return quadraticRootsValidT(B, C, D, t);
111    }
112    double a, b, c;
113    {
114        double invA = 1 / A;
115        a = B * invA;
116        b = C * invA;
117        c = D * invA;
118    }
119    double a2 = a * a;
120    double Q = (a2 - b * 3) / 9;
121    double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
122    double Q3 = Q * Q * Q;
123    double R2MinusQ3 = R * R - Q3;
124    double adiv3 = a / 3;
125    double* roots = t;
126    double r;
127
128    if (R2MinusQ3 < 0)   // we have 3 real roots
129    {
130        double theta = acos(R / sqrt(Q3));
131        double neg2RootQ = -2 * sqrt(Q);
132
133        r = neg2RootQ * cos(theta / 3) - adiv3;
134        if (is_unit_interval(r))
135            *roots++ = r;
136
137        r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
138        if (is_unit_interval(r))
139            *roots++ = r;
140
141        r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
142        if (is_unit_interval(r))
143            *roots++ = r;
144    }
145    else                // we have 1 real root
146    {
147        double A = fabs(R) + sqrt(R2MinusQ3);
148        A = cube_root(A);
149        if (R > 0) {
150            A = -A;
151        }
152        if (A != 0) {
153            A += Q / A;
154        }
155        r = A - adiv3;
156        if (is_unit_interval(r))
157            *roots++ = r;
158    }
159    return (int)(roots - t);
160#else
161    double s[3];
162    int realRoots = cubicRootsReal(A, B, C, D, s);
163    int foundRoots = add_valid_ts(s, realRoots, t);
164    return foundRoots;
165#endif
166}
167
168int cubicRootsReal(double A, double B, double C, double D, double s[3]) {
169#if SK_DEBUG
170    // create a string mathematica understands
171    // GDB set print repe 15 # if repeated digits is a bother
172    //     set print elements 400 # if line doesn't fit
173    char str[1024];
174    bzero(str, sizeof(str));
175    sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A, B, C, D);
176    mathematica_ize(str, sizeof(str));
177#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA
178    SkDebugf("%s\n", str);
179#endif
180#endif
181    if (approximately_zero(A)
182            && approximately_zero_when_compared_to(A, B)
183            && approximately_zero_when_compared_to(A, C)
184            && approximately_zero_when_compared_to(A, D)) {  // we're just a quadratic
185        return quadraticRootsReal(B, C, D, s);
186    }
187    if (approximately_zero_when_compared_to(D, A)
188            && approximately_zero_when_compared_to(D, B)
189            && approximately_zero_when_compared_to(D, C)) { // 0 is one root
190        int num = quadraticRootsReal(A, B, C, s);
191        for (int i = 0; i < num; ++i) {
192            if (approximately_zero(s[i])) {
193                return num;
194            }
195        }
196        s[num++] = 0;
197        return num;
198    }
199    if (approximately_zero(A + B + C + D)) { // 1 is one root
200        int num = quadraticRootsReal(A, A + B, -D, s);
201        for (int i = 0; i < num; ++i) {
202            if (AlmostEqualUlps(s[i], 1)) {
203                return num;
204            }
205        }
206        s[num++] = 1;
207        return num;
208    }
209    double a, b, c;
210    {
211        double invA = 1 / A;
212        a = B * invA;
213        b = C * invA;
214        c = D * invA;
215    }
216    double a2 = a * a;
217    double Q = (a2 - b * 3) / 9;
218    double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
219    double R2 = R * R;
220    double Q3 = Q * Q * Q;
221    double R2MinusQ3 = R2 - Q3;
222    double adiv3 = a / 3;
223    double r;
224    double* roots = s;
225#if 0
226    if (approximately_zero_squared(R2MinusQ3) && AlmostEqualUlps(R2, Q3)) {
227        if (approximately_zero_squared(R)) {/* one triple solution */
228            *roots++ = -adiv3;
229        } else { /* one single and one double solution */
230
231            double u = cube_root(-R);
232            *roots++ = 2 * u - adiv3;
233            *roots++ = -u - adiv3;
234        }
235    }
236    else
237#endif
238    if (R2MinusQ3 < 0)   // we have 3 real roots
239    {
240        double theta = acos(R / sqrt(Q3));
241        double neg2RootQ = -2 * sqrt(Q);
242
243        r = neg2RootQ * cos(theta / 3) - adiv3;
244        *roots++ = r;
245
246        r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
247        if (!AlmostEqualUlps(s[0], r)) {
248            *roots++ = r;
249        }
250        r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
251        if (!AlmostEqualUlps(s[0], r) && (roots - s == 1 || !AlmostEqualUlps(s[1], r))) {
252            *roots++ = r;
253        }
254    }
255    else                // we have 1 real root
256    {
257        double sqrtR2MinusQ3 = sqrt(R2MinusQ3);
258        double A = fabs(R) + sqrtR2MinusQ3;
259        A = cube_root(A);
260        if (R > 0) {
261            A = -A;
262        }
263        if (A != 0) {
264            A += Q / A;
265        }
266        r = A - adiv3;
267        *roots++ = r;
268        if (AlmostEqualUlps(R2, Q3)) {
269            r = -A / 2 - adiv3;
270            if (!AlmostEqualUlps(s[0], r)) {
271                *roots++ = r;
272            }
273        }
274    }
275    return (int)(roots - s);
276}
277
278// from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf
279// c(t)  = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3
280// c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2
281//       = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2
282static double derivativeAtT(const double* cubic, double t) {
283    double one_t = 1 - t;
284    double a = cubic[0];
285    double b = cubic[2];
286    double c = cubic[4];
287    double d = cubic[6];
288    return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t);
289}
290
291double dx_at_t(const Cubic& cubic, double t) {
292    return derivativeAtT(&cubic[0].x, t);
293}
294
295double dy_at_t(const Cubic& cubic, double t) {
296    return derivativeAtT(&cubic[0].y, t);
297}
298
299// OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t?
300_Vector dxdy_at_t(const Cubic& cubic, double t) {
301    _Vector result = { derivativeAtT(&cubic[0].x, t), derivativeAtT(&cubic[0].y, t) };
302    return result;
303}
304
305// OPTIMIZE? share code with formulate_F1DotF2
306int find_cubic_inflections(const Cubic& src, double tValues[])
307{
308    double Ax = src[1].x - src[0].x;
309    double Ay = src[1].y - src[0].y;
310    double Bx = src[2].x - 2 * src[1].x + src[0].x;
311    double By = src[2].y - 2 * src[1].y + src[0].y;
312    double Cx = src[3].x + 3 * (src[1].x - src[2].x) - src[0].x;
313    double Cy = src[3].y + 3 * (src[1].y - src[2].y) - src[0].y;
314    return quadraticRootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues);
315}
316
317static void formulate_F1DotF2(const double src[], double coeff[4])
318{
319    double a = src[2] - src[0];
320    double b = src[4] - 2 * src[2] + src[0];
321    double c = src[6] + 3 * (src[2] - src[4]) - src[0];
322    coeff[0] = c * c;
323    coeff[1] = 3 * b * c;
324    coeff[2] = 2 * b * b + c * a;
325    coeff[3] = a * b;
326}
327
328/*  from SkGeometry.cpp
329    Looking for F' dot F'' == 0
330
331    A = b - a
332    B = c - 2b + a
333    C = d - 3c + 3b - a
334
335    F' = 3Ct^2 + 6Bt + 3A
336    F'' = 6Ct + 6B
337
338    F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
339*/
340int find_cubic_max_curvature(const Cubic& src, double tValues[])
341{
342    double coeffX[4], coeffY[4];
343    int i;
344    formulate_F1DotF2(&src[0].x, coeffX);
345    formulate_F1DotF2(&src[0].y, coeffY);
346    for (i = 0; i < 4; i++) {
347        coeffX[i] = coeffX[i] + coeffY[i];
348    }
349    return cubicRootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues);
350}
351
352
353bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath) {
354    double dy = cubic[index].y - cubic[zero].y;
355    double dx = cubic[index].x - cubic[zero].x;
356    if (approximately_zero(dy)) {
357        if (approximately_zero(dx)) {
358            return false;
359        }
360        memcpy(rotPath, cubic, sizeof(Cubic));
361        return true;
362    }
363    for (int index = 0; index < 4; ++index) {
364        rotPath[index].x = cubic[index].x * dx + cubic[index].y * dy;
365        rotPath[index].y = cubic[index].y * dx - cubic[index].x * dy;
366    }
367    return true;
368}
369
370#if 0 // unused for now
371double secondDerivativeAtT(const double* cubic, double t) {
372    double a = cubic[0];
373    double b = cubic[2];
374    double c = cubic[4];
375    double d = cubic[6];
376    return (c - 2 * b + a) * (1 - t) + (d - 2 * c + b) * t;
377}
378#endif
379
380_Point top(const Cubic& cubic, double startT, double endT) {
381    Cubic sub;
382    sub_divide(cubic, startT, endT, sub);
383    _Point topPt = sub[0];
384    if (topPt.y > sub[3].y || (topPt.y == sub[3].y && topPt.x > sub[3].x)) {
385        topPt = sub[3];
386    }
387    double extremeTs[2];
388    if (!monotonic_in_y(sub)) {
389        int roots = findExtrema(sub[0].y, sub[1].y, sub[2].y, sub[3].y, extremeTs);
390        for (int index = 0; index < roots; ++index) {
391            _Point mid;
392            double t = startT + (endT - startT) * extremeTs[index];
393            xy_at_t(cubic, t, mid.x, mid.y);
394            if (topPt.y > mid.y || (topPt.y == mid.y && topPt.x > mid.x)) {
395                topPt = mid;
396            }
397        }
398    }
399    return topPt;
400}
401
402// OPTIMIZE: avoid computing the unused half
403void xy_at_t(const Cubic& cubic, double t, double& x, double& y) {
404    _Point xy = xy_at_t(cubic, t);
405    if (&x) {
406        x = xy.x;
407    }
408    if (&y) {
409        y = xy.y;
410    }
411}
412
413_Point xy_at_t(const Cubic& cubic, double t) {
414    double one_t = 1 - t;
415    double one_t2 = one_t * one_t;
416    double a = one_t2 * one_t;
417    double b = 3 * one_t2 * t;
418    double t2 = t * t;
419    double c = 3 * one_t * t2;
420    double d = t2 * t;
421    _Point result = {a * cubic[0].x + b * cubic[1].x + c * cubic[2].x + d * cubic[3].x,
422            a * cubic[0].y + b * cubic[1].y + c * cubic[2].y + d * cubic[3].y};
423    return result;
424}
425