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