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
nearestT(const Quadratic & quad,const _Point & pt)13 double 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
point_in_hull(const Quadratic & quad,const _Point & pt)45 bool point_in_hull(const Quadratic& quad, const _Point& pt) {
46 return pointInTriangle((const Triangle&) quad, pt);
47 }
48
top(const Quadratic & quad,double startT,double endT)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 /*
71 Numeric Solutions (5.6) suggests to solve the quadratic by computing
72 Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
73 and using the roots
74 t1 = Q / A
75 t2 = C / Q
76 */
add_valid_ts(double s[],int realRoots,double * t)77 int 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 }
94 nextRoot:
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
quadraticRootsValidT(double A,double B,double C,double t[2])105 int 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
quadraticRootsReal(const double A,const double B,const double C,double s[2])151 int 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
toCubic(const Quadratic & quad,Cubic & cubic)200 void 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
derivativeAtT(const double * quad,double t)210 static 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
dx_at_t(const Quadratic & quad,double t)217 double dx_at_t(const Quadratic& quad, double t) {
218 return derivativeAtT(&quad[0].x, t);
219 }
220
dy_at_t(const Quadratic & quad,double t)221 double dy_at_t(const Quadratic& quad, double t) {
222 return derivativeAtT(&quad[0].y, t);
223 }
224
dxdy_at_t(const Quadratic & quad,double t)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
xy_at_t(const Quadratic & quad,double t,double & x,double & y)234 void 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
xy_at_t(const Quadratic & quad,double t)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