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 "SkGeometry.h" 8 #include "SkLineParameters.h" 9 #include "SkPathOpsConic.h" 10 #include "SkPathOpsCubic.h" 11 #include "SkPathOpsCurve.h" 12 #include "SkPathOpsLine.h" 13 #include "SkPathOpsQuad.h" 14 #include "SkPathOpsRect.h" 15 #include "SkTSort.h" 16 17 const int SkDCubic::gPrecisionUnit = 256; // FIXME: test different values in test framework 18 19 void SkDCubic::align(int endIndex, int ctrlIndex, SkDPoint* dstPt) const { 20 if (fPts[endIndex].fX == fPts[ctrlIndex].fX) { 21 dstPt->fX = fPts[endIndex].fX; 22 } 23 if (fPts[endIndex].fY == fPts[ctrlIndex].fY) { 24 dstPt->fY = fPts[endIndex].fY; 25 } 26 } 27 28 // give up when changing t no longer moves point 29 // also, copy point rather than recompute it when it does change 30 double SkDCubic::binarySearch(double min, double max, double axisIntercept, 31 SearchAxis xAxis) const { 32 double t = (min + max) / 2; 33 double step = (t - min) / 2; 34 SkDPoint cubicAtT = ptAtT(t); 35 double calcPos = (&cubicAtT.fX)[xAxis]; 36 double calcDist = calcPos - axisIntercept; 37 do { 38 double priorT = t - step; 39 SkOPASSERT(priorT >= min); 40 SkDPoint lessPt = ptAtT(priorT); 41 if (approximately_equal_half(lessPt.fX, cubicAtT.fX) 42 && approximately_equal_half(lessPt.fY, cubicAtT.fY)) { 43 return -1; // binary search found no point at this axis intercept 44 } 45 double lessDist = (&lessPt.fX)[xAxis] - axisIntercept; 46 #if DEBUG_CUBIC_BINARY_SEARCH 47 SkDebugf("t=%1.9g calc=%1.9g dist=%1.9g step=%1.9g less=%1.9g\n", t, calcPos, calcDist, 48 step, lessDist); 49 #endif 50 double lastStep = step; 51 step /= 2; 52 if (calcDist > 0 ? calcDist > lessDist : calcDist < lessDist) { 53 t = priorT; 54 } else { 55 double nextT = t + lastStep; 56 if (nextT > max) { 57 return -1; 58 } 59 SkDPoint morePt = ptAtT(nextT); 60 if (approximately_equal_half(morePt.fX, cubicAtT.fX) 61 && approximately_equal_half(morePt.fY, cubicAtT.fY)) { 62 return -1; // binary search found no point at this axis intercept 63 } 64 double moreDist = (&morePt.fX)[xAxis] - axisIntercept; 65 if (calcDist > 0 ? calcDist <= moreDist : calcDist >= moreDist) { 66 continue; 67 } 68 t = nextT; 69 } 70 SkDPoint testAtT = ptAtT(t); 71 cubicAtT = testAtT; 72 calcPos = (&cubicAtT.fX)[xAxis]; 73 calcDist = calcPos - axisIntercept; 74 } while (!approximately_equal(calcPos, axisIntercept)); 75 return t; 76 } 77 78 // get the rough scale of the cubic; used to determine if curvature is extreme 79 double SkDCubic::calcPrecision() const { 80 return ((fPts[1] - fPts[0]).length() 81 + (fPts[2] - fPts[1]).length() 82 + (fPts[3] - fPts[2]).length()) / gPrecisionUnit; 83 } 84 85 /* classic one t subdivision */ 86 static void interp_cubic_coords(const double* src, double* dst, double t) { 87 double ab = SkDInterp(src[0], src[2], t); 88 double bc = SkDInterp(src[2], src[4], t); 89 double cd = SkDInterp(src[4], src[6], t); 90 double abc = SkDInterp(ab, bc, t); 91 double bcd = SkDInterp(bc, cd, t); 92 double abcd = SkDInterp(abc, bcd, t); 93 94 dst[0] = src[0]; 95 dst[2] = ab; 96 dst[4] = abc; 97 dst[6] = abcd; 98 dst[8] = bcd; 99 dst[10] = cd; 100 dst[12] = src[6]; 101 } 102 103 SkDCubicPair SkDCubic::chopAt(double t) const { 104 SkDCubicPair dst; 105 if (t == 0.5) { 106 dst.pts[0] = fPts[0]; 107 dst.pts[1].fX = (fPts[0].fX + fPts[1].fX) / 2; 108 dst.pts[1].fY = (fPts[0].fY + fPts[1].fY) / 2; 109 dst.pts[2].fX = (fPts[0].fX + 2 * fPts[1].fX + fPts[2].fX) / 4; 110 dst.pts[2].fY = (fPts[0].fY + 2 * fPts[1].fY + fPts[2].fY) / 4; 111 dst.pts[3].fX = (fPts[0].fX + 3 * (fPts[1].fX + fPts[2].fX) + fPts[3].fX) / 8; 112 dst.pts[3].fY = (fPts[0].fY + 3 * (fPts[1].fY + fPts[2].fY) + fPts[3].fY) / 8; 113 dst.pts[4].fX = (fPts[1].fX + 2 * fPts[2].fX + fPts[3].fX) / 4; 114 dst.pts[4].fY = (fPts[1].fY + 2 * fPts[2].fY + fPts[3].fY) / 4; 115 dst.pts[5].fX = (fPts[2].fX + fPts[3].fX) / 2; 116 dst.pts[5].fY = (fPts[2].fY + fPts[3].fY) / 2; 117 dst.pts[6] = fPts[3]; 118 return dst; 119 } 120 interp_cubic_coords(&fPts[0].fX, &dst.pts[0].fX, t); 121 interp_cubic_coords(&fPts[0].fY, &dst.pts[0].fY, t); 122 return dst; 123 } 124 125 void SkDCubic::Coefficients(const double* src, double* A, double* B, double* C, double* D) { 126 *A = src[6]; // d 127 *B = src[4] * 3; // 3*c 128 *C = src[2] * 3; // 3*b 129 *D = src[0]; // a 130 *A -= *D - *C + *B; // A = -a + 3*b - 3*c + d 131 *B += 3 * *D - 2 * *C; // B = 3*a - 6*b + 3*c 132 *C -= 3 * *D; // C = -3*a + 3*b 133 } 134 135 bool SkDCubic::endsAreExtremaInXOrY() const { 136 return (between(fPts[0].fX, fPts[1].fX, fPts[3].fX) 137 && between(fPts[0].fX, fPts[2].fX, fPts[3].fX)) 138 || (between(fPts[0].fY, fPts[1].fY, fPts[3].fY) 139 && between(fPts[0].fY, fPts[2].fY, fPts[3].fY)); 140 } 141 142 // Do a quick reject by rotating all points relative to a line formed by 143 // a pair of one cubic's points. If the 2nd cubic's points 144 // are on the line or on the opposite side from the 1st cubic's 'odd man', the 145 // curves at most intersect at the endpoints. 146 /* if returning true, check contains true if cubic's hull collapsed, making the cubic linear 147 if returning false, check contains true if the the cubic pair have only the end point in common 148 */ 149 bool SkDCubic::hullIntersects(const SkDPoint* pts, int ptCount, bool* isLinear) const { 150 bool linear = true; 151 char hullOrder[4]; 152 int hullCount = convexHull(hullOrder); 153 int end1 = hullOrder[0]; 154 int hullIndex = 0; 155 const SkDPoint* endPt[2]; 156 endPt[0] = &fPts[end1]; 157 do { 158 hullIndex = (hullIndex + 1) % hullCount; 159 int end2 = hullOrder[hullIndex]; 160 endPt[1] = &fPts[end2]; 161 double origX = endPt[0]->fX; 162 double origY = endPt[0]->fY; 163 double adj = endPt[1]->fX - origX; 164 double opp = endPt[1]->fY - origY; 165 int oddManMask = other_two(end1, end2); 166 int oddMan = end1 ^ oddManMask; 167 double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp; 168 int oddMan2 = end2 ^ oddManMask; 169 double sign2 = (fPts[oddMan2].fY - origY) * adj - (fPts[oddMan2].fX - origX) * opp; 170 if (sign * sign2 < 0) { 171 continue; 172 } 173 if (approximately_zero(sign)) { 174 sign = sign2; 175 if (approximately_zero(sign)) { 176 continue; 177 } 178 } 179 linear = false; 180 bool foundOutlier = false; 181 for (int n = 0; n < ptCount; ++n) { 182 double test = (pts[n].fY - origY) * adj - (pts[n].fX - origX) * opp; 183 if (test * sign > 0 && !precisely_zero(test)) { 184 foundOutlier = true; 185 break; 186 } 187 } 188 if (!foundOutlier) { 189 return false; 190 } 191 endPt[0] = endPt[1]; 192 end1 = end2; 193 } while (hullIndex); 194 *isLinear = linear; 195 return true; 196 } 197 198 bool SkDCubic::hullIntersects(const SkDCubic& c2, bool* isLinear) const { 199 return hullIntersects(c2.fPts, c2.kPointCount, isLinear); 200 } 201 202 bool SkDCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const { 203 return hullIntersects(quad.fPts, quad.kPointCount, isLinear); 204 } 205 206 bool SkDCubic::hullIntersects(const SkDConic& conic, bool* isLinear) const { 207 208 return hullIntersects(conic.fPts, isLinear); 209 } 210 211 bool SkDCubic::isLinear(int startIndex, int endIndex) const { 212 if (fPts[0].approximatelyDEqual(fPts[3])) { 213 return ((const SkDQuad *) this)->isLinear(0, 2); 214 } 215 SkLineParameters lineParameters; 216 lineParameters.cubicEndPoints(*this, startIndex, endIndex); 217 // FIXME: maybe it's possible to avoid this and compare non-normalized 218 lineParameters.normalize(); 219 double tiniest = SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(fPts[0].fX, fPts[0].fY), 220 fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY), fPts[3].fX), fPts[3].fY); 221 double largest = SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(fPts[0].fX, fPts[0].fY), 222 fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY), fPts[3].fX), fPts[3].fY); 223 largest = SkTMax(largest, -tiniest); 224 double distance = lineParameters.controlPtDistance(*this, 1); 225 if (!approximately_zero_when_compared_to(distance, largest)) { 226 return false; 227 } 228 distance = lineParameters.controlPtDistance(*this, 2); 229 return approximately_zero_when_compared_to(distance, largest); 230 } 231 232 // from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf 233 // c(t) = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3 234 // c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2 235 // = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2 236 static double derivative_at_t(const double* src, double t) { 237 double one_t = 1 - t; 238 double a = src[0]; 239 double b = src[2]; 240 double c = src[4]; 241 double d = src[6]; 242 return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t); 243 } 244 245 int SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) { 246 SkDCubic cubic; 247 cubic.set(pointsPtr); 248 if (cubic.monotonicInX() && cubic.monotonicInY()) { 249 return 0; 250 } 251 double tt[2], ss[2]; 252 SkCubicType cubicType = SkClassifyCubic(pointsPtr, tt, ss); 253 switch (cubicType) { 254 case SkCubicType::kLoop: { 255 const double &td = tt[0], &te = tt[1], &sd = ss[0], &se = ss[1]; 256 if (roughly_between(0, td, sd) && roughly_between(0, te, se)) { 257 SkASSERT(roughly_between(0, td/sd, 1) && roughly_between(0, te/se, 1)); 258 t[0] = static_cast<SkScalar>((td * se + te * sd) / (2 * sd * se)); 259 SkASSERT(roughly_between(0, *t, 1)); 260 return (int) (t[0] > 0 && t[0] < 1); 261 } 262 } 263 // fall through if no t value found 264 case SkCubicType::kSerpentine: 265 case SkCubicType::kLocalCusp: 266 case SkCubicType::kCuspAtInfinity: { 267 double inflectionTs[2]; 268 int infTCount = cubic.findInflections(inflectionTs); 269 double maxCurvature[3]; 270 int roots = cubic.findMaxCurvature(maxCurvature); 271 #if DEBUG_CUBIC_SPLIT 272 SkDebugf("%s\n", __FUNCTION__); 273 cubic.dump(); 274 for (int index = 0; index < infTCount; ++index) { 275 SkDebugf("inflectionsTs[%d]=%1.9g ", index, inflectionTs[index]); 276 SkDPoint pt = cubic.ptAtT(inflectionTs[index]); 277 SkDVector dPt = cubic.dxdyAtT(inflectionTs[index]); 278 SkDLine perp = {{pt - dPt, pt + dPt}}; 279 perp.dump(); 280 } 281 for (int index = 0; index < roots; ++index) { 282 SkDebugf("maxCurvature[%d]=%1.9g ", index, maxCurvature[index]); 283 SkDPoint pt = cubic.ptAtT(maxCurvature[index]); 284 SkDVector dPt = cubic.dxdyAtT(maxCurvature[index]); 285 SkDLine perp = {{pt - dPt, pt + dPt}}; 286 perp.dump(); 287 } 288 #endif 289 if (infTCount == 2) { 290 for (int index = 0; index < roots; ++index) { 291 if (between(inflectionTs[0], maxCurvature[index], inflectionTs[1])) { 292 t[0] = maxCurvature[index]; 293 return (int) (t[0] > 0 && t[0] < 1); 294 } 295 } 296 } else { 297 int resultCount = 0; 298 // FIXME: constant found through experimentation -- maybe there's a better way.... 299 double precision = cubic.calcPrecision() * 2; 300 for (int index = 0; index < roots; ++index) { 301 double testT = maxCurvature[index]; 302 if (0 >= testT || testT >= 1) { 303 continue; 304 } 305 // don't call dxdyAtT since we want (0,0) results 306 SkDVector dPt = { derivative_at_t(&cubic.fPts[0].fX, testT), 307 derivative_at_t(&cubic.fPts[0].fY, testT) }; 308 double dPtLen = dPt.length(); 309 if (dPtLen < precision) { 310 t[resultCount++] = testT; 311 } 312 } 313 if (!resultCount && infTCount == 1) { 314 t[0] = inflectionTs[0]; 315 resultCount = (int) (t[0] > 0 && t[0] < 1); 316 } 317 return resultCount; 318 } 319 } 320 default: 321 ; 322 } 323 return 0; 324 } 325 326 bool SkDCubic::monotonicInX() const { 327 return precisely_between(fPts[0].fX, fPts[1].fX, fPts[3].fX) 328 && precisely_between(fPts[0].fX, fPts[2].fX, fPts[3].fX); 329 } 330 331 bool SkDCubic::monotonicInY() const { 332 return precisely_between(fPts[0].fY, fPts[1].fY, fPts[3].fY) 333 && precisely_between(fPts[0].fY, fPts[2].fY, fPts[3].fY); 334 } 335 336 void SkDCubic::otherPts(int index, const SkDPoint* o1Pts[kPointCount - 1]) const { 337 int offset = (int) !SkToBool(index); 338 o1Pts[0] = &fPts[offset]; 339 o1Pts[1] = &fPts[++offset]; 340 o1Pts[2] = &fPts[++offset]; 341 } 342 343 int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept, 344 SearchAxis xAxis, double* validRoots) const { 345 extrema += findInflections(&extremeTs[extrema]); 346 extremeTs[extrema++] = 0; 347 extremeTs[extrema] = 1; 348 SkASSERT(extrema < 6); 349 SkTQSort(extremeTs, extremeTs + extrema); 350 int validCount = 0; 351 for (int index = 0; index < extrema; ) { 352 double min = extremeTs[index]; 353 double max = extremeTs[++index]; 354 if (min == max) { 355 continue; 356 } 357 double newT = binarySearch(min, max, axisIntercept, xAxis); 358 if (newT >= 0) { 359 if (validCount >= 3) { 360 return 0; 361 } 362 validRoots[validCount++] = newT; 363 } 364 } 365 return validCount; 366 } 367 368 // cubic roots 369 370 static const double PI = 3.141592653589793; 371 372 // from SkGeometry.cpp (and Numeric Solutions, 5.6) 373 int SkDCubic::RootsValidT(double A, double B, double C, double D, double t[3]) { 374 double s[3]; 375 int realRoots = RootsReal(A, B, C, D, s); 376 int foundRoots = SkDQuad::AddValidTs(s, realRoots, t); 377 for (int index = 0; index < realRoots; ++index) { 378 double tValue = s[index]; 379 if (!approximately_one_or_less(tValue) && between(1, tValue, 1.00005)) { 380 for (int idx2 = 0; idx2 < foundRoots; ++idx2) { 381 if (approximately_equal(t[idx2], 1)) { 382 goto nextRoot; 383 } 384 } 385 SkASSERT(foundRoots < 3); 386 t[foundRoots++] = 1; 387 } else if (!approximately_zero_or_more(tValue) && between(-0.00005, tValue, 0)) { 388 for (int idx2 = 0; idx2 < foundRoots; ++idx2) { 389 if (approximately_equal(t[idx2], 0)) { 390 goto nextRoot; 391 } 392 } 393 SkASSERT(foundRoots < 3); 394 t[foundRoots++] = 0; 395 } 396 nextRoot: 397 ; 398 } 399 return foundRoots; 400 } 401 402 int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) { 403 #ifdef SK_DEBUG 404 // create a string mathematica understands 405 // GDB set print repe 15 # if repeated digits is a bother 406 // set print elements 400 # if line doesn't fit 407 char str[1024]; 408 sk_bzero(str, sizeof(str)); 409 SK_SNPRINTF(str, sizeof(str), "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", 410 A, B, C, D); 411 SkPathOpsDebug::MathematicaIze(str, sizeof(str)); 412 #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA 413 SkDebugf("%s\n", str); 414 #endif 415 #endif 416 if (approximately_zero(A) 417 && approximately_zero_when_compared_to(A, B) 418 && approximately_zero_when_compared_to(A, C) 419 && approximately_zero_when_compared_to(A, D)) { // we're just a quadratic 420 return SkDQuad::RootsReal(B, C, D, s); 421 } 422 if (approximately_zero_when_compared_to(D, A) 423 && approximately_zero_when_compared_to(D, B) 424 && approximately_zero_when_compared_to(D, C)) { // 0 is one root 425 int num = SkDQuad::RootsReal(A, B, C, s); 426 for (int i = 0; i < num; ++i) { 427 if (approximately_zero(s[i])) { 428 return num; 429 } 430 } 431 s[num++] = 0; 432 return num; 433 } 434 if (approximately_zero(A + B + C + D)) { // 1 is one root 435 int num = SkDQuad::RootsReal(A, A + B, -D, s); 436 for (int i = 0; i < num; ++i) { 437 if (AlmostDequalUlps(s[i], 1)) { 438 return num; 439 } 440 } 441 s[num++] = 1; 442 return num; 443 } 444 double a, b, c; 445 { 446 double invA = 1 / A; 447 a = B * invA; 448 b = C * invA; 449 c = D * invA; 450 } 451 double a2 = a * a; 452 double Q = (a2 - b * 3) / 9; 453 double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; 454 double R2 = R * R; 455 double Q3 = Q * Q * Q; 456 double R2MinusQ3 = R2 - Q3; 457 double adiv3 = a / 3; 458 double r; 459 double* roots = s; 460 if (R2MinusQ3 < 0) { // we have 3 real roots 461 // the divide/root can, due to finite precisions, be slightly outside of -1...1 462 double theta = acos(SkTPin(R / sqrt(Q3), -1., 1.)); 463 double neg2RootQ = -2 * sqrt(Q); 464 465 r = neg2RootQ * cos(theta / 3) - adiv3; 466 *roots++ = r; 467 468 r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; 469 if (!AlmostDequalUlps(s[0], r)) { 470 *roots++ = r; 471 } 472 r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; 473 if (!AlmostDequalUlps(s[0], r) && (roots - s == 1 || !AlmostDequalUlps(s[1], r))) { 474 *roots++ = r; 475 } 476 } else { // we have 1 real root 477 double sqrtR2MinusQ3 = sqrt(R2MinusQ3); 478 double A = fabs(R) + sqrtR2MinusQ3; 479 A = SkDCubeRoot(A); 480 if (R > 0) { 481 A = -A; 482 } 483 if (A != 0) { 484 A += Q / A; 485 } 486 r = A - adiv3; 487 *roots++ = r; 488 if (AlmostDequalUlps((double) R2, (double) Q3)) { 489 r = -A / 2 - adiv3; 490 if (!AlmostDequalUlps(s[0], r)) { 491 *roots++ = r; 492 } 493 } 494 } 495 return static_cast<int>(roots - s); 496 } 497 498 // OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t? 499 SkDVector SkDCubic::dxdyAtT(double t) const { 500 SkDVector result = { derivative_at_t(&fPts[0].fX, t), derivative_at_t(&fPts[0].fY, t) }; 501 if (result.fX == 0 && result.fY == 0) { 502 if (t == 0) { 503 result = fPts[2] - fPts[0]; 504 } else if (t == 1) { 505 result = fPts[3] - fPts[1]; 506 } else { 507 // incomplete 508 SkDebugf("!c"); 509 } 510 if (result.fX == 0 && result.fY == 0 && zero_or_one(t)) { 511 result = fPts[3] - fPts[0]; 512 } 513 } 514 return result; 515 } 516 517 // OPTIMIZE? share code with formulate_F1DotF2 518 int SkDCubic::findInflections(double tValues[]) const { 519 double Ax = fPts[1].fX - fPts[0].fX; 520 double Ay = fPts[1].fY - fPts[0].fY; 521 double Bx = fPts[2].fX - 2 * fPts[1].fX + fPts[0].fX; 522 double By = fPts[2].fY - 2 * fPts[1].fY + fPts[0].fY; 523 double Cx = fPts[3].fX + 3 * (fPts[1].fX - fPts[2].fX) - fPts[0].fX; 524 double Cy = fPts[3].fY + 3 * (fPts[1].fY - fPts[2].fY) - fPts[0].fY; 525 return SkDQuad::RootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues); 526 } 527 528 static void formulate_F1DotF2(const double src[], double coeff[4]) { 529 double a = src[2] - src[0]; 530 double b = src[4] - 2 * src[2] + src[0]; 531 double c = src[6] + 3 * (src[2] - src[4]) - src[0]; 532 coeff[0] = c * c; 533 coeff[1] = 3 * b * c; 534 coeff[2] = 2 * b * b + c * a; 535 coeff[3] = a * b; 536 } 537 538 /** SkDCubic'(t) = At^2 + Bt + C, where 539 A = 3(-a + 3(b - c) + d) 540 B = 6(a - 2b + c) 541 C = 3(b - a) 542 Solve for t, keeping only those that fit between 0 < t < 1 543 */ 544 int SkDCubic::FindExtrema(const double src[], double tValues[2]) { 545 // we divide A,B,C by 3 to simplify 546 double a = src[0]; 547 double b = src[2]; 548 double c = src[4]; 549 double d = src[6]; 550 double A = d - a + 3 * (b - c); 551 double B = 2 * (a - b - b + c); 552 double C = b - a; 553 554 return SkDQuad::RootsValidT(A, B, C, tValues); 555 } 556 557 /* from SkGeometry.cpp 558 Looking for F' dot F'' == 0 559 560 A = b - a 561 B = c - 2b + a 562 C = d - 3c + 3b - a 563 564 F' = 3Ct^2 + 6Bt + 3A 565 F'' = 6Ct + 6B 566 567 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 568 */ 569 int SkDCubic::findMaxCurvature(double tValues[]) const { 570 double coeffX[4], coeffY[4]; 571 int i; 572 formulate_F1DotF2(&fPts[0].fX, coeffX); 573 formulate_F1DotF2(&fPts[0].fY, coeffY); 574 for (i = 0; i < 4; i++) { 575 coeffX[i] = coeffX[i] + coeffY[i]; 576 } 577 return RootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues); 578 } 579 580 SkDPoint SkDCubic::ptAtT(double t) const { 581 if (0 == t) { 582 return fPts[0]; 583 } 584 if (1 == t) { 585 return fPts[3]; 586 } 587 double one_t = 1 - t; 588 double one_t2 = one_t * one_t; 589 double a = one_t2 * one_t; 590 double b = 3 * one_t2 * t; 591 double t2 = t * t; 592 double c = 3 * one_t * t2; 593 double d = t2 * t; 594 SkDPoint result = {a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX + d * fPts[3].fX, 595 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY + d * fPts[3].fY}; 596 return result; 597 } 598 599 /* 600 Given a cubic c, t1, and t2, find a small cubic segment. 601 602 The new cubic is defined as points A, B, C, and D, where 603 s1 = 1 - t1 604 s2 = 1 - t2 605 A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1 606 D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2 607 608 We don't have B or C. So We define two equations to isolate them. 609 First, compute two reference T values 1/3 and 2/3 from t1 to t2: 610 611 c(at (2*t1 + t2)/3) == E 612 c(at (t1 + 2*t2)/3) == F 613 614 Next, compute where those values must be if we know the values of B and C: 615 616 _12 = A*2/3 + B*1/3 617 12_ = A*1/3 + B*2/3 618 _23 = B*2/3 + C*1/3 619 23_ = B*1/3 + C*2/3 620 _34 = C*2/3 + D*1/3 621 34_ = C*1/3 + D*2/3 622 _123 = (A*2/3 + B*1/3)*2/3 + (B*2/3 + C*1/3)*1/3 = A*4/9 + B*4/9 + C*1/9 623 123_ = (A*1/3 + B*2/3)*1/3 + (B*1/3 + C*2/3)*2/3 = A*1/9 + B*4/9 + C*4/9 624 _234 = (B*2/3 + C*1/3)*2/3 + (C*2/3 + D*1/3)*1/3 = B*4/9 + C*4/9 + D*1/9 625 234_ = (B*1/3 + C*2/3)*1/3 + (C*1/3 + D*2/3)*2/3 = B*1/9 + C*4/9 + D*4/9 626 _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3 627 = A*8/27 + B*12/27 + C*6/27 + D*1/27 628 = E 629 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3 630 = A*1/27 + B*6/27 + C*12/27 + D*8/27 631 = F 632 E*27 = A*8 + B*12 + C*6 + D 633 F*27 = A + B*6 + C*12 + D*8 634 635 Group the known values on one side: 636 637 M = E*27 - A*8 - D = B*12 + C* 6 638 N = F*27 - A - D*8 = B* 6 + C*12 639 M*2 - N = B*18 640 N*2 - M = C*18 641 B = (M*2 - N)/18 642 C = (N*2 - M)/18 643 */ 644 645 static double interp_cubic_coords(const double* src, double t) { 646 double ab = SkDInterp(src[0], src[2], t); 647 double bc = SkDInterp(src[2], src[4], t); 648 double cd = SkDInterp(src[4], src[6], t); 649 double abc = SkDInterp(ab, bc, t); 650 double bcd = SkDInterp(bc, cd, t); 651 double abcd = SkDInterp(abc, bcd, t); 652 return abcd; 653 } 654 655 SkDCubic SkDCubic::subDivide(double t1, double t2) const { 656 if (t1 == 0 || t2 == 1) { 657 if (t1 == 0 && t2 == 1) { 658 return *this; 659 } 660 SkDCubicPair pair = chopAt(t1 == 0 ? t2 : t1); 661 SkDCubic dst = t1 == 0 ? pair.first() : pair.second(); 662 return dst; 663 } 664 SkDCubic dst; 665 double ax = dst[0].fX = interp_cubic_coords(&fPts[0].fX, t1); 666 double ay = dst[0].fY = interp_cubic_coords(&fPts[0].fY, t1); 667 double ex = interp_cubic_coords(&fPts[0].fX, (t1*2+t2)/3); 668 double ey = interp_cubic_coords(&fPts[0].fY, (t1*2+t2)/3); 669 double fx = interp_cubic_coords(&fPts[0].fX, (t1+t2*2)/3); 670 double fy = interp_cubic_coords(&fPts[0].fY, (t1+t2*2)/3); 671 double dx = dst[3].fX = interp_cubic_coords(&fPts[0].fX, t2); 672 double dy = dst[3].fY = interp_cubic_coords(&fPts[0].fY, t2); 673 double mx = ex * 27 - ax * 8 - dx; 674 double my = ey * 27 - ay * 8 - dy; 675 double nx = fx * 27 - ax - dx * 8; 676 double ny = fy * 27 - ay - dy * 8; 677 /* bx = */ dst[1].fX = (mx * 2 - nx) / 18; 678 /* by = */ dst[1].fY = (my * 2 - ny) / 18; 679 /* cx = */ dst[2].fX = (nx * 2 - mx) / 18; 680 /* cy = */ dst[2].fY = (ny * 2 - my) / 18; 681 // FIXME: call align() ? 682 return dst; 683 } 684 685 void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d, 686 double t1, double t2, SkDPoint dst[2]) const { 687 SkASSERT(t1 != t2); 688 // this approach assumes that the control points computed directly are accurate enough 689 SkDCubic sub = subDivide(t1, t2); 690 dst[0] = sub[1] + (a - sub[0]); 691 dst[1] = sub[2] + (d - sub[3]); 692 if (t1 == 0 || t2 == 0) { 693 align(0, 1, t1 == 0 ? &dst[0] : &dst[1]); 694 } 695 if (t1 == 1 || t2 == 1) { 696 align(3, 2, t1 == 1 ? &dst[0] : &dst[1]); 697 } 698 if (AlmostBequalUlps(dst[0].fX, a.fX)) { 699 dst[0].fX = a.fX; 700 } 701 if (AlmostBequalUlps(dst[0].fY, a.fY)) { 702 dst[0].fY = a.fY; 703 } 704 if (AlmostBequalUlps(dst[1].fX, d.fX)) { 705 dst[1].fX = d.fX; 706 } 707 if (AlmostBequalUlps(dst[1].fY, d.fY)) { 708 dst[1].fY = d.fY; 709 } 710 } 711 712 bool SkDCubic::toFloatPoints(SkPoint* pts) const { 713 const double* dCubic = &fPts[0].fX; 714 SkScalar* cubic = &pts[0].fX; 715 for (int index = 0; index < kPointCount * 2; ++index) { 716 cubic[index] = SkDoubleToScalar(dCubic[index]); 717 if (SkScalarAbs(cubic[index]) < FLT_EPSILON_ORDERABLE_ERR) { 718 cubic[index] = 0; 719 } 720 } 721 return SkScalarsAreFinite(&pts->fX, kPointCount * 2); 722 } 723 724 double SkDCubic::top(const SkDCubic& dCurve, double startT, double endT, SkDPoint*topPt) const { 725 double extremeTs[2]; 726 double topT = -1; 727 int roots = SkDCubic::FindExtrema(&fPts[0].fY, extremeTs); 728 for (int index = 0; index < roots; ++index) { 729 double t = startT + (endT - startT) * extremeTs[index]; 730 SkDPoint mid = dCurve.ptAtT(t); 731 if (topPt->fY > mid.fY || (topPt->fY == mid.fY && topPt->fX > mid.fX)) { 732 topT = t; 733 *topPt = mid; 734 } 735 } 736 return topT; 737 } 738