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