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 "SkIntersections.h" 8 #include "SkLineParameters.h" 9 #include "SkPathOpsCubic.h" 10 #include "SkPathOpsQuad.h" 11 #include "SkPathOpsTriangle.h" 12 13 // from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html 14 // (currently only used by testing) 15 double SkDQuad::nearestT(const SkDPoint& pt) const { 16 SkDVector pos = fPts[0] - pt; 17 // search points P of bezier curve with PM.(dP / dt) = 0 18 // a calculus leads to a 3d degree equation : 19 SkDVector A = fPts[1] - fPts[0]; 20 SkDVector B = fPts[2] - fPts[1]; 21 B -= A; 22 double a = B.dot(B); 23 double b = 3 * A.dot(B); 24 double c = 2 * A.dot(A) + pos.dot(B); 25 double d = pos.dot(A); 26 double ts[3]; 27 int roots = SkDCubic::RootsValidT(a, b, c, d, ts); 28 double d0 = pt.distanceSquared(fPts[0]); 29 double d2 = pt.distanceSquared(fPts[2]); 30 double distMin = SkTMin(d0, d2); 31 int bestIndex = -1; 32 for (int index = 0; index < roots; ++index) { 33 SkDPoint onQuad = ptAtT(ts[index]); 34 double dist = pt.distanceSquared(onQuad); 35 if (distMin > dist) { 36 distMin = dist; 37 bestIndex = index; 38 } 39 } 40 if (bestIndex >= 0) { 41 return ts[bestIndex]; 42 } 43 return d0 < d2 ? 0 : 1; 44 } 45 46 bool SkDQuad::pointInHull(const SkDPoint& pt) const { 47 return ((const SkDTriangle&) fPts).contains(pt); 48 } 49 50 SkDPoint SkDQuad::top(double startT, double endT) const { 51 SkDQuad sub = subDivide(startT, endT); 52 SkDPoint topPt = sub[0]; 53 if (topPt.fY > sub[2].fY || (topPt.fY == sub[2].fY && topPt.fX > sub[2].fX)) { 54 topPt = sub[2]; 55 } 56 if (!between(sub[0].fY, sub[1].fY, sub[2].fY)) { 57 double extremeT; 58 if (FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, &extremeT)) { 59 extremeT = startT + (endT - startT) * extremeT; 60 SkDPoint test = ptAtT(extremeT); 61 if (topPt.fY > test.fY || (topPt.fY == test.fY && topPt.fX > test.fX)) { 62 topPt = test; 63 } 64 } 65 } 66 return topPt; 67 } 68 69 int SkDQuad::AddValidTs(double s[], int realRoots, double* t) { 70 int foundRoots = 0; 71 for (int index = 0; index < realRoots; ++index) { 72 double tValue = s[index]; 73 if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) { 74 if (approximately_less_than_zero(tValue)) { 75 tValue = 0; 76 } else if (approximately_greater_than_one(tValue)) { 77 tValue = 1; 78 } 79 for (int idx2 = 0; idx2 < foundRoots; ++idx2) { 80 if (approximately_equal(t[idx2], tValue)) { 81 goto nextRoot; 82 } 83 } 84 t[foundRoots++] = tValue; 85 } 86 nextRoot: 87 {} 88 } 89 return foundRoots; 90 } 91 92 // note: caller expects multiple results to be sorted smaller first 93 // note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting 94 // analysis of the quadratic equation, suggesting why the following looks at 95 // the sign of B -- and further suggesting that the greatest loss of precision 96 // is in b squared less two a c 97 int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) { 98 double s[2]; 99 int realRoots = RootsReal(A, B, C, s); 100 int foundRoots = AddValidTs(s, realRoots, t); 101 return foundRoots; 102 } 103 104 /* 105 Numeric Solutions (5.6) suggests to solve the quadratic by computing 106 Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C)) 107 and using the roots 108 t1 = Q / A 109 t2 = C / Q 110 */ 111 // this does not discard real roots <= 0 or >= 1 112 int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) { 113 const double p = B / (2 * A); 114 const double q = C / A; 115 if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) { 116 if (approximately_zero(B)) { 117 s[0] = 0; 118 return C == 0; 119 } 120 s[0] = -C / B; 121 return 1; 122 } 123 /* normal form: x^2 + px + q = 0 */ 124 const double p2 = p * p; 125 if (!AlmostEqualUlps(p2, q) && p2 < q) { 126 return 0; 127 } 128 double sqrt_D = 0; 129 if (p2 > q) { 130 sqrt_D = sqrt(p2 - q); 131 } 132 s[0] = sqrt_D - p; 133 s[1] = -sqrt_D - p; 134 return 1 + !AlmostEqualUlps(s[0], s[1]); 135 } 136 137 bool SkDQuad::isLinear(int startIndex, int endIndex) const { 138 SkLineParameters lineParameters; 139 lineParameters.quadEndPoints(*this, startIndex, endIndex); 140 // FIXME: maybe it's possible to avoid this and compare non-normalized 141 lineParameters.normalize(); 142 double distance = lineParameters.controlPtDistance(*this); 143 return approximately_zero(distance); 144 } 145 146 SkDCubic SkDQuad::toCubic() const { 147 SkDCubic cubic; 148 cubic[0] = fPts[0]; 149 cubic[2] = fPts[1]; 150 cubic[3] = fPts[2]; 151 cubic[1].fX = (cubic[0].fX + cubic[2].fX * 2) / 3; 152 cubic[1].fY = (cubic[0].fY + cubic[2].fY * 2) / 3; 153 cubic[2].fX = (cubic[3].fX + cubic[2].fX * 2) / 3; 154 cubic[2].fY = (cubic[3].fY + cubic[2].fY * 2) / 3; 155 return cubic; 156 } 157 158 SkDVector SkDQuad::dxdyAtT(double t) const { 159 double a = t - 1; 160 double b = 1 - 2 * t; 161 double c = t; 162 SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX, 163 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY }; 164 return result; 165 } 166 167 // OPTIMIZE: assert if caller passes in t == 0 / t == 1 ? 168 SkDPoint SkDQuad::ptAtT(double t) const { 169 if (0 == t) { 170 return fPts[0]; 171 } 172 if (1 == t) { 173 return fPts[2]; 174 } 175 double one_t = 1 - t; 176 double a = one_t * one_t; 177 double b = 2 * one_t * t; 178 double c = t * t; 179 SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX, 180 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY }; 181 return result; 182 } 183 184 /* 185 Given a quadratic q, t1, and t2, find a small quadratic segment. 186 187 The new quadratic is defined by A, B, and C, where 188 A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1 189 C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1 190 191 To find B, compute the point halfway between t1 and t2: 192 193 q(at (t1 + t2)/2) == D 194 195 Next, compute where D must be if we know the value of B: 196 197 _12 = A/2 + B/2 198 12_ = B/2 + C/2 199 123 = A/4 + B/2 + C/4 200 = D 201 202 Group the known values on one side: 203 204 B = D*2 - A/2 - C/2 205 */ 206 207 static double interp_quad_coords(const double* src, double t) { 208 double ab = SkDInterp(src[0], src[2], t); 209 double bc = SkDInterp(src[2], src[4], t); 210 double abc = SkDInterp(ab, bc, t); 211 return abc; 212 } 213 214 bool SkDQuad::monotonicInY() const { 215 return between(fPts[0].fY, fPts[1].fY, fPts[2].fY); 216 } 217 218 SkDQuad SkDQuad::subDivide(double t1, double t2) const { 219 SkDQuad dst; 220 double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1); 221 double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1); 222 double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2); 223 double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2); 224 double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2); 225 double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2); 226 /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2; 227 /* by = */ dst[1].fY = 2*dy - (ay + cy)/2; 228 return dst; 229 } 230 231 void SkDQuad::align(int endIndex, SkDPoint* dstPt) const { 232 if (fPts[endIndex].fX == fPts[1].fX) { 233 dstPt->fX = fPts[endIndex].fX; 234 } 235 if (fPts[endIndex].fY == fPts[1].fY) { 236 dstPt->fY = fPts[endIndex].fY; 237 } 238 } 239 240 SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const { 241 SkASSERT(t1 != t2); 242 SkDPoint b; 243 #if 0 244 // this approach assumes that the control point computed directly is accurate enough 245 double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2); 246 double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2); 247 b.fX = 2 * dx - (a.fX + c.fX) / 2; 248 b.fY = 2 * dy - (a.fY + c.fY) / 2; 249 #else 250 SkDQuad sub = subDivide(t1, t2); 251 SkDLine b0 = {{a, sub[1] + (a - sub[0])}}; 252 SkDLine b1 = {{c, sub[1] + (c - sub[2])}}; 253 SkIntersections i; 254 i.intersectRay(b0, b1); 255 if (i.used() == 1) { 256 b = i.pt(0); 257 } else { 258 SkASSERT(i.used() == 2 || i.used() == 0); 259 b = SkDPoint::Mid(b0[1], b1[1]); 260 } 261 #endif 262 if (t1 == 0 || t2 == 0) { 263 align(0, &b); 264 } 265 if (t1 == 1 || t2 == 1) { 266 align(2, &b); 267 } 268 if (precisely_subdivide_equal(b.fX, a.fX)) { 269 b.fX = a.fX; 270 } else if (precisely_subdivide_equal(b.fX, c.fX)) { 271 b.fX = c.fX; 272 } 273 if (precisely_subdivide_equal(b.fY, a.fY)) { 274 b.fY = a.fY; 275 } else if (precisely_subdivide_equal(b.fY, c.fY)) { 276 b.fY = c.fY; 277 } 278 return b; 279 } 280 281 /* classic one t subdivision */ 282 static void interp_quad_coords(const double* src, double* dst, double t) { 283 double ab = SkDInterp(src[0], src[2], t); 284 double bc = SkDInterp(src[2], src[4], t); 285 dst[0] = src[0]; 286 dst[2] = ab; 287 dst[4] = SkDInterp(ab, bc, t); 288 dst[6] = bc; 289 dst[8] = src[4]; 290 } 291 292 SkDQuadPair SkDQuad::chopAt(double t) const 293 { 294 SkDQuadPair dst; 295 interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t); 296 interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t); 297 return dst; 298 } 299 300 static int valid_unit_divide(double numer, double denom, double* ratio) 301 { 302 if (numer < 0) { 303 numer = -numer; 304 denom = -denom; 305 } 306 if (denom == 0 || numer == 0 || numer >= denom) { 307 return 0; 308 } 309 double r = numer / denom; 310 if (r == 0) { // catch underflow if numer <<<< denom 311 return 0; 312 } 313 *ratio = r; 314 return 1; 315 } 316 317 /** Quad'(t) = At + B, where 318 A = 2(a - 2b + c) 319 B = 2(b - a) 320 Solve for t, only if it fits between 0 < t < 1 321 */ 322 int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) { 323 /* At + B == 0 324 t = -B / A 325 */ 326 return valid_unit_divide(a - b, a - b - b + c, tValue); 327 } 328 329 /* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t) 330 * 331 * a = A - 2*B + C 332 * b = 2*B - 2*C 333 * c = C 334 */ 335 void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) { 336 *a = quad[0]; // a = A 337 *b = 2 * quad[2]; // b = 2*B 338 *c = quad[4]; // c = C 339 *b -= *c; // b = 2*B - C 340 *a -= *b; // a = A - 2*B + C 341 *b -= *c; // b = 2*B - 2*C 342 } 343