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 SkScalar d[3]; 252 SkCubicType cubicType = SkClassifyCubic(pointsPtr, d); 253 switch (cubicType) { 254 case kLoop_SkCubicType: { 255 // crib code from gpu path utils that finds t values where loop self-intersects 256 // use it to find mid of t values which should be a friendly place to chop 257 SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]); 258 SkScalar ls = d[1] - tempSqrt; 259 SkScalar lt = 2.f * d[0]; 260 SkScalar ms = d[1] + tempSqrt; 261 SkScalar mt = 2.f * d[0]; 262 if (roughly_between(0, ls, lt) && roughly_between(0, ms, mt)) { 263 ls = ls / lt; 264 ms = ms / mt; 265 SkASSERT(roughly_between(0, ls, 1) && roughly_between(0, ms, 1)); 266 t[0] = (ls + ms) / 2; 267 SkASSERT(roughly_between(0, *t, 1)); 268 return (int) (t[0] > 0 && t[0] < 1); 269 } 270 } 271 // fall through if no t value found 272 case kSerpentine_SkCubicType: 273 case kCusp_SkCubicType: { 274 double inflectionTs[2]; 275 int infTCount = cubic.findInflections(inflectionTs); 276 double maxCurvature[3]; 277 int roots = cubic.findMaxCurvature(maxCurvature); 278 #if DEBUG_CUBIC_SPLIT 279 SkDebugf("%s\n", __FUNCTION__); 280 cubic.dump(); 281 for (int index = 0; index < infTCount; ++index) { 282 SkDebugf("inflectionsTs[%d]=%1.9g ", index, inflectionTs[index]); 283 SkDPoint pt = cubic.ptAtT(inflectionTs[index]); 284 SkDVector dPt = cubic.dxdyAtT(inflectionTs[index]); 285 SkDLine perp = {{pt - dPt, pt + dPt}}; 286 perp.dump(); 287 } 288 for (int index = 0; index < roots; ++index) { 289 SkDebugf("maxCurvature[%d]=%1.9g ", index, maxCurvature[index]); 290 SkDPoint pt = cubic.ptAtT(maxCurvature[index]); 291 SkDVector dPt = cubic.dxdyAtT(maxCurvature[index]); 292 SkDLine perp = {{pt - dPt, pt + dPt}}; 293 perp.dump(); 294 } 295 #endif 296 if (infTCount == 2) { 297 for (int index = 0; index < roots; ++index) { 298 if (between(inflectionTs[0], maxCurvature[index], inflectionTs[1])) { 299 t[0] = maxCurvature[index]; 300 return (int) (t[0] > 0 && t[0] < 1); 301 } 302 } 303 } else { 304 int resultCount = 0; 305 // FIXME: constant found through experimentation -- maybe there's a better way.... 306 double precision = cubic.calcPrecision() * 2; 307 for (int index = 0; index < roots; ++index) { 308 double testT = maxCurvature[index]; 309 if (0 >= testT || testT >= 1) { 310 continue; 311 } 312 // don't call dxdyAtT since we want (0,0) results 313 SkDVector dPt = { derivative_at_t(&cubic.fPts[0].fX, testT), 314 derivative_at_t(&cubic.fPts[0].fY, testT) }; 315 double dPtLen = dPt.length(); 316 if (dPtLen < precision) { 317 t[resultCount++] = testT; 318 } 319 } 320 if (!resultCount && infTCount == 1) { 321 t[0] = inflectionTs[0]; 322 resultCount = (int) (t[0] > 0 && t[0] < 1); 323 } 324 return resultCount; 325 } 326 } 327 default: 328 ; 329 } 330 return 0; 331 } 332 333 bool SkDCubic::monotonicInX() const { 334 return precisely_between(fPts[0].fX, fPts[1].fX, fPts[3].fX) 335 && precisely_between(fPts[0].fX, fPts[2].fX, fPts[3].fX); 336 } 337 338 bool SkDCubic::monotonicInY() const { 339 return precisely_between(fPts[0].fY, fPts[1].fY, fPts[3].fY) 340 && precisely_between(fPts[0].fY, fPts[2].fY, fPts[3].fY); 341 } 342 343 void SkDCubic::otherPts(int index, const SkDPoint* o1Pts[kPointCount - 1]) const { 344 int offset = (int) !SkToBool(index); 345 o1Pts[0] = &fPts[offset]; 346 o1Pts[1] = &fPts[++offset]; 347 o1Pts[2] = &fPts[++offset]; 348 } 349 350 int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept, 351 SearchAxis xAxis, double* validRoots) const { 352 extrema += findInflections(&extremeTs[extrema]); 353 extremeTs[extrema++] = 0; 354 extremeTs[extrema] = 1; 355 SkASSERT(extrema < 6); 356 SkTQSort(extremeTs, extremeTs + extrema); 357 int validCount = 0; 358 for (int index = 0; index < extrema; ) { 359 double min = extremeTs[index]; 360 double max = extremeTs[++index]; 361 if (min == max) { 362 continue; 363 } 364 double newT = binarySearch(min, max, axisIntercept, xAxis); 365 if (newT >= 0) { 366 if (validCount >= 3) { 367 return 0; 368 } 369 validRoots[validCount++] = newT; 370 } 371 } 372 return validCount; 373 } 374 375 // cubic roots 376 377 static const double PI = 3.141592653589793; 378 379 // from SkGeometry.cpp (and Numeric Solutions, 5.6) 380 int SkDCubic::RootsValidT(double A, double B, double C, double D, double t[3]) { 381 double s[3]; 382 int realRoots = RootsReal(A, B, C, D, s); 383 int foundRoots = SkDQuad::AddValidTs(s, realRoots, t); 384 for (int index = 0; index < realRoots; ++index) { 385 double tValue = s[index]; 386 if (!approximately_one_or_less(tValue) && between(1, tValue, 1.00005)) { 387 for (int idx2 = 0; idx2 < foundRoots; ++idx2) { 388 if (approximately_equal(t[idx2], 1)) { 389 goto nextRoot; 390 } 391 } 392 SkASSERT(foundRoots < 3); 393 t[foundRoots++] = 1; 394 } else if (!approximately_zero_or_more(tValue) && between(-0.00005, tValue, 0)) { 395 for (int idx2 = 0; idx2 < foundRoots; ++idx2) { 396 if (approximately_equal(t[idx2], 0)) { 397 goto nextRoot; 398 } 399 } 400 SkASSERT(foundRoots < 3); 401 t[foundRoots++] = 0; 402 } 403 nextRoot: 404 ; 405 } 406 return foundRoots; 407 } 408 409 int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) { 410 #ifdef SK_DEBUG 411 // create a string mathematica understands 412 // GDB set print repe 15 # if repeated digits is a bother 413 // set print elements 400 # if line doesn't fit 414 char str[1024]; 415 sk_bzero(str, sizeof(str)); 416 SK_SNPRINTF(str, sizeof(str), "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", 417 A, B, C, D); 418 SkPathOpsDebug::MathematicaIze(str, sizeof(str)); 419 #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA 420 SkDebugf("%s\n", str); 421 #endif 422 #endif 423 if (approximately_zero(A) 424 && approximately_zero_when_compared_to(A, B) 425 && approximately_zero_when_compared_to(A, C) 426 && approximately_zero_when_compared_to(A, D)) { // we're just a quadratic 427 return SkDQuad::RootsReal(B, C, D, s); 428 } 429 if (approximately_zero_when_compared_to(D, A) 430 && approximately_zero_when_compared_to(D, B) 431 && approximately_zero_when_compared_to(D, C)) { // 0 is one root 432 int num = SkDQuad::RootsReal(A, B, C, s); 433 for (int i = 0; i < num; ++i) { 434 if (approximately_zero(s[i])) { 435 return num; 436 } 437 } 438 s[num++] = 0; 439 return num; 440 } 441 if (approximately_zero(A + B + C + D)) { // 1 is one root 442 int num = SkDQuad::RootsReal(A, A + B, -D, s); 443 for (int i = 0; i < num; ++i) { 444 if (AlmostDequalUlps(s[i], 1)) { 445 return num; 446 } 447 } 448 s[num++] = 1; 449 return num; 450 } 451 double a, b, c; 452 { 453 double invA = 1 / A; 454 a = B * invA; 455 b = C * invA; 456 c = D * invA; 457 } 458 double a2 = a * a; 459 double Q = (a2 - b * 3) / 9; 460 double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; 461 double R2 = R * R; 462 double Q3 = Q * Q * Q; 463 double R2MinusQ3 = R2 - Q3; 464 double adiv3 = a / 3; 465 double r; 466 double* roots = s; 467 if (R2MinusQ3 < 0) { // we have 3 real roots 468 // the divide/root can, due to finite precisions, be slightly outside of -1...1 469 double theta = acos(SkTPin(R / sqrt(Q3), -1., 1.)); 470 double neg2RootQ = -2 * sqrt(Q); 471 472 r = neg2RootQ * cos(theta / 3) - adiv3; 473 *roots++ = r; 474 475 r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; 476 if (!AlmostDequalUlps(s[0], r)) { 477 *roots++ = r; 478 } 479 r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; 480 if (!AlmostDequalUlps(s[0], r) && (roots - s == 1 || !AlmostDequalUlps(s[1], r))) { 481 *roots++ = r; 482 } 483 } else { // we have 1 real root 484 double sqrtR2MinusQ3 = sqrt(R2MinusQ3); 485 double A = fabs(R) + sqrtR2MinusQ3; 486 A = SkDCubeRoot(A); 487 if (R > 0) { 488 A = -A; 489 } 490 if (A != 0) { 491 A += Q / A; 492 } 493 r = A - adiv3; 494 *roots++ = r; 495 if (AlmostDequalUlps((double) R2, (double) Q3)) { 496 r = -A / 2 - adiv3; 497 if (!AlmostDequalUlps(s[0], r)) { 498 *roots++ = r; 499 } 500 } 501 } 502 return static_cast<int>(roots - s); 503 } 504 505 // OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t? 506 SkDVector SkDCubic::dxdyAtT(double t) const { 507 SkDVector result = { derivative_at_t(&fPts[0].fX, t), derivative_at_t(&fPts[0].fY, t) }; 508 if (result.fX == 0 && result.fY == 0) { 509 if (t == 0) { 510 result = fPts[2] - fPts[0]; 511 } else if (t == 1) { 512 result = fPts[3] - fPts[1]; 513 } else { 514 // incomplete 515 SkDebugf("!c"); 516 } 517 if (result.fX == 0 && result.fY == 0 && zero_or_one(t)) { 518 result = fPts[3] - fPts[0]; 519 } 520 } 521 return result; 522 } 523 524 // OPTIMIZE? share code with formulate_F1DotF2 525 int SkDCubic::findInflections(double tValues[]) const { 526 double Ax = fPts[1].fX - fPts[0].fX; 527 double Ay = fPts[1].fY - fPts[0].fY; 528 double Bx = fPts[2].fX - 2 * fPts[1].fX + fPts[0].fX; 529 double By = fPts[2].fY - 2 * fPts[1].fY + fPts[0].fY; 530 double Cx = fPts[3].fX + 3 * (fPts[1].fX - fPts[2].fX) - fPts[0].fX; 531 double Cy = fPts[3].fY + 3 * (fPts[1].fY - fPts[2].fY) - fPts[0].fY; 532 return SkDQuad::RootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues); 533 } 534 535 static void formulate_F1DotF2(const double src[], double coeff[4]) { 536 double a = src[2] - src[0]; 537 double b = src[4] - 2 * src[2] + src[0]; 538 double c = src[6] + 3 * (src[2] - src[4]) - src[0]; 539 coeff[0] = c * c; 540 coeff[1] = 3 * b * c; 541 coeff[2] = 2 * b * b + c * a; 542 coeff[3] = a * b; 543 } 544 545 /** SkDCubic'(t) = At^2 + Bt + C, where 546 A = 3(-a + 3(b - c) + d) 547 B = 6(a - 2b + c) 548 C = 3(b - a) 549 Solve for t, keeping only those that fit between 0 < t < 1 550 */ 551 int SkDCubic::FindExtrema(const double src[], double tValues[2]) { 552 // we divide A,B,C by 3 to simplify 553 double a = src[0]; 554 double b = src[2]; 555 double c = src[4]; 556 double d = src[6]; 557 double A = d - a + 3 * (b - c); 558 double B = 2 * (a - b - b + c); 559 double C = b - a; 560 561 return SkDQuad::RootsValidT(A, B, C, tValues); 562 } 563 564 /* from SkGeometry.cpp 565 Looking for F' dot F'' == 0 566 567 A = b - a 568 B = c - 2b + a 569 C = d - 3c + 3b - a 570 571 F' = 3Ct^2 + 6Bt + 3A 572 F'' = 6Ct + 6B 573 574 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 575 */ 576 int SkDCubic::findMaxCurvature(double tValues[]) const { 577 double coeffX[4], coeffY[4]; 578 int i; 579 formulate_F1DotF2(&fPts[0].fX, coeffX); 580 formulate_F1DotF2(&fPts[0].fY, coeffY); 581 for (i = 0; i < 4; i++) { 582 coeffX[i] = coeffX[i] + coeffY[i]; 583 } 584 return RootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues); 585 } 586 587 SkDPoint SkDCubic::ptAtT(double t) const { 588 if (0 == t) { 589 return fPts[0]; 590 } 591 if (1 == t) { 592 return fPts[3]; 593 } 594 double one_t = 1 - t; 595 double one_t2 = one_t * one_t; 596 double a = one_t2 * one_t; 597 double b = 3 * one_t2 * t; 598 double t2 = t * t; 599 double c = 3 * one_t * t2; 600 double d = t2 * t; 601 SkDPoint result = {a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX + d * fPts[3].fX, 602 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY + d * fPts[3].fY}; 603 return result; 604 } 605 606 /* 607 Given a cubic c, t1, and t2, find a small cubic segment. 608 609 The new cubic is defined as points A, B, C, and D, where 610 s1 = 1 - t1 611 s2 = 1 - t2 612 A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1 613 D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2 614 615 We don't have B or C. So We define two equations to isolate them. 616 First, compute two reference T values 1/3 and 2/3 from t1 to t2: 617 618 c(at (2*t1 + t2)/3) == E 619 c(at (t1 + 2*t2)/3) == F 620 621 Next, compute where those values must be if we know the values of B and C: 622 623 _12 = A*2/3 + B*1/3 624 12_ = A*1/3 + B*2/3 625 _23 = B*2/3 + C*1/3 626 23_ = B*1/3 + C*2/3 627 _34 = C*2/3 + D*1/3 628 34_ = C*1/3 + D*2/3 629 _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 630 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 631 _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 632 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 633 _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3 634 = A*8/27 + B*12/27 + C*6/27 + D*1/27 635 = E 636 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3 637 = A*1/27 + B*6/27 + C*12/27 + D*8/27 638 = F 639 E*27 = A*8 + B*12 + C*6 + D 640 F*27 = A + B*6 + C*12 + D*8 641 642 Group the known values on one side: 643 644 M = E*27 - A*8 - D = B*12 + C* 6 645 N = F*27 - A - D*8 = B* 6 + C*12 646 M*2 - N = B*18 647 N*2 - M = C*18 648 B = (M*2 - N)/18 649 C = (N*2 - M)/18 650 */ 651 652 static double interp_cubic_coords(const double* src, double t) { 653 double ab = SkDInterp(src[0], src[2], t); 654 double bc = SkDInterp(src[2], src[4], t); 655 double cd = SkDInterp(src[4], src[6], t); 656 double abc = SkDInterp(ab, bc, t); 657 double bcd = SkDInterp(bc, cd, t); 658 double abcd = SkDInterp(abc, bcd, t); 659 return abcd; 660 } 661 662 SkDCubic SkDCubic::subDivide(double t1, double t2) const { 663 if (t1 == 0 || t2 == 1) { 664 if (t1 == 0 && t2 == 1) { 665 return *this; 666 } 667 SkDCubicPair pair = chopAt(t1 == 0 ? t2 : t1); 668 SkDCubic dst = t1 == 0 ? pair.first() : pair.second(); 669 return dst; 670 } 671 SkDCubic dst; 672 double ax = dst[0].fX = interp_cubic_coords(&fPts[0].fX, t1); 673 double ay = dst[0].fY = interp_cubic_coords(&fPts[0].fY, t1); 674 double ex = interp_cubic_coords(&fPts[0].fX, (t1*2+t2)/3); 675 double ey = interp_cubic_coords(&fPts[0].fY, (t1*2+t2)/3); 676 double fx = interp_cubic_coords(&fPts[0].fX, (t1+t2*2)/3); 677 double fy = interp_cubic_coords(&fPts[0].fY, (t1+t2*2)/3); 678 double dx = dst[3].fX = interp_cubic_coords(&fPts[0].fX, t2); 679 double dy = dst[3].fY = interp_cubic_coords(&fPts[0].fY, t2); 680 double mx = ex * 27 - ax * 8 - dx; 681 double my = ey * 27 - ay * 8 - dy; 682 double nx = fx * 27 - ax - dx * 8; 683 double ny = fy * 27 - ay - dy * 8; 684 /* bx = */ dst[1].fX = (mx * 2 - nx) / 18; 685 /* by = */ dst[1].fY = (my * 2 - ny) / 18; 686 /* cx = */ dst[2].fX = (nx * 2 - mx) / 18; 687 /* cy = */ dst[2].fY = (ny * 2 - my) / 18; 688 // FIXME: call align() ? 689 return dst; 690 } 691 692 void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d, 693 double t1, double t2, SkDPoint dst[2]) const { 694 SkASSERT(t1 != t2); 695 // this approach assumes that the control points computed directly are accurate enough 696 SkDCubic sub = subDivide(t1, t2); 697 dst[0] = sub[1] + (a - sub[0]); 698 dst[1] = sub[2] + (d - sub[3]); 699 if (t1 == 0 || t2 == 0) { 700 align(0, 1, t1 == 0 ? &dst[0] : &dst[1]); 701 } 702 if (t1 == 1 || t2 == 1) { 703 align(3, 2, t1 == 1 ? &dst[0] : &dst[1]); 704 } 705 if (AlmostBequalUlps(dst[0].fX, a.fX)) { 706 dst[0].fX = a.fX; 707 } 708 if (AlmostBequalUlps(dst[0].fY, a.fY)) { 709 dst[0].fY = a.fY; 710 } 711 if (AlmostBequalUlps(dst[1].fX, d.fX)) { 712 dst[1].fX = d.fX; 713 } 714 if (AlmostBequalUlps(dst[1].fY, d.fY)) { 715 dst[1].fY = d.fY; 716 } 717 } 718 719 bool SkDCubic::toFloatPoints(SkPoint* pts) const { 720 const double* dCubic = &fPts[0].fX; 721 SkScalar* cubic = &pts[0].fX; 722 for (int index = 0; index < kPointCount * 2; ++index) { 723 *cubic++ = SkDoubleToScalar(*dCubic++); 724 } 725 return SkScalarsAreFinite(&pts->fX, kPointCount * 2); 726 } 727 728 double SkDCubic::top(const SkDCubic& dCurve, double startT, double endT, SkDPoint*topPt) const { 729 double extremeTs[2]; 730 double topT = -1; 731 int roots = SkDCubic::FindExtrema(&fPts[0].fY, extremeTs); 732 for (int index = 0; index < roots; ++index) { 733 double t = startT + (endT - startT) * extremeTs[index]; 734 SkDPoint mid = dCurve.ptAtT(t); 735 if (topPt->fY > mid.fY || (topPt->fY == mid.fY && topPt->fX > mid.fX)) { 736 topT = t; 737 *topPt = mid; 738 } 739 } 740 return topT; 741 } 742
