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 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 45 bool 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 /* 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 */ 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 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 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 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 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 217 double dx_at_t(const Quadratic& quad, double t) { 218 return derivativeAtT(&quad[0].x, t); 219 } 220 221 double 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 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 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
