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 17const int SkDCubic::gPrecisionUnit = 256; // FIXME: test different values in test framework 18 19void 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 30double 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 SkASSERT(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// FIXME: cache keep the bounds and/or precision with the caller? 79double SkDCubic::calcPrecision() const { 80 SkDRect dRect; 81 dRect.setBounds(*this); // OPTIMIZATION: just use setRawBounds ? 82 double width = dRect.fRight - dRect.fLeft; 83 double height = dRect.fBottom - dRect.fTop; 84 return (width > height ? width : height) / gPrecisionUnit; 85} 86 87 88/* classic one t subdivision */ 89static void interp_cubic_coords(const double* src, double* dst, double t) { 90 double ab = SkDInterp(src[0], src[2], t); 91 double bc = SkDInterp(src[2], src[4], t); 92 double cd = SkDInterp(src[4], src[6], t); 93 double abc = SkDInterp(ab, bc, t); 94 double bcd = SkDInterp(bc, cd, t); 95 double abcd = SkDInterp(abc, bcd, t); 96 97 dst[0] = src[0]; 98 dst[2] = ab; 99 dst[4] = abc; 100 dst[6] = abcd; 101 dst[8] = bcd; 102 dst[10] = cd; 103 dst[12] = src[6]; 104} 105 106SkDCubicPair SkDCubic::chopAt(double t) const { 107 SkDCubicPair dst; 108 if (t == 0.5) { 109 dst.pts[0] = fPts[0]; 110 dst.pts[1].fX = (fPts[0].fX + fPts[1].fX) / 2; 111 dst.pts[1].fY = (fPts[0].fY + fPts[1].fY) / 2; 112 dst.pts[2].fX = (fPts[0].fX + 2 * fPts[1].fX + fPts[2].fX) / 4; 113 dst.pts[2].fY = (fPts[0].fY + 2 * fPts[1].fY + fPts[2].fY) / 4; 114 dst.pts[3].fX = (fPts[0].fX + 3 * (fPts[1].fX + fPts[2].fX) + fPts[3].fX) / 8; 115 dst.pts[3].fY = (fPts[0].fY + 3 * (fPts[1].fY + fPts[2].fY) + fPts[3].fY) / 8; 116 dst.pts[4].fX = (fPts[1].fX + 2 * fPts[2].fX + fPts[3].fX) / 4; 117 dst.pts[4].fY = (fPts[1].fY + 2 * fPts[2].fY + fPts[3].fY) / 4; 118 dst.pts[5].fX = (fPts[2].fX + fPts[3].fX) / 2; 119 dst.pts[5].fY = (fPts[2].fY + fPts[3].fY) / 2; 120 dst.pts[6] = fPts[3]; 121 return dst; 122 } 123 interp_cubic_coords(&fPts[0].fX, &dst.pts[0].fX, t); 124 interp_cubic_coords(&fPts[0].fY, &dst.pts[0].fY, t); 125 return dst; 126} 127 128void SkDCubic::Coefficients(const double* src, double* A, double* B, double* C, double* D) { 129 *A = src[6]; // d 130 *B = src[4] * 3; // 3*c 131 *C = src[2] * 3; // 3*b 132 *D = src[0]; // a 133 *A -= *D - *C + *B; // A = -a + 3*b - 3*c + d 134 *B += 3 * *D - 2 * *C; // B = 3*a - 6*b + 3*c 135 *C -= 3 * *D; // C = -3*a + 3*b 136} 137 138bool SkDCubic::endsAreExtremaInXOrY() const { 139 return (between(fPts[0].fX, fPts[1].fX, fPts[3].fX) 140 && between(fPts[0].fX, fPts[2].fX, fPts[3].fX)) 141 || (between(fPts[0].fY, fPts[1].fY, fPts[3].fY) 142 && between(fPts[0].fY, fPts[2].fY, fPts[3].fY)); 143} 144 145// Do a quick reject by rotating all points relative to a line formed by 146// a pair of one cubic's points. If the 2nd cubic's points 147// are on the line or on the opposite side from the 1st cubic's 'odd man', the 148// curves at most intersect at the endpoints. 149/* if returning true, check contains true if cubic's hull collapsed, making the cubic linear 150 if returning false, check contains true if the the cubic pair have only the end point in common 151*/ 152bool SkDCubic::hullIntersects(const SkDPoint* pts, int ptCount, bool* isLinear) const { 153 bool linear = true; 154 char hullOrder[4]; 155 int hullCount = convexHull(hullOrder); 156 int end1 = hullOrder[0]; 157 int hullIndex = 0; 158 const SkDPoint* endPt[2]; 159 endPt[0] = &fPts[end1]; 160 do { 161 hullIndex = (hullIndex + 1) % hullCount; 162 int end2 = hullOrder[hullIndex]; 163 endPt[1] = &fPts[end2]; 164 double origX = endPt[0]->fX; 165 double origY = endPt[0]->fY; 166 double adj = endPt[1]->fX - origX; 167 double opp = endPt[1]->fY - origY; 168 int oddManMask = other_two(end1, end2); 169 int oddMan = end1 ^ oddManMask; 170 double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp; 171 int oddMan2 = end2 ^ oddManMask; 172 double sign2 = (fPts[oddMan2].fY - origY) * adj - (fPts[oddMan2].fX - origX) * opp; 173 if (sign * sign2 < 0) { 174 continue; 175 } 176 if (approximately_zero(sign)) { 177 sign = sign2; 178 if (approximately_zero(sign)) { 179 continue; 180 } 181 } 182 linear = false; 183 bool foundOutlier = false; 184 for (int n = 0; n < ptCount; ++n) { 185 double test = (pts[n].fY - origY) * adj - (pts[n].fX - origX) * opp; 186 if (test * sign > 0 && !precisely_zero(test)) { 187 foundOutlier = true; 188 break; 189 } 190 } 191 if (!foundOutlier) { 192 return false; 193 } 194 endPt[0] = endPt[1]; 195 end1 = end2; 196 } while (hullIndex); 197 *isLinear = linear; 198 return true; 199} 200 201bool SkDCubic::hullIntersects(const SkDCubic& c2, bool* isLinear) const { 202 return hullIntersects(c2.fPts, c2.kPointCount, isLinear); 203} 204 205bool SkDCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const { 206 return hullIntersects(quad.fPts, quad.kPointCount, isLinear); 207} 208 209bool SkDCubic::hullIntersects(const SkDConic& conic, bool* isLinear) const { 210 211 return hullIntersects(conic.fPts, isLinear); 212} 213 214bool SkDCubic::isLinear(int startIndex, int endIndex) const { 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 232bool SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t, CubicType* resultType) { 233 SkScalar d[3]; 234 SkCubicType cubicType = SkClassifyCubic(pointsPtr, d); 235 if (cubicType == kLoop_SkCubicType) { 236 // crib code from gpu path utils that finds t values where loop self-intersects 237 // use it to find mid of t values which should be a friendly place to chop 238 SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]); 239 SkScalar ls = d[1] - tempSqrt; 240 SkScalar lt = 2.f * d[0]; 241 SkScalar ms = d[1] + tempSqrt; 242 SkScalar mt = 2.f * d[0]; 243 if (between(0, ls, lt) || between(0, ms, mt)) { 244 ls = ls / lt; 245 ms = ms / mt; 246 SkScalar smaller = SkTMax(0.f, SkTMin(ls, ms)); 247 SkScalar larger = SkTMin(1.f, SkTMax(ls, ms)); 248 *t = (smaller + larger) / 2; 249 *resultType = kSplitAtLoop_SkDCubicType; 250 return *t > 0 && *t < 1; 251 } 252 } else if (kSerpentine_SkCubicType == cubicType || kCusp_SkCubicType == cubicType) { 253 SkDCubic cubic; 254 cubic.set(pointsPtr); 255 double inflectionTs[2]; 256 int infTCount = cubic.findInflections(inflectionTs); 257 if (infTCount == 2) { 258 double maxCurvature[3]; 259 int roots = cubic.findMaxCurvature(maxCurvature); 260#if DEBUG_CUBIC_SPLIT 261 SkDebugf("%s\n", __FUNCTION__); 262 cubic.dump(); 263 for (int index = 0; index < infTCount; ++index) { 264 SkDebugf("inflectionsTs[%d]=%1.9g ", index, inflectionTs[index]); 265 SkDPoint pt = cubic.ptAtT(inflectionTs[index]); 266 SkDVector dPt = cubic.dxdyAtT(inflectionTs[index]); 267 SkDLine perp = {{pt - dPt, pt + dPt}}; 268 perp.dump(); 269 } 270 for (int index = 0; index < roots; ++index) { 271 SkDebugf("maxCurvature[%d]=%1.9g ", index, maxCurvature[index]); 272 SkDPoint pt = cubic.ptAtT(maxCurvature[index]); 273 SkDVector dPt = cubic.dxdyAtT(maxCurvature[index]); 274 SkDLine perp = {{pt - dPt, pt + dPt}}; 275 perp.dump(); 276 } 277#endif 278 for (int index = 0; index < roots; ++index) { 279 if (between(inflectionTs[0], maxCurvature[index], inflectionTs[1])) { 280 *t = maxCurvature[index]; 281 *resultType = kSplitAtMaxCurvature_SkDCubicType; 282 return true; 283 } 284 } 285 } else if (infTCount == 1) { 286 *t = inflectionTs[0]; 287 *resultType = kSplitAtInflection_SkDCubicType; 288 return *t > 0 && *t < 1; 289 } 290 } 291 return false; 292} 293 294bool SkDCubic::monotonicInX() const { 295 return precisely_between(fPts[0].fX, fPts[1].fX, fPts[3].fX) 296 && precisely_between(fPts[0].fX, fPts[2].fX, fPts[3].fX); 297} 298 299bool SkDCubic::monotonicInY() const { 300 return precisely_between(fPts[0].fY, fPts[1].fY, fPts[3].fY) 301 && precisely_between(fPts[0].fY, fPts[2].fY, fPts[3].fY); 302} 303 304void SkDCubic::otherPts(int index, const SkDPoint* o1Pts[kPointCount - 1]) const { 305 int offset = (int) !SkToBool(index); 306 o1Pts[0] = &fPts[offset]; 307 o1Pts[1] = &fPts[++offset]; 308 o1Pts[2] = &fPts[++offset]; 309} 310 311int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept, 312 SearchAxis xAxis, double* validRoots) const { 313 extrema += findInflections(&extremeTs[extrema]); 314 extremeTs[extrema++] = 0; 315 extremeTs[extrema] = 1; 316 SkTQSort(extremeTs, extremeTs + extrema); 317 int validCount = 0; 318 for (int index = 0; index < extrema; ) { 319 double min = extremeTs[index]; 320 double max = extremeTs[++index]; 321 if (min == max) { 322 continue; 323 } 324 double newT = binarySearch(min, max, axisIntercept, xAxis); 325 if (newT >= 0) { 326 validRoots[validCount++] = newT; 327 } 328 } 329 return validCount; 330} 331 332// cubic roots 333 334static const double PI = 3.141592653589793; 335 336// from SkGeometry.cpp (and Numeric Solutions, 5.6) 337int SkDCubic::RootsValidT(double A, double B, double C, double D, double t[3]) { 338 double s[3]; 339 int realRoots = RootsReal(A, B, C, D, s); 340 int foundRoots = SkDQuad::AddValidTs(s, realRoots, t); 341 for (int index = 0; index < realRoots; ++index) { 342 double tValue = s[index]; 343 if (!approximately_one_or_less(tValue) && between(1, tValue, 1.00005)) { 344 for (int idx2 = 0; idx2 < foundRoots; ++idx2) { 345 if (approximately_equal(t[idx2], 1)) { 346 goto nextRoot; 347 } 348 } 349 SkASSERT(foundRoots < 3); 350 t[foundRoots++] = 1; 351 } else if (!approximately_zero_or_more(tValue) && between(-0.00005, tValue, 0)) { 352 for (int idx2 = 0; idx2 < foundRoots; ++idx2) { 353 if (approximately_equal(t[idx2], 0)) { 354 goto nextRoot; 355 } 356 } 357 SkASSERT(foundRoots < 3); 358 t[foundRoots++] = 0; 359 } 360nextRoot: 361 ; 362 } 363 return foundRoots; 364} 365 366int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) { 367#ifdef SK_DEBUG 368 // create a string mathematica understands 369 // GDB set print repe 15 # if repeated digits is a bother 370 // set print elements 400 # if line doesn't fit 371 char str[1024]; 372 sk_bzero(str, sizeof(str)); 373 SK_SNPRINTF(str, sizeof(str), "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", 374 A, B, C, D); 375 SkPathOpsDebug::MathematicaIze(str, sizeof(str)); 376#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA 377 SkDebugf("%s\n", str); 378#endif 379#endif 380 if (approximately_zero(A) 381 && approximately_zero_when_compared_to(A, B) 382 && approximately_zero_when_compared_to(A, C) 383 && approximately_zero_when_compared_to(A, D)) { // we're just a quadratic 384 return SkDQuad::RootsReal(B, C, D, s); 385 } 386 if (approximately_zero_when_compared_to(D, A) 387 && approximately_zero_when_compared_to(D, B) 388 && approximately_zero_when_compared_to(D, C)) { // 0 is one root 389 int num = SkDQuad::RootsReal(A, B, C, s); 390 for (int i = 0; i < num; ++i) { 391 if (approximately_zero(s[i])) { 392 return num; 393 } 394 } 395 s[num++] = 0; 396 return num; 397 } 398 if (approximately_zero(A + B + C + D)) { // 1 is one root 399 int num = SkDQuad::RootsReal(A, A + B, -D, s); 400 for (int i = 0; i < num; ++i) { 401 if (AlmostDequalUlps(s[i], 1)) { 402 return num; 403 } 404 } 405 s[num++] = 1; 406 return num; 407 } 408 double a, b, c; 409 { 410 double invA = 1 / A; 411 a = B * invA; 412 b = C * invA; 413 c = D * invA; 414 } 415 double a2 = a * a; 416 double Q = (a2 - b * 3) / 9; 417 double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; 418 double R2 = R * R; 419 double Q3 = Q * Q * Q; 420 double R2MinusQ3 = R2 - Q3; 421 double adiv3 = a / 3; 422 double r; 423 double* roots = s; 424 if (R2MinusQ3 < 0) { // we have 3 real roots 425 double theta = acos(R / sqrt(Q3)); 426 double neg2RootQ = -2 * sqrt(Q); 427 428 r = neg2RootQ * cos(theta / 3) - adiv3; 429 *roots++ = r; 430 431 r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; 432 if (!AlmostDequalUlps(s[0], r)) { 433 *roots++ = r; 434 } 435 r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; 436 if (!AlmostDequalUlps(s[0], r) && (roots - s == 1 || !AlmostDequalUlps(s[1], r))) { 437 *roots++ = r; 438 } 439 } else { // we have 1 real root 440 double sqrtR2MinusQ3 = sqrt(R2MinusQ3); 441 double A = fabs(R) + sqrtR2MinusQ3; 442 A = SkDCubeRoot(A); 443 if (R > 0) { 444 A = -A; 445 } 446 if (A != 0) { 447 A += Q / A; 448 } 449 r = A - adiv3; 450 *roots++ = r; 451 if (AlmostDequalUlps((double) R2, (double) Q3)) { 452 r = -A / 2 - adiv3; 453 if (!AlmostDequalUlps(s[0], r)) { 454 *roots++ = r; 455 } 456 } 457 } 458 return static_cast<int>(roots - s); 459} 460 461// from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf 462// c(t) = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3 463// c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2 464// = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2 465static double derivative_at_t(const double* src, double t) { 466 double one_t = 1 - t; 467 double a = src[0]; 468 double b = src[2]; 469 double c = src[4]; 470 double d = src[6]; 471 return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t); 472} 473 474// OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t? 475SkDVector SkDCubic::dxdyAtT(double t) const { 476 SkDVector result = { derivative_at_t(&fPts[0].fX, t), derivative_at_t(&fPts[0].fY, t) }; 477 return result; 478} 479 480// OPTIMIZE? share code with formulate_F1DotF2 481int SkDCubic::findInflections(double tValues[]) const { 482 double Ax = fPts[1].fX - fPts[0].fX; 483 double Ay = fPts[1].fY - fPts[0].fY; 484 double Bx = fPts[2].fX - 2 * fPts[1].fX + fPts[0].fX; 485 double By = fPts[2].fY - 2 * fPts[1].fY + fPts[0].fY; 486 double Cx = fPts[3].fX + 3 * (fPts[1].fX - fPts[2].fX) - fPts[0].fX; 487 double Cy = fPts[3].fY + 3 * (fPts[1].fY - fPts[2].fY) - fPts[0].fY; 488 return SkDQuad::RootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues); 489} 490 491static void formulate_F1DotF2(const double src[], double coeff[4]) { 492 double a = src[2] - src[0]; 493 double b = src[4] - 2 * src[2] + src[0]; 494 double c = src[6] + 3 * (src[2] - src[4]) - src[0]; 495 coeff[0] = c * c; 496 coeff[1] = 3 * b * c; 497 coeff[2] = 2 * b * b + c * a; 498 coeff[3] = a * b; 499} 500 501/** SkDCubic'(t) = At^2 + Bt + C, where 502 A = 3(-a + 3(b - c) + d) 503 B = 6(a - 2b + c) 504 C = 3(b - a) 505 Solve for t, keeping only those that fit between 0 < t < 1 506*/ 507int SkDCubic::FindExtrema(const double src[], double tValues[2]) { 508 // we divide A,B,C by 3 to simplify 509 double a = src[0]; 510 double b = src[2]; 511 double c = src[4]; 512 double d = src[6]; 513 double A = d - a + 3 * (b - c); 514 double B = 2 * (a - b - b + c); 515 double C = b - a; 516 517 return SkDQuad::RootsValidT(A, B, C, tValues); 518} 519 520/* from SkGeometry.cpp 521 Looking for F' dot F'' == 0 522 523 A = b - a 524 B = c - 2b + a 525 C = d - 3c + 3b - a 526 527 F' = 3Ct^2 + 6Bt + 3A 528 F'' = 6Ct + 6B 529 530 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 531*/ 532int SkDCubic::findMaxCurvature(double tValues[]) const { 533 double coeffX[4], coeffY[4]; 534 int i; 535 formulate_F1DotF2(&fPts[0].fX, coeffX); 536 formulate_F1DotF2(&fPts[0].fY, coeffY); 537 for (i = 0; i < 4; i++) { 538 coeffX[i] = coeffX[i] + coeffY[i]; 539 } 540 return RootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues); 541} 542 543SkDPoint SkDCubic::ptAtT(double t) const { 544 if (0 == t) { 545 return fPts[0]; 546 } 547 if (1 == t) { 548 return fPts[3]; 549 } 550 double one_t = 1 - t; 551 double one_t2 = one_t * one_t; 552 double a = one_t2 * one_t; 553 double b = 3 * one_t2 * t; 554 double t2 = t * t; 555 double c = 3 * one_t * t2; 556 double d = t2 * t; 557 SkDPoint result = {a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX + d * fPts[3].fX, 558 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY + d * fPts[3].fY}; 559 return result; 560} 561 562/* 563 Given a cubic c, t1, and t2, find a small cubic segment. 564 565 The new cubic is defined as points A, B, C, and D, where 566 s1 = 1 - t1 567 s2 = 1 - t2 568 A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1 569 D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2 570 571 We don't have B or C. So We define two equations to isolate them. 572 First, compute two reference T values 1/3 and 2/3 from t1 to t2: 573 574 c(at (2*t1 + t2)/3) == E 575 c(at (t1 + 2*t2)/3) == F 576 577 Next, compute where those values must be if we know the values of B and C: 578 579 _12 = A*2/3 + B*1/3 580 12_ = A*1/3 + B*2/3 581 _23 = B*2/3 + C*1/3 582 23_ = B*1/3 + C*2/3 583 _34 = C*2/3 + D*1/3 584 34_ = C*1/3 + D*2/3 585 _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 586 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 587 _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 588 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 589 _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3 590 = A*8/27 + B*12/27 + C*6/27 + D*1/27 591 = E 592 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3 593 = A*1/27 + B*6/27 + C*12/27 + D*8/27 594 = F 595 E*27 = A*8 + B*12 + C*6 + D 596 F*27 = A + B*6 + C*12 + D*8 597 598Group the known values on one side: 599 600 M = E*27 - A*8 - D = B*12 + C* 6 601 N = F*27 - A - D*8 = B* 6 + C*12 602 M*2 - N = B*18 603 N*2 - M = C*18 604 B = (M*2 - N)/18 605 C = (N*2 - M)/18 606 */ 607 608static double interp_cubic_coords(const double* src, double t) { 609 double ab = SkDInterp(src[0], src[2], t); 610 double bc = SkDInterp(src[2], src[4], t); 611 double cd = SkDInterp(src[4], src[6], t); 612 double abc = SkDInterp(ab, bc, t); 613 double bcd = SkDInterp(bc, cd, t); 614 double abcd = SkDInterp(abc, bcd, t); 615 return abcd; 616} 617 618SkDCubic SkDCubic::subDivide(double t1, double t2) const { 619 if (t1 == 0 || t2 == 1) { 620 if (t1 == 0 && t2 == 1) { 621 return *this; 622 } 623 SkDCubicPair pair = chopAt(t1 == 0 ? t2 : t1); 624 SkDCubic dst = t1 == 0 ? pair.first() : pair.second(); 625 return dst; 626 } 627 SkDCubic dst; 628 double ax = dst[0].fX = interp_cubic_coords(&fPts[0].fX, t1); 629 double ay = dst[0].fY = interp_cubic_coords(&fPts[0].fY, t1); 630 double ex = interp_cubic_coords(&fPts[0].fX, (t1*2+t2)/3); 631 double ey = interp_cubic_coords(&fPts[0].fY, (t1*2+t2)/3); 632 double fx = interp_cubic_coords(&fPts[0].fX, (t1+t2*2)/3); 633 double fy = interp_cubic_coords(&fPts[0].fY, (t1+t2*2)/3); 634 double dx = dst[3].fX = interp_cubic_coords(&fPts[0].fX, t2); 635 double dy = dst[3].fY = interp_cubic_coords(&fPts[0].fY, t2); 636 double mx = ex * 27 - ax * 8 - dx; 637 double my = ey * 27 - ay * 8 - dy; 638 double nx = fx * 27 - ax - dx * 8; 639 double ny = fy * 27 - ay - dy * 8; 640 /* bx = */ dst[1].fX = (mx * 2 - nx) / 18; 641 /* by = */ dst[1].fY = (my * 2 - ny) / 18; 642 /* cx = */ dst[2].fX = (nx * 2 - mx) / 18; 643 /* cy = */ dst[2].fY = (ny * 2 - my) / 18; 644 // FIXME: call align() ? 645 return dst; 646} 647 648void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d, 649 double t1, double t2, SkDPoint dst[2]) const { 650 SkASSERT(t1 != t2); 651 // this approach assumes that the control points computed directly are accurate enough 652 SkDCubic sub = subDivide(t1, t2); 653 dst[0] = sub[1] + (a - sub[0]); 654 dst[1] = sub[2] + (d - sub[3]); 655 if (t1 == 0 || t2 == 0) { 656 align(0, 1, t1 == 0 ? &dst[0] : &dst[1]); 657 } 658 if (t1 == 1 || t2 == 1) { 659 align(3, 2, t1 == 1 ? &dst[0] : &dst[1]); 660 } 661 if (AlmostBequalUlps(dst[0].fX, a.fX)) { 662 dst[0].fX = a.fX; 663 } 664 if (AlmostBequalUlps(dst[0].fY, a.fY)) { 665 dst[0].fY = a.fY; 666 } 667 if (AlmostBequalUlps(dst[1].fX, d.fX)) { 668 dst[1].fX = d.fX; 669 } 670 if (AlmostBequalUlps(dst[1].fY, d.fY)) { 671 dst[1].fY = d.fY; 672 } 673} 674 675double SkDCubic::top(const SkDCubic& dCurve, double startT, double endT, SkDPoint*topPt) const { 676 double extremeTs[2]; 677 double topT = -1; 678 int roots = SkDCubic::FindExtrema(&fPts[0].fY, extremeTs); 679 for (int index = 0; index < roots; ++index) { 680 double t = startT + (endT - startT) * extremeTs[index]; 681 SkDPoint mid = dCurve.ptAtT(t); 682 if (topPt->fY > mid.fY || (topPt->fY == mid.fY && topPt->fX > mid.fX)) { 683 topT = t; 684 *topPt = mid; 685 } 686 } 687 return topT; 688} 689