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 "QuadraticUtilities.h"
10#include "TriangleUtilities.h"
11
12// from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
13double nearestT(const Quadratic& quad, const _Point& pt) {
14    _Vector pos = quad[0] - pt;
15    // search points P of bezier curve with PM.(dP / dt) = 0
16    // a calculus leads to a 3d degree equation :
17    _Vector A = quad[1] - quad[0];
18    _Vector B = quad[2] - quad[1];
19    B -= A;
20    double a = B.dot(B);
21    double b = 3 * A.dot(B);
22    double c = 2 * A.dot(A) + pos.dot(B);
23    double d = pos.dot(A);
24    double ts[3];
25    int roots = cubicRootsValidT(a, b, c, d, ts);
26    double d0 = pt.distanceSquared(quad[0]);
27    double d2 = pt.distanceSquared(quad[2]);
28    double distMin = SkTMin(d0, d2);
29    int bestIndex = -1;
30    for (int index = 0; index < roots; ++index) {
31        _Point onQuad;
32        xy_at_t(quad, ts[index], onQuad.x, onQuad.y);
33        double dist = pt.distanceSquared(onQuad);
34        if (distMin > dist) {
35            distMin = dist;
36            bestIndex = index;
37        }
38    }
39    if (bestIndex >= 0) {
40        return ts[bestIndex];
41    }
42    return d0 < d2 ? 0 : 1;
43}
44
45bool point_in_hull(const Quadratic& quad, const _Point& pt) {
46    return pointInTriangle((const Triangle&) quad, pt);
47}
48
49_Point top(const Quadratic& quad, double startT, double endT) {
50    Quadratic sub;
51    sub_divide(quad, startT, endT, sub);
52    _Point topPt = sub[0];
53    if (topPt.y > sub[2].y || (topPt.y == sub[2].y && topPt.x > sub[2].x)) {
54        topPt = sub[2];
55    }
56    if (!between(sub[0].y, sub[1].y, sub[2].y)) {
57        double extremeT;
58        if (findExtrema(sub[0].y, sub[1].y, sub[2].y, &extremeT)) {
59            extremeT = startT + (endT - startT) * extremeT;
60            _Point test;
61            xy_at_t(quad, extremeT, test.x, test.y);
62            if (topPt.y > test.y || (topPt.y == test.y && topPt.x > test.x)) {
63                topPt = test;
64            }
65        }
66    }
67    return topPt;
68}
69
70/*
71Numeric Solutions (5.6) suggests to solve the quadratic by computing
72       Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
73and using the roots
74      t1 = Q / A
75      t2 = C / Q
76*/
77int add_valid_ts(double s[], int realRoots, double* t) {
78    int foundRoots = 0;
79    for (int index = 0; index < realRoots; ++index) {
80        double tValue = s[index];
81        if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
82            if (approximately_less_than_zero(tValue)) {
83                tValue = 0;
84            } else if (approximately_greater_than_one(tValue)) {
85                tValue = 1;
86            }
87            for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
88                if (approximately_equal(t[idx2], tValue)) {
89                    goto nextRoot;
90                }
91            }
92            t[foundRoots++] = tValue;
93        }
94nextRoot:
95        ;
96    }
97    return foundRoots;
98}
99
100// note: caller expects multiple results to be sorted smaller first
101// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
102//  analysis of the quadratic equation, suggesting why the following looks at
103//  the sign of B -- and further suggesting that the greatest loss of precision
104//  is in b squared less two a c
105int quadraticRootsValidT(double A, double B, double C, double t[2]) {
106#if 0
107    B *= 2;
108    double square = B * B - 4 * A * C;
109    if (approximately_negative(square)) {
110        if (!approximately_positive(square)) {
111            return 0;
112        }
113        square = 0;
114    }
115    double squareRt = sqrt(square);
116    double Q = (B + (B < 0 ? -squareRt : squareRt)) / -2;
117    int foundRoots = 0;
118    double ratio = Q / A;
119    if (approximately_zero_or_more(ratio) && approximately_one_or_less(ratio)) {
120        if (approximately_less_than_zero(ratio)) {
121            ratio = 0;
122        } else if (approximately_greater_than_one(ratio)) {
123            ratio = 1;
124        }
125        t[0] = ratio;
126        ++foundRoots;
127    }
128    ratio = C / Q;
129    if (approximately_zero_or_more(ratio) && approximately_one_or_less(ratio)) {
130        if (approximately_less_than_zero(ratio)) {
131            ratio = 0;
132        } else if (approximately_greater_than_one(ratio)) {
133            ratio = 1;
134        }
135        if (foundRoots == 0 || !approximately_negative(ratio - t[0])) {
136            t[foundRoots++] = ratio;
137        } else if (!approximately_negative(t[0] - ratio)) {
138            t[foundRoots++] = t[0];
139            t[0] = ratio;
140        }
141    }
142#else
143    double s[2];
144    int realRoots = quadraticRootsReal(A, B, C, s);
145    int foundRoots = add_valid_ts(s, realRoots, t);
146#endif
147    return foundRoots;
148}
149
150// unlike quadratic roots, this does not discard real roots <= 0 or >= 1
151int quadraticRootsReal(const double A, const double B, const double C, double s[2]) {
152    const double p = B / (2 * A);
153    const double q = C / A;
154    if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) {
155        if (approximately_zero(B)) {
156            s[0] = 0;
157            return C == 0;
158        }
159        s[0] = -C / B;
160        return 1;
161    }
162    /* normal form: x^2 + px + q = 0 */
163    const double p2 = p * p;
164#if 0
165    double D = AlmostEqualUlps(p2, q) ? 0 : p2 - q;
166    if (D <= 0) {
167        if (D < 0) {
168            return 0;
169        }
170        s[0] = -p;
171        SkDebugf("[%d] %1.9g\n", 1, s[0]);
172        return 1;
173    }
174    double sqrt_D = sqrt(D);
175    s[0] = sqrt_D - p;
176    s[1] = -sqrt_D - p;
177    SkDebugf("[%d] %1.9g %1.9g\n", 2, s[0], s[1]);
178    return 2;
179#else
180    if (!AlmostEqualUlps(p2, q) && p2 < q) {
181        return 0;
182    }
183    double sqrt_D = 0;
184    if (p2 > q) {
185        sqrt_D = sqrt(p2 - q);
186    }
187    s[0] = sqrt_D - p;
188    s[1] = -sqrt_D - p;
189#if 0
190    if (AlmostEqualUlps(s[0], s[1])) {
191        SkDebugf("[%d] %1.9g\n", 1, s[0]);
192    } else {
193        SkDebugf("[%d] %1.9g %1.9g\n", 2, s[0], s[1]);
194    }
195#endif
196    return 1 + !AlmostEqualUlps(s[0], s[1]);
197#endif
198}
199
200void toCubic(const Quadratic& quad, Cubic& cubic) {
201    cubic[0] = quad[0];
202    cubic[2] = quad[1];
203    cubic[3] = quad[2];
204    cubic[1].x = (cubic[0].x + cubic[2].x * 2) / 3;
205    cubic[1].y = (cubic[0].y + cubic[2].y * 2) / 3;
206    cubic[2].x = (cubic[3].x + cubic[2].x * 2) / 3;
207    cubic[2].y = (cubic[3].y + cubic[2].y * 2) / 3;
208}
209
210static double derivativeAtT(const double* quad, double t) {
211    double a = t - 1;
212    double b = 1 - 2 * t;
213    double c = t;
214    return a * quad[0] + b * quad[2] + c * quad[4];
215}
216
217double dx_at_t(const Quadratic& quad, double t) {
218    return derivativeAtT(&quad[0].x, t);
219}
220
221double dy_at_t(const Quadratic& quad, double t) {
222    return derivativeAtT(&quad[0].y, t);
223}
224
225_Vector dxdy_at_t(const Quadratic& quad, double t) {
226    double a = t - 1;
227    double b = 1 - 2 * t;
228    double c = t;
229    _Vector result = { a * quad[0].x + b * quad[1].x + c * quad[2].x,
230            a * quad[0].y + b * quad[1].y + c * quad[2].y };
231    return result;
232}
233
234void xy_at_t(const Quadratic& quad, double t, double& x, double& y) {
235    double one_t = 1 - t;
236    double a = one_t * one_t;
237    double b = 2 * one_t * t;
238    double c = t * t;
239    if (&x) {
240        x = a * quad[0].x + b * quad[1].x + c * quad[2].x;
241    }
242    if (&y) {
243        y = a * quad[0].y + b * quad[1].y + c * quad[2].y;
244    }
245}
246
247_Point xy_at_t(const Quadratic& quad, double t) {
248    double one_t = 1 - t;
249    double a = one_t * one_t;
250    double b = 2 * one_t * t;
251    double c = t * t;
252    _Point result = { a * quad[0].x + b * quad[1].x + c * quad[2].x,
253            a * quad[0].y + b * quad[1].y + c * quad[2].y };
254    return result;
255}
256