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 "SkPathOpsCurve.h" 11 #include "SkPathOpsQuad.h" 12 #include "SkPathOpsRect.h" 13 14 // from blackpawn.com/texts/pointinpoly 15 static bool pointInTriangle(const SkDPoint fPts[3], const SkDPoint& test) { 16 SkDVector v0 = fPts[2] - fPts[0]; 17 SkDVector v1 = fPts[1] - fPts[0]; 18 SkDVector v2 = test - fPts[0]; 19 double dot00 = v0.dot(v0); 20 double dot01 = v0.dot(v1); 21 double dot02 = v0.dot(v2); 22 double dot11 = v1.dot(v1); 23 double dot12 = v1.dot(v2); 24 // Compute barycentric coordinates 25 double denom = dot00 * dot11 - dot01 * dot01; 26 double u = dot11 * dot02 - dot01 * dot12; 27 double v = dot00 * dot12 - dot01 * dot02; 28 // Check if point is in triangle 29 if (denom >= 0) { 30 return u >= 0 && v >= 0 && u + v < denom; 31 } 32 return u <= 0 && v <= 0 && u + v > denom; 33 } 34 35 static bool matchesEnd(const SkDPoint fPts[3], const SkDPoint& test) { 36 return fPts[0] == test || fPts[2] == test; 37 } 38 39 /* started with at_most_end_pts_in_common from SkDQuadIntersection.cpp */ 40 // Do a quick reject by rotating all points relative to a line formed by 41 // a pair of one quad's points. If the 2nd quad's points 42 // are on the line or on the opposite side from the 1st quad's 'odd man', the 43 // curves at most intersect at the endpoints. 44 /* if returning true, check contains true if quad's hull collapsed, making the cubic linear 45 if returning false, check contains true if the the quad pair have only the end point in common 46 */ 47 bool SkDQuad::hullIntersects(const SkDQuad& q2, bool* isLinear) const { 48 bool linear = true; 49 for (int oddMan = 0; oddMan < kPointCount; ++oddMan) { 50 const SkDPoint* endPt[2]; 51 this->otherPts(oddMan, endPt); 52 double origX = endPt[0]->fX; 53 double origY = endPt[0]->fY; 54 double adj = endPt[1]->fX - origX; 55 double opp = endPt[1]->fY - origY; 56 double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp; 57 if (approximately_zero(sign)) { 58 continue; 59 } 60 linear = false; 61 bool foundOutlier = false; 62 for (int n = 0; n < kPointCount; ++n) { 63 double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp; 64 if (test * sign > 0 && !precisely_zero(test)) { 65 foundOutlier = true; 66 break; 67 } 68 } 69 if (!foundOutlier) { 70 return false; 71 } 72 } 73 if (linear && !matchesEnd(fPts, q2.fPts[0]) && !matchesEnd(fPts, q2.fPts[2])) { 74 // if the end point of the opposite quad is inside the hull that is nearly a line, 75 // then representing the quad as a line may cause the intersection to be missed. 76 // Check to see if the endpoint is in the triangle. 77 if (pointInTriangle(fPts, q2.fPts[0]) || pointInTriangle(fPts, q2.fPts[2])) { 78 linear = false; 79 } 80 } 81 *isLinear = linear; 82 return true; 83 } 84 85 bool SkDQuad::hullIntersects(const SkDConic& conic, bool* isLinear) const { 86 return conic.hullIntersects(*this, isLinear); 87 } 88 89 bool SkDQuad::hullIntersects(const SkDCubic& cubic, bool* isLinear) const { 90 return cubic.hullIntersects(*this, isLinear); 91 } 92 93 /* bit twiddling for finding the off curve index (x&~m is the pair in [0,1,2] excluding oddMan) 94 oddMan opp x=oddMan^opp x=x-oddMan m=x>>2 x&~m 95 0 1 1 1 0 1 96 2 2 2 0 2 97 1 1 0 -1 -1 0 98 2 3 2 0 2 99 2 1 3 1 0 1 100 2 0 -2 -1 0 101 */ 102 void SkDQuad::otherPts(int oddMan, const SkDPoint* endPt[2]) const { 103 for (int opp = 1; opp < kPointCount; ++opp) { 104 int end = (oddMan ^ opp) - oddMan; // choose a value not equal to oddMan 105 end &= ~(end >> 2); // if the value went negative, set it to zero 106 endPt[opp - 1] = &fPts[end]; 107 } 108 } 109 110 int SkDQuad::AddValidTs(double s[], int realRoots, double* t) { 111 int foundRoots = 0; 112 for (int index = 0; index < realRoots; ++index) { 113 double tValue = s[index]; 114 if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) { 115 if (approximately_less_than_zero(tValue)) { 116 tValue = 0; 117 } else if (approximately_greater_than_one(tValue)) { 118 tValue = 1; 119 } 120 for (int idx2 = 0; idx2 < foundRoots; ++idx2) { 121 if (approximately_equal(t[idx2], tValue)) { 122 goto nextRoot; 123 } 124 } 125 t[foundRoots++] = tValue; 126 } 127 nextRoot: 128 {} 129 } 130 return foundRoots; 131 } 132 133 // note: caller expects multiple results to be sorted smaller first 134 // note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting 135 // analysis of the quadratic equation, suggesting why the following looks at 136 // the sign of B -- and further suggesting that the greatest loss of precision 137 // is in b squared less two a c 138 int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) { 139 double s[2]; 140 int realRoots = RootsReal(A, B, C, s); 141 int foundRoots = AddValidTs(s, realRoots, t); 142 return foundRoots; 143 } 144 145 static int handle_zero(const double B, const double C, double s[2]) { 146 if (approximately_zero(B)) { 147 s[0] = 0; 148 return C == 0; 149 } 150 s[0] = -C / B; 151 return 1; 152 } 153 154 /* 155 Numeric Solutions (5.6) suggests to solve the quadratic by computing 156 Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C)) 157 and using the roots 158 t1 = Q / A 159 t2 = C / Q 160 */ 161 // this does not discard real roots <= 0 or >= 1 162 int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) { 163 if (!A) { 164 return handle_zero(B, C, s); 165 } 166 const double p = B / (2 * A); 167 const double q = C / A; 168 if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) { 169 return handle_zero(B, C, s); 170 } 171 /* normal form: x^2 + px + q = 0 */ 172 const double p2 = p * p; 173 if (!AlmostDequalUlps(p2, q) && p2 < q) { 174 return 0; 175 } 176 double sqrt_D = 0; 177 if (p2 > q) { 178 sqrt_D = sqrt(p2 - q); 179 } 180 s[0] = sqrt_D - p; 181 s[1] = -sqrt_D - p; 182 return 1 + !AlmostDequalUlps(s[0], s[1]); 183 } 184 185 bool SkDQuad::isLinear(int startIndex, int endIndex) const { 186 SkLineParameters lineParameters; 187 lineParameters.quadEndPoints(*this, startIndex, endIndex); 188 // FIXME: maybe it's possible to avoid this and compare non-normalized 189 lineParameters.normalize(); 190 double distance = lineParameters.controlPtDistance(*this); 191 double tiniest = SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(fPts[0].fX, fPts[0].fY), 192 fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY); 193 double largest = SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(fPts[0].fX, fPts[0].fY), 194 fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY); 195 largest = SkTMax(largest, -tiniest); 196 return approximately_zero_when_compared_to(distance, largest); 197 } 198 199 SkDVector SkDQuad::dxdyAtT(double t) const { 200 double a = t - 1; 201 double b = 1 - 2 * t; 202 double c = t; 203 SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX, 204 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY }; 205 if (result.fX == 0 && result.fY == 0) { 206 if (zero_or_one(t)) { 207 result = fPts[2] - fPts[0]; 208 } else { 209 // incomplete 210 SkDebugf("!q"); 211 } 212 } 213 return result; 214 } 215 216 // OPTIMIZE: assert if caller passes in t == 0 / t == 1 ? 217 SkDPoint SkDQuad::ptAtT(double t) const { 218 if (0 == t) { 219 return fPts[0]; 220 } 221 if (1 == t) { 222 return fPts[2]; 223 } 224 double one_t = 1 - t; 225 double a = one_t * one_t; 226 double b = 2 * one_t * t; 227 double c = t * t; 228 SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX, 229 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY }; 230 return result; 231 } 232 233 static double interp_quad_coords(const double* src, double t) { 234 if (0 == t) { 235 return src[0]; 236 } 237 if (1 == t) { 238 return src[4]; 239 } 240 double ab = SkDInterp(src[0], src[2], t); 241 double bc = SkDInterp(src[2], src[4], t); 242 double abc = SkDInterp(ab, bc, t); 243 return abc; 244 } 245 246 bool SkDQuad::monotonicInX() const { 247 return between(fPts[0].fX, fPts[1].fX, fPts[2].fX); 248 } 249 250 bool SkDQuad::monotonicInY() const { 251 return between(fPts[0].fY, fPts[1].fY, fPts[2].fY); 252 } 253 254 /* 255 Given a quadratic q, t1, and t2, find a small quadratic segment. 256 257 The new quadratic is defined by A, B, and C, where 258 A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1 259 C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1 260 261 To find B, compute the point halfway between t1 and t2: 262 263 q(at (t1 + t2)/2) == D 264 265 Next, compute where D must be if we know the value of B: 266 267 _12 = A/2 + B/2 268 12_ = B/2 + C/2 269 123 = A/4 + B/2 + C/4 270 = D 271 272 Group the known values on one side: 273 274 B = D*2 - A/2 - C/2 275 */ 276 277 // OPTIMIZE? : special case t1 = 1 && t2 = 0 278 SkDQuad SkDQuad::subDivide(double t1, double t2) const { 279 if (0 == t1 && 1 == t2) { 280 return *this; 281 } 282 SkDQuad dst; 283 double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1); 284 double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1); 285 double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2); 286 double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2); 287 double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2); 288 double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2); 289 /* bx = */ dst[1].fX = 2 * dx - (ax + cx) / 2; 290 /* by = */ dst[1].fY = 2 * dy - (ay + cy) / 2; 291 return dst; 292 } 293 294 void SkDQuad::align(int endIndex, SkDPoint* dstPt) const { 295 if (fPts[endIndex].fX == fPts[1].fX) { 296 dstPt->fX = fPts[endIndex].fX; 297 } 298 if (fPts[endIndex].fY == fPts[1].fY) { 299 dstPt->fY = fPts[endIndex].fY; 300 } 301 } 302 303 SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const { 304 SkASSERT(t1 != t2); 305 SkDPoint b; 306 SkDQuad sub = subDivide(t1, t2); 307 SkDLine b0 = {{a, sub[1] + (a - sub[0])}}; 308 SkDLine b1 = {{c, sub[1] + (c - sub[2])}}; 309 SkIntersections i; 310 i.intersectRay(b0, b1); 311 if (i.used() == 1 && i[0][0] >= 0 && i[1][0] >= 0) { 312 b = i.pt(0); 313 } else { 314 SkASSERT(i.used() <= 2); 315 return SkDPoint::Mid(b0[1], b1[1]); 316 } 317 if (t1 == 0 || t2 == 0) { 318 align(0, &b); 319 } 320 if (t1 == 1 || t2 == 1) { 321 align(2, &b); 322 } 323 if (AlmostBequalUlps(b.fX, a.fX)) { 324 b.fX = a.fX; 325 } else if (AlmostBequalUlps(b.fX, c.fX)) { 326 b.fX = c.fX; 327 } 328 if (AlmostBequalUlps(b.fY, a.fY)) { 329 b.fY = a.fY; 330 } else if (AlmostBequalUlps(b.fY, c.fY)) { 331 b.fY = c.fY; 332 } 333 return b; 334 } 335 336 /* classic one t subdivision */ 337 static void interp_quad_coords(const double* src, double* dst, double t) { 338 double ab = SkDInterp(src[0], src[2], t); 339 double bc = SkDInterp(src[2], src[4], t); 340 dst[0] = src[0]; 341 dst[2] = ab; 342 dst[4] = SkDInterp(ab, bc, t); 343 dst[6] = bc; 344 dst[8] = src[4]; 345 } 346 347 SkDQuadPair SkDQuad::chopAt(double t) const 348 { 349 SkDQuadPair dst; 350 interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t); 351 interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t); 352 return dst; 353 } 354 355 static int valid_unit_divide(double numer, double denom, double* ratio) 356 { 357 if (numer < 0) { 358 numer = -numer; 359 denom = -denom; 360 } 361 if (denom == 0 || numer == 0 || numer >= denom) { 362 return 0; 363 } 364 double r = numer / denom; 365 if (r == 0) { // catch underflow if numer <<<< denom 366 return 0; 367 } 368 *ratio = r; 369 return 1; 370 } 371 372 /** Quad'(t) = At + B, where 373 A = 2(a - 2b + c) 374 B = 2(b - a) 375 Solve for t, only if it fits between 0 < t < 1 376 */ 377 int SkDQuad::FindExtrema(const double src[], double tValue[1]) { 378 /* At + B == 0 379 t = -B / A 380 */ 381 double a = src[0]; 382 double b = src[2]; 383 double c = src[4]; 384 return valid_unit_divide(a - b, a - b - b + c, tValue); 385 } 386 387 /* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t) 388 * 389 * a = A - 2*B + C 390 * b = 2*B - 2*C 391 * c = C 392 */ 393 void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) { 394 *a = quad[0]; // a = A 395 *b = 2 * quad[2]; // b = 2*B 396 *c = quad[4]; // c = C 397 *b -= *c; // b = 2*B - C 398 *a -= *b; // a = A - 2*B + C 399 *b -= *c; // b = 2*B - 2*C 400 } 401 402 int SkTQuad::intersectRay(SkIntersections* i, const SkDLine& line) const { 403 return i->intersectRay(fQuad, line); 404 } 405 406 bool SkTQuad::hullIntersects(const SkDConic& conic, bool* isLinear) const { 407 return conic.hullIntersects(fQuad, isLinear); 408 } 409 410 bool SkTQuad::hullIntersects(const SkDCubic& cubic, bool* isLinear) const { 411 return cubic.hullIntersects(fQuad, isLinear); 412 } 413 414 void SkTQuad::setBounds(SkDRect* rect) const { 415 rect->setBounds(fQuad); 416 } 417
