1 /* 2 * Copyright 2011 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 8 #include "GrPathUtils.h" 9 10 #include "GrTypes.h" 11 #include "SkGeometry.h" 12 #include "SkMathPriv.h" 13 14 SkScalar GrPathUtils::scaleToleranceToSrc(SkScalar devTol, 15 const SkMatrix& viewM, 16 const SkRect& pathBounds) { 17 // In order to tesselate the path we get a bound on how much the matrix can 18 // scale when mapping to screen coordinates. 19 SkScalar stretch = viewM.getMaxScale(); 20 SkScalar srcTol = devTol; 21 22 if (stretch < 0) { 23 // take worst case mapRadius amoung four corners. 24 // (less than perfect) 25 for (int i = 0; i < 4; ++i) { 26 SkMatrix mat; 27 mat.setTranslate((i % 2) ? pathBounds.fLeft : pathBounds.fRight, 28 (i < 2) ? pathBounds.fTop : pathBounds.fBottom); 29 mat.postConcat(viewM); 30 stretch = SkMaxScalar(stretch, mat.mapRadius(SK_Scalar1)); 31 } 32 } 33 return srcTol / stretch; 34 } 35 36 static const int MAX_POINTS_PER_CURVE = 1 << 10; 37 static const SkScalar gMinCurveTol = 0.0001f; 38 39 uint32_t GrPathUtils::quadraticPointCount(const SkPoint points[], 40 SkScalar tol) { 41 if (tol < gMinCurveTol) { 42 tol = gMinCurveTol; 43 } 44 SkASSERT(tol > 0); 45 46 SkScalar d = points[1].distanceToLineSegmentBetween(points[0], points[2]); 47 if (!SkScalarIsFinite(d)) { 48 return MAX_POINTS_PER_CURVE; 49 } else if (d <= tol) { 50 return 1; 51 } else { 52 // Each time we subdivide, d should be cut in 4. So we need to 53 // subdivide x = log4(d/tol) times. x subdivisions creates 2^(x) 54 // points. 55 // 2^(log4(x)) = sqrt(x); 56 SkScalar divSqrt = SkScalarSqrt(d / tol); 57 if (((SkScalar)SK_MaxS32) <= divSqrt) { 58 return MAX_POINTS_PER_CURVE; 59 } else { 60 int temp = SkScalarCeilToInt(divSqrt); 61 int pow2 = GrNextPow2(temp); 62 // Because of NaNs & INFs we can wind up with a degenerate temp 63 // such that pow2 comes out negative. Also, our point generator 64 // will always output at least one pt. 65 if (pow2 < 1) { 66 pow2 = 1; 67 } 68 return SkTMin(pow2, MAX_POINTS_PER_CURVE); 69 } 70 } 71 } 72 73 uint32_t GrPathUtils::generateQuadraticPoints(const SkPoint& p0, 74 const SkPoint& p1, 75 const SkPoint& p2, 76 SkScalar tolSqd, 77 SkPoint** points, 78 uint32_t pointsLeft) { 79 if (pointsLeft < 2 || 80 (p1.distanceToLineSegmentBetweenSqd(p0, p2)) < tolSqd) { 81 (*points)[0] = p2; 82 *points += 1; 83 return 1; 84 } 85 86 SkPoint q[] = { 87 { SkScalarAve(p0.fX, p1.fX), SkScalarAve(p0.fY, p1.fY) }, 88 { SkScalarAve(p1.fX, p2.fX), SkScalarAve(p1.fY, p2.fY) }, 89 }; 90 SkPoint r = { SkScalarAve(q[0].fX, q[1].fX), SkScalarAve(q[0].fY, q[1].fY) }; 91 92 pointsLeft >>= 1; 93 uint32_t a = generateQuadraticPoints(p0, q[0], r, tolSqd, points, pointsLeft); 94 uint32_t b = generateQuadraticPoints(r, q[1], p2, tolSqd, points, pointsLeft); 95 return a + b; 96 } 97 98 uint32_t GrPathUtils::cubicPointCount(const SkPoint points[], 99 SkScalar tol) { 100 if (tol < gMinCurveTol) { 101 tol = gMinCurveTol; 102 } 103 SkASSERT(tol > 0); 104 105 SkScalar d = SkTMax( 106 points[1].distanceToLineSegmentBetweenSqd(points[0], points[3]), 107 points[2].distanceToLineSegmentBetweenSqd(points[0], points[3])); 108 d = SkScalarSqrt(d); 109 if (!SkScalarIsFinite(d)) { 110 return MAX_POINTS_PER_CURVE; 111 } else if (d <= tol) { 112 return 1; 113 } else { 114 SkScalar divSqrt = SkScalarSqrt(d / tol); 115 if (((SkScalar)SK_MaxS32) <= divSqrt) { 116 return MAX_POINTS_PER_CURVE; 117 } else { 118 int temp = SkScalarCeilToInt(SkScalarSqrt(d / tol)); 119 int pow2 = GrNextPow2(temp); 120 // Because of NaNs & INFs we can wind up with a degenerate temp 121 // such that pow2 comes out negative. Also, our point generator 122 // will always output at least one pt. 123 if (pow2 < 1) { 124 pow2 = 1; 125 } 126 return SkTMin(pow2, MAX_POINTS_PER_CURVE); 127 } 128 } 129 } 130 131 uint32_t GrPathUtils::generateCubicPoints(const SkPoint& p0, 132 const SkPoint& p1, 133 const SkPoint& p2, 134 const SkPoint& p3, 135 SkScalar tolSqd, 136 SkPoint** points, 137 uint32_t pointsLeft) { 138 if (pointsLeft < 2 || 139 (p1.distanceToLineSegmentBetweenSqd(p0, p3) < tolSqd && 140 p2.distanceToLineSegmentBetweenSqd(p0, p3) < tolSqd)) { 141 (*points)[0] = p3; 142 *points += 1; 143 return 1; 144 } 145 SkPoint q[] = { 146 { SkScalarAve(p0.fX, p1.fX), SkScalarAve(p0.fY, p1.fY) }, 147 { SkScalarAve(p1.fX, p2.fX), SkScalarAve(p1.fY, p2.fY) }, 148 { SkScalarAve(p2.fX, p3.fX), SkScalarAve(p2.fY, p3.fY) } 149 }; 150 SkPoint r[] = { 151 { SkScalarAve(q[0].fX, q[1].fX), SkScalarAve(q[0].fY, q[1].fY) }, 152 { SkScalarAve(q[1].fX, q[2].fX), SkScalarAve(q[1].fY, q[2].fY) } 153 }; 154 SkPoint s = { SkScalarAve(r[0].fX, r[1].fX), SkScalarAve(r[0].fY, r[1].fY) }; 155 pointsLeft >>= 1; 156 uint32_t a = generateCubicPoints(p0, q[0], r[0], s, tolSqd, points, pointsLeft); 157 uint32_t b = generateCubicPoints(s, r[1], q[2], p3, tolSqd, points, pointsLeft); 158 return a + b; 159 } 160 161 int GrPathUtils::worstCasePointCount(const SkPath& path, int* subpaths, 162 SkScalar tol) { 163 if (tol < gMinCurveTol) { 164 tol = gMinCurveTol; 165 } 166 SkASSERT(tol > 0); 167 168 int pointCount = 0; 169 *subpaths = 1; 170 171 bool first = true; 172 173 SkPath::Iter iter(path, false); 174 SkPath::Verb verb; 175 176 SkPoint pts[4]; 177 while ((verb = iter.next(pts)) != SkPath::kDone_Verb) { 178 179 switch (verb) { 180 case SkPath::kLine_Verb: 181 pointCount += 1; 182 break; 183 case SkPath::kConic_Verb: { 184 SkScalar weight = iter.conicWeight(); 185 SkAutoConicToQuads converter; 186 const SkPoint* quadPts = converter.computeQuads(pts, weight, 0.25f); 187 for (int i = 0; i < converter.countQuads(); ++i) { 188 pointCount += quadraticPointCount(quadPts + 2*i, tol); 189 } 190 } 191 case SkPath::kQuad_Verb: 192 pointCount += quadraticPointCount(pts, tol); 193 break; 194 case SkPath::kCubic_Verb: 195 pointCount += cubicPointCount(pts, tol); 196 break; 197 case SkPath::kMove_Verb: 198 pointCount += 1; 199 if (!first) { 200 ++(*subpaths); 201 } 202 break; 203 default: 204 break; 205 } 206 first = false; 207 } 208 return pointCount; 209 } 210 211 void GrPathUtils::QuadUVMatrix::set(const SkPoint qPts[3]) { 212 SkMatrix m; 213 // We want M such that M * xy_pt = uv_pt 214 // We know M * control_pts = [0 1/2 1] 215 // [0 0 1] 216 // [1 1 1] 217 // And control_pts = [x0 x1 x2] 218 // [y0 y1 y2] 219 // [1 1 1 ] 220 // We invert the control pt matrix and post concat to both sides to get M. 221 // Using the known form of the control point matrix and the result, we can 222 // optimize and improve precision. 223 224 double x0 = qPts[0].fX; 225 double y0 = qPts[0].fY; 226 double x1 = qPts[1].fX; 227 double y1 = qPts[1].fY; 228 double x2 = qPts[2].fX; 229 double y2 = qPts[2].fY; 230 double det = x0*y1 - y0*x1 + x2*y0 - y2*x0 + x1*y2 - y1*x2; 231 232 if (!sk_float_isfinite(det) 233 || SkScalarNearlyZero((float)det, SK_ScalarNearlyZero * SK_ScalarNearlyZero)) { 234 // The quad is degenerate. Hopefully this is rare. Find the pts that are 235 // farthest apart to compute a line (unless it is really a pt). 236 SkScalar maxD = qPts[0].distanceToSqd(qPts[1]); 237 int maxEdge = 0; 238 SkScalar d = qPts[1].distanceToSqd(qPts[2]); 239 if (d > maxD) { 240 maxD = d; 241 maxEdge = 1; 242 } 243 d = qPts[2].distanceToSqd(qPts[0]); 244 if (d > maxD) { 245 maxD = d; 246 maxEdge = 2; 247 } 248 // We could have a tolerance here, not sure if it would improve anything 249 if (maxD > 0) { 250 // Set the matrix to give (u = 0, v = distance_to_line) 251 SkVector lineVec = qPts[(maxEdge + 1)%3] - qPts[maxEdge]; 252 // when looking from the point 0 down the line we want positive 253 // distances to be to the left. This matches the non-degenerate 254 // case. 255 lineVec.setOrthog(lineVec, SkPoint::kLeft_Side); 256 // first row 257 fM[0] = 0; 258 fM[1] = 0; 259 fM[2] = 0; 260 // second row 261 fM[3] = lineVec.fX; 262 fM[4] = lineVec.fY; 263 fM[5] = -lineVec.dot(qPts[maxEdge]); 264 } else { 265 // It's a point. It should cover zero area. Just set the matrix such 266 // that (u, v) will always be far away from the quad. 267 fM[0] = 0; fM[1] = 0; fM[2] = 100.f; 268 fM[3] = 0; fM[4] = 0; fM[5] = 100.f; 269 } 270 } else { 271 double scale = 1.0/det; 272 273 // compute adjugate matrix 274 double a2, a3, a4, a5, a6, a7, a8; 275 a2 = x1*y2-x2*y1; 276 277 a3 = y2-y0; 278 a4 = x0-x2; 279 a5 = x2*y0-x0*y2; 280 281 a6 = y0-y1; 282 a7 = x1-x0; 283 a8 = x0*y1-x1*y0; 284 285 // this performs the uv_pts*adjugate(control_pts) multiply, 286 // then does the scale by 1/det afterwards to improve precision 287 m[SkMatrix::kMScaleX] = (float)((0.5*a3 + a6)*scale); 288 m[SkMatrix::kMSkewX] = (float)((0.5*a4 + a7)*scale); 289 m[SkMatrix::kMTransX] = (float)((0.5*a5 + a8)*scale); 290 291 m[SkMatrix::kMSkewY] = (float)(a6*scale); 292 m[SkMatrix::kMScaleY] = (float)(a7*scale); 293 m[SkMatrix::kMTransY] = (float)(a8*scale); 294 295 // kMPersp0 & kMPersp1 should algebraically be zero 296 m[SkMatrix::kMPersp0] = 0.0f; 297 m[SkMatrix::kMPersp1] = 0.0f; 298 m[SkMatrix::kMPersp2] = (float)((a2 + a5 + a8)*scale); 299 300 // It may not be normalized to have 1.0 in the bottom right 301 float m33 = m.get(SkMatrix::kMPersp2); 302 if (1.f != m33) { 303 m33 = 1.f / m33; 304 fM[0] = m33 * m.get(SkMatrix::kMScaleX); 305 fM[1] = m33 * m.get(SkMatrix::kMSkewX); 306 fM[2] = m33 * m.get(SkMatrix::kMTransX); 307 fM[3] = m33 * m.get(SkMatrix::kMSkewY); 308 fM[4] = m33 * m.get(SkMatrix::kMScaleY); 309 fM[5] = m33 * m.get(SkMatrix::kMTransY); 310 } else { 311 fM[0] = m.get(SkMatrix::kMScaleX); 312 fM[1] = m.get(SkMatrix::kMSkewX); 313 fM[2] = m.get(SkMatrix::kMTransX); 314 fM[3] = m.get(SkMatrix::kMSkewY); 315 fM[4] = m.get(SkMatrix::kMScaleY); 316 fM[5] = m.get(SkMatrix::kMTransY); 317 } 318 } 319 } 320 321 //////////////////////////////////////////////////////////////////////////////// 322 323 // k = (y2 - y0, x0 - x2, x2*y0 - x0*y2) 324 // l = (y1 - y0, x0 - x1, x1*y0 - x0*y1) * 2*w 325 // m = (y2 - y1, x1 - x2, x2*y1 - x1*y2) * 2*w 326 void GrPathUtils::getConicKLM(const SkPoint p[3], const SkScalar weight, SkMatrix* out) { 327 SkMatrix& klm = *out; 328 const SkScalar w2 = 2.f * weight; 329 klm[0] = p[2].fY - p[0].fY; 330 klm[1] = p[0].fX - p[2].fX; 331 klm[2] = p[2].fX * p[0].fY - p[0].fX * p[2].fY; 332 333 klm[3] = w2 * (p[1].fY - p[0].fY); 334 klm[4] = w2 * (p[0].fX - p[1].fX); 335 klm[5] = w2 * (p[1].fX * p[0].fY - p[0].fX * p[1].fY); 336 337 klm[6] = w2 * (p[2].fY - p[1].fY); 338 klm[7] = w2 * (p[1].fX - p[2].fX); 339 klm[8] = w2 * (p[2].fX * p[1].fY - p[1].fX * p[2].fY); 340 341 // scale the max absolute value of coeffs to 10 342 SkScalar scale = 0.f; 343 for (int i = 0; i < 9; ++i) { 344 scale = SkMaxScalar(scale, SkScalarAbs(klm[i])); 345 } 346 SkASSERT(scale > 0.f); 347 scale = 10.f / scale; 348 for (int i = 0; i < 9; ++i) { 349 klm[i] *= scale; 350 } 351 } 352 353 //////////////////////////////////////////////////////////////////////////////// 354 355 namespace { 356 357 // a is the first control point of the cubic. 358 // ab is the vector from a to the second control point. 359 // dc is the vector from the fourth to the third control point. 360 // d is the fourth control point. 361 // p is the candidate quadratic control point. 362 // this assumes that the cubic doesn't inflect and is simple 363 bool is_point_within_cubic_tangents(const SkPoint& a, 364 const SkVector& ab, 365 const SkVector& dc, 366 const SkPoint& d, 367 SkPathPriv::FirstDirection dir, 368 const SkPoint p) { 369 SkVector ap = p - a; 370 SkScalar apXab = ap.cross(ab); 371 if (SkPathPriv::kCW_FirstDirection == dir) { 372 if (apXab > 0) { 373 return false; 374 } 375 } else { 376 SkASSERT(SkPathPriv::kCCW_FirstDirection == dir); 377 if (apXab < 0) { 378 return false; 379 } 380 } 381 382 SkVector dp = p - d; 383 SkScalar dpXdc = dp.cross(dc); 384 if (SkPathPriv::kCW_FirstDirection == dir) { 385 if (dpXdc < 0) { 386 return false; 387 } 388 } else { 389 SkASSERT(SkPathPriv::kCCW_FirstDirection == dir); 390 if (dpXdc > 0) { 391 return false; 392 } 393 } 394 return true; 395 } 396 397 void convert_noninflect_cubic_to_quads(const SkPoint p[4], 398 SkScalar toleranceSqd, 399 bool constrainWithinTangents, 400 SkPathPriv::FirstDirection dir, 401 SkTArray<SkPoint, true>* quads, 402 int sublevel = 0) { 403 404 // Notation: Point a is always p[0]. Point b is p[1] unless p[1] == p[0], in which case it is 405 // p[2]. Point d is always p[3]. Point c is p[2] unless p[2] == p[3], in which case it is p[1]. 406 407 SkVector ab = p[1] - p[0]; 408 SkVector dc = p[2] - p[3]; 409 410 if (ab.lengthSqd() < SK_ScalarNearlyZero) { 411 if (dc.lengthSqd() < SK_ScalarNearlyZero) { 412 SkPoint* degQuad = quads->push_back_n(3); 413 degQuad[0] = p[0]; 414 degQuad[1] = p[0]; 415 degQuad[2] = p[3]; 416 return; 417 } 418 ab = p[2] - p[0]; 419 } 420 if (dc.lengthSqd() < SK_ScalarNearlyZero) { 421 dc = p[1] - p[3]; 422 } 423 424 // When the ab and cd tangents are degenerate or nearly parallel with vector from d to a the 425 // constraint that the quad point falls between the tangents becomes hard to enforce and we are 426 // likely to hit the max subdivision count. However, in this case the cubic is approaching a 427 // line and the accuracy of the quad point isn't so important. We check if the two middle cubic 428 // control points are very close to the baseline vector. If so then we just pick quadratic 429 // points on the control polygon. 430 431 if (constrainWithinTangents) { 432 SkVector da = p[0] - p[3]; 433 bool doQuads = dc.lengthSqd() < SK_ScalarNearlyZero || 434 ab.lengthSqd() < SK_ScalarNearlyZero; 435 if (!doQuads) { 436 SkScalar invDALengthSqd = da.lengthSqd(); 437 if (invDALengthSqd > SK_ScalarNearlyZero) { 438 invDALengthSqd = SkScalarInvert(invDALengthSqd); 439 // cross(ab, da)^2/length(da)^2 == sqd distance from b to line from d to a. 440 // same goes for point c using vector cd. 441 SkScalar detABSqd = ab.cross(da); 442 detABSqd = SkScalarSquare(detABSqd); 443 SkScalar detDCSqd = dc.cross(da); 444 detDCSqd = SkScalarSquare(detDCSqd); 445 if (detABSqd * invDALengthSqd < toleranceSqd && 446 detDCSqd * invDALengthSqd < toleranceSqd) 447 { 448 doQuads = true; 449 } 450 } 451 } 452 if (doQuads) { 453 SkPoint b = p[0] + ab; 454 SkPoint c = p[3] + dc; 455 SkPoint mid = b + c; 456 mid.scale(SK_ScalarHalf); 457 // Insert two quadratics to cover the case when ab points away from d and/or dc 458 // points away from a. 459 if (SkVector::DotProduct(da, dc) < 0 || SkVector::DotProduct(ab,da) > 0) { 460 SkPoint* qpts = quads->push_back_n(6); 461 qpts[0] = p[0]; 462 qpts[1] = b; 463 qpts[2] = mid; 464 qpts[3] = mid; 465 qpts[4] = c; 466 qpts[5] = p[3]; 467 } else { 468 SkPoint* qpts = quads->push_back_n(3); 469 qpts[0] = p[0]; 470 qpts[1] = mid; 471 qpts[2] = p[3]; 472 } 473 return; 474 } 475 } 476 477 static const SkScalar kLengthScale = 3 * SK_Scalar1 / 2; 478 static const int kMaxSubdivs = 10; 479 480 ab.scale(kLengthScale); 481 dc.scale(kLengthScale); 482 483 // e0 and e1 are extrapolations along vectors ab and dc. 484 SkVector c0 = p[0]; 485 c0 += ab; 486 SkVector c1 = p[3]; 487 c1 += dc; 488 489 SkScalar dSqd = sublevel > kMaxSubdivs ? 0 : c0.distanceToSqd(c1); 490 if (dSqd < toleranceSqd) { 491 SkPoint cAvg = c0; 492 cAvg += c1; 493 cAvg.scale(SK_ScalarHalf); 494 495 bool subdivide = false; 496 497 if (constrainWithinTangents && 498 !is_point_within_cubic_tangents(p[0], ab, dc, p[3], dir, cAvg)) { 499 // choose a new cAvg that is the intersection of the two tangent lines. 500 ab.setOrthog(ab); 501 SkScalar z0 = -ab.dot(p[0]); 502 dc.setOrthog(dc); 503 SkScalar z1 = -dc.dot(p[3]); 504 cAvg.fX = ab.fY * z1 - z0 * dc.fY; 505 cAvg.fY = z0 * dc.fX - ab.fX * z1; 506 SkScalar z = ab.fX * dc.fY - ab.fY * dc.fX; 507 z = SkScalarInvert(z); 508 cAvg.fX *= z; 509 cAvg.fY *= z; 510 if (sublevel <= kMaxSubdivs) { 511 SkScalar d0Sqd = c0.distanceToSqd(cAvg); 512 SkScalar d1Sqd = c1.distanceToSqd(cAvg); 513 // We need to subdivide if d0 + d1 > tolerance but we have the sqd values. We know 514 // the distances and tolerance can't be negative. 515 // (d0 + d1)^2 > toleranceSqd 516 // d0Sqd + 2*d0*d1 + d1Sqd > toleranceSqd 517 SkScalar d0d1 = SkScalarSqrt(d0Sqd * d1Sqd); 518 subdivide = 2 * d0d1 + d0Sqd + d1Sqd > toleranceSqd; 519 } 520 } 521 if (!subdivide) { 522 SkPoint* pts = quads->push_back_n(3); 523 pts[0] = p[0]; 524 pts[1] = cAvg; 525 pts[2] = p[3]; 526 return; 527 } 528 } 529 SkPoint choppedPts[7]; 530 SkChopCubicAtHalf(p, choppedPts); 531 convert_noninflect_cubic_to_quads(choppedPts + 0, 532 toleranceSqd, 533 constrainWithinTangents, 534 dir, 535 quads, 536 sublevel + 1); 537 convert_noninflect_cubic_to_quads(choppedPts + 3, 538 toleranceSqd, 539 constrainWithinTangents, 540 dir, 541 quads, 542 sublevel + 1); 543 } 544 } 545 546 void GrPathUtils::convertCubicToQuads(const SkPoint p[4], 547 SkScalar tolScale, 548 SkTArray<SkPoint, true>* quads) { 549 SkPoint chopped[10]; 550 int count = SkChopCubicAtInflections(p, chopped); 551 552 const SkScalar tolSqd = SkScalarSquare(tolScale); 553 554 for (int i = 0; i < count; ++i) { 555 SkPoint* cubic = chopped + 3*i; 556 // The direction param is ignored if the third param is false. 557 convert_noninflect_cubic_to_quads(cubic, tolSqd, false, 558 SkPathPriv::kCCW_FirstDirection, quads); 559 } 560 } 561 562 void GrPathUtils::convertCubicToQuadsConstrainToTangents(const SkPoint p[4], 563 SkScalar tolScale, 564 SkPathPriv::FirstDirection dir, 565 SkTArray<SkPoint, true>* quads) { 566 SkPoint chopped[10]; 567 int count = SkChopCubicAtInflections(p, chopped); 568 569 const SkScalar tolSqd = SkScalarSquare(tolScale); 570 571 for (int i = 0; i < count; ++i) { 572 SkPoint* cubic = chopped + 3*i; 573 convert_noninflect_cubic_to_quads(cubic, tolSqd, true, dir, quads); 574 } 575 } 576 577 //////////////////////////////////////////////////////////////////////////////// 578 579 /** 580 * Computes an SkMatrix that can find the cubic KLM functionals as follows: 581 * 582 * | ..K.. | | ..kcoeffs.. | 583 * | ..L.. | = | ..lcoeffs.. | * inverse_transpose_power_basis_matrix 584 * | ..M.. | | ..mcoeffs.. | 585 * 586 * 'kcoeffs' are the power basis coefficients to a scalar valued cubic function that returns the 587 * signed distance to line K from a given point on the curve: 588 * 589 * k(t,s) = C(t,s) * K [C(t,s) is defined in the following comment] 590 * 591 * The same applies for lcoeffs and mcoeffs. These are found separately, depending on the type of 592 * curve. There are 4 coefficients but 3 rows in the matrix, so in order to do this calculation the 593 * caller must first remove a specific column of coefficients. 594 * 595 * @return which column of klm coefficients to exclude from the calculation. 596 */ 597 static int calc_inverse_transpose_power_basis_matrix(const SkPoint pts[4], SkMatrix* out) { 598 using SkScalar4 = SkNx<4, SkScalar>; 599 600 // First we convert the bezier coordinates 'pts' to power basis coefficients X,Y,W=[0 0 0 1]. 601 // M3 is the matrix that does this conversion. The homogeneous equation for the cubic becomes: 602 // 603 // | X Y 0 | 604 // C(t,s) = [t^3 t^2*s t*s^2 s^3] * | . . 0 | 605 // | . . 0 | 606 // | . . 1 | 607 // 608 const SkScalar4 M3[3] = {SkScalar4(-1, 3, -3, 1), 609 SkScalar4(3, -6, 3, 0), 610 SkScalar4(-3, 3, 0, 0)}; 611 // 4th column of M3 = SkScalar4(1, 0, 0, 0)}; 612 SkScalar4 X(pts[3].x(), 0, 0, 0); 613 SkScalar4 Y(pts[3].y(), 0, 0, 0); 614 for (int i = 2; i >= 0; --i) { 615 X += M3[i] * pts[i].x(); 616 Y += M3[i] * pts[i].y(); 617 } 618 619 // The matrix is 3x4. In order to invert it, we first need to make it square by throwing out one 620 // of the top three rows. We toss the row that leaves us with the largest determinant. Since the 621 // right column will be [0 0 1], the determinant reduces to x0*y1 - y0*x1. 622 SkScalar det[4]; 623 SkScalar4 DETX1 = SkNx_shuffle<1,0,0,3>(X), DETY1 = SkNx_shuffle<1,0,0,3>(Y); 624 SkScalar4 DETX2 = SkNx_shuffle<2,2,1,3>(X), DETY2 = SkNx_shuffle<2,2,1,3>(Y); 625 (DETX1 * DETY2 - DETY1 * DETX2).store(det); 626 const int skipRow = det[0] > det[2] ? (det[0] > det[1] ? 0 : 1) 627 : (det[1] > det[2] ? 1 : 2); 628 const SkScalar rdet = 1 / det[skipRow]; 629 const int row0 = (0 != skipRow) ? 0 : 1; 630 const int row1 = (2 == skipRow) ? 1 : 2; 631 632 // Compute the inverse-transpose of the power basis matrix with the 'skipRow'th row removed. 633 // Since W=[0 0 0 1], it follows that our corresponding solution will be equal to: 634 // 635 // | y1 -x1 x1*y2 - y1*x2 | 636 // 1/det * | -y0 x0 -x0*y2 + y0*x2 | 637 // | 0 0 det | 638 // 639 const SkScalar4 R(rdet, rdet, rdet, 1); 640 X *= R; 641 Y *= R; 642 643 SkScalar x[4], y[4], z[4]; 644 X.store(x); 645 Y.store(y); 646 (X * SkNx_shuffle<3,3,3,3>(Y) - Y * SkNx_shuffle<3,3,3,3>(X)).store(z); 647 648 out->setAll( y[row1], -x[row1], z[row1], 649 -y[row0], x[row0], -z[row0], 650 0, 0, 1); 651 652 return skipRow; 653 } 654 655 static void negate_kl(SkMatrix* klm) { 656 // We could use klm->postScale(-1, -1), but it ends up doing a full matrix multiply. 657 for (int i = 0; i < 6; ++i) { 658 (*klm)[i] = -(*klm)[i]; 659 } 660 } 661 662 static void calc_serp_klm(const SkPoint pts[4], const SkScalar d[3], SkMatrix* klm) { 663 SkMatrix CIT; 664 int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT); 665 666 const SkScalar root = SkScalarSqrt(9 * d[1] * d[1] - 12 * d[0] * d[2]); 667 668 const SkScalar tl = 3 * d[1] + root; 669 const SkScalar sl = 6 * d[0]; 670 const SkScalar tm = 3 * d[1] - root; 671 const SkScalar sm = 6 * d[0]; 672 673 SkMatrix klmCoeffs; 674 int col = 0; 675 if (0 != skipCol) { 676 klmCoeffs[0] = 0; 677 klmCoeffs[3] = -sl * sl * sl; 678 klmCoeffs[6] = -sm * sm * sm; 679 ++col; 680 } 681 if (1 != skipCol) { 682 klmCoeffs[col + 0] = sl * sm; 683 klmCoeffs[col + 3] = 3 * sl * sl * tl; 684 klmCoeffs[col + 6] = 3 * sm * sm * tm; 685 ++col; 686 } 687 if (2 != skipCol) { 688 klmCoeffs[col + 0] = -tl * sm - tm * sl; 689 klmCoeffs[col + 3] = -3 * sl * tl * tl; 690 klmCoeffs[col + 6] = -3 * sm * tm * tm; 691 ++col; 692 } 693 694 SkASSERT(2 == col); 695 klmCoeffs[2] = tl * tm; 696 klmCoeffs[5] = tl * tl * tl; 697 klmCoeffs[8] = tm * tm * tm; 698 699 klm->setConcat(klmCoeffs, CIT); 700 701 // If d0 > 0 we need to flip the orientation of our curve 702 // This is done by negating the k and l values 703 // We want negative distance values to be on the inside 704 if (d[0] > 0) { 705 negate_kl(klm); 706 } 707 } 708 709 static void calc_loop_klm(const SkPoint pts[4], SkScalar d1, SkScalar td, SkScalar sd, 710 SkScalar te, SkScalar se, SkMatrix* klm) { 711 SkMatrix CIT; 712 int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT); 713 714 const SkScalar tesd = te * sd; 715 const SkScalar tdse = td * se; 716 717 SkMatrix klmCoeffs; 718 int col = 0; 719 if (0 != skipCol) { 720 klmCoeffs[0] = 0; 721 klmCoeffs[3] = -sd * sd * se; 722 klmCoeffs[6] = -se * se * sd; 723 ++col; 724 } 725 if (1 != skipCol) { 726 klmCoeffs[col + 0] = sd * se; 727 klmCoeffs[col + 3] = sd * (2 * tdse + tesd); 728 klmCoeffs[col + 6] = se * (2 * tesd + tdse); 729 ++col; 730 } 731 if (2 != skipCol) { 732 klmCoeffs[col + 0] = -tdse - tesd; 733 klmCoeffs[col + 3] = -td * (tdse + 2 * tesd); 734 klmCoeffs[col + 6] = -te * (tesd + 2 * tdse); 735 ++col; 736 } 737 738 SkASSERT(2 == col); 739 klmCoeffs[2] = td * te; 740 klmCoeffs[5] = td * td * te; 741 klmCoeffs[8] = te * te * td; 742 743 klm->setConcat(klmCoeffs, CIT); 744 745 // For the general loop curve, we flip the orientation in the same pattern as the serp case 746 // above. Thus we only check d1. Technically we should check the value of the hessian as well 747 // cause we care about the sign of d1*Hessian. However, the Hessian is always negative outside 748 // the loop section and positive inside. We take care of the flipping for the loop sections 749 // later on. 750 if (d1 > 0) { 751 negate_kl(klm); 752 } 753 } 754 755 // For the case when we have a cusp at a parameter value of infinity (discr == 0, d1 == 0). 756 static void calc_inf_cusp_klm(const SkPoint pts[4], SkScalar d2, SkScalar d3, SkMatrix* klm) { 757 SkMatrix CIT; 758 int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT); 759 760 const SkScalar tn = d3; 761 const SkScalar sn = 3 * d2; 762 763 SkMatrix klmCoeffs; 764 int col = 0; 765 if (0 != skipCol) { 766 klmCoeffs[0] = 0; 767 klmCoeffs[3] = -sn * sn * sn; 768 ++col; 769 } 770 if (1 != skipCol) { 771 klmCoeffs[col + 0] = 0; 772 klmCoeffs[col + 3] = 3 * sn * sn * tn; 773 ++col; 774 } 775 if (2 != skipCol) { 776 klmCoeffs[col + 0] = -sn; 777 klmCoeffs[col + 3] = -3 * sn * tn * tn; 778 ++col; 779 } 780 781 SkASSERT(2 == col); 782 klmCoeffs[2] = tn; 783 klmCoeffs[5] = tn * tn * tn; 784 785 klmCoeffs[6] = 0; 786 klmCoeffs[7] = 0; 787 klmCoeffs[8] = 1; 788 789 klm->setConcat(klmCoeffs, CIT); 790 } 791 792 // For the case when a cubic bezier is actually a quadratic. We duplicate k in l so that the 793 // implicit becomes: 794 // 795 // k^3 - l*m == k^3 - l*k == k * (k^2 - l) 796 // 797 // In the quadratic case we can simply assign fixed values at each control point: 798 // 799 // | ..K.. | | pts[0] pts[1] pts[2] pts[3] | | 0 1/3 2/3 1 | 800 // | ..L.. | * | . . . . | == | 0 0 1/3 1 | 801 // | ..K.. | | 1 1 1 1 | | 0 1/3 2/3 1 | 802 // 803 static void calc_quadratic_klm(const SkPoint pts[4], SkScalar d3, SkMatrix* klm) { 804 SkMatrix klmAtPts; 805 klmAtPts.setAll(0, 1.f/3, 1, 806 0, 0, 1, 807 0, 1.f/3, 1); 808 809 SkMatrix inversePts; 810 inversePts.setAll(pts[0].x(), pts[1].x(), pts[3].x(), 811 pts[0].y(), pts[1].y(), pts[3].y(), 812 1, 1, 1); 813 SkAssertResult(inversePts.invert(&inversePts)); 814 815 klm->setConcat(klmAtPts, inversePts); 816 817 // If d3 > 0 we need to flip the orientation of our curve 818 // This is done by negating the k and l values 819 if (d3 > 0) { 820 negate_kl(klm); 821 } 822 } 823 824 // For the case when a cubic bezier is actually a line. We set K=0, L=1, M=-line, which results in 825 // the following implicit: 826 // 827 // k^3 - l*m == 0^3 - 1*(-line) == -(-line) == line 828 // 829 static void calc_line_klm(const SkPoint pts[4], SkMatrix* klm) { 830 SkScalar ny = pts[0].x() - pts[3].x(); 831 SkScalar nx = pts[3].y() - pts[0].y(); 832 SkScalar k = nx * pts[0].x() + ny * pts[0].y(); 833 klm->setAll( 0, 0, 0, 834 0, 0, 1, 835 -nx, -ny, k); 836 } 837 838 int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[10], SkMatrix* klm, 839 int* loopIndex) { 840 // Variables to store the two parametric values at the loop double point. 841 SkScalar t1 = 0, t2 = 0; 842 843 // Homogeneous parametric values at the loop double point. 844 SkScalar td, sd, te, se; 845 846 SkScalar d[3]; 847 SkCubicType cType = SkClassifyCubic(src, d); 848 849 int chop_count = 0; 850 if (kLoop_SkCubicType == cType) { 851 SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]); 852 td = d[1] + tempSqrt; 853 sd = 2.f * d[0]; 854 te = d[1] - tempSqrt; 855 se = 2.f * d[0]; 856 857 t1 = td / sd; 858 t2 = te / se; 859 // need to have t values sorted since this is what is expected by SkChopCubicAt 860 if (t1 > t2) { 861 SkTSwap(t1, t2); 862 } 863 864 SkScalar chop_ts[2]; 865 if (t1 > 0.f && t1 < 1.f) { 866 chop_ts[chop_count++] = t1; 867 } 868 if (t2 > 0.f && t2 < 1.f) { 869 chop_ts[chop_count++] = t2; 870 } 871 if(dst) { 872 SkChopCubicAt(src, dst, chop_ts, chop_count); 873 } 874 } else { 875 if (dst) { 876 memcpy(dst, src, sizeof(SkPoint) * 4); 877 } 878 } 879 880 if (loopIndex) { 881 if (2 == chop_count) { 882 *loopIndex = 1; 883 } else if (1 == chop_count) { 884 if (t1 < 0.f) { 885 *loopIndex = 0; 886 } else { 887 *loopIndex = 1; 888 } 889 } else { 890 if (t1 < 0.f && t2 > 1.f) { 891 *loopIndex = 0; 892 } else { 893 *loopIndex = -1; 894 } 895 } 896 } 897 898 if (klm) { 899 switch (cType) { 900 case kSerpentine_SkCubicType: 901 calc_serp_klm(src, d, klm); 902 break; 903 case kLoop_SkCubicType: 904 calc_loop_klm(src, d[0], td, sd, te, se, klm); 905 break; 906 case kCusp_SkCubicType: 907 if (0 != d[0]) { 908 // FIXME: SkClassifyCubic has a tolerance, but we need an exact classification 909 // here to be sure we won't get a negative in the square root. 910 calc_serp_klm(src, d, klm); 911 } else { 912 calc_inf_cusp_klm(src, d[1], d[2], klm); 913 } 914 break; 915 case kQuadratic_SkCubicType: 916 calc_quadratic_klm(src, d[2], klm); 917 break; 918 case kLine_SkCubicType: 919 case kPoint_SkCubicType: 920 calc_line_klm(src, klm); 921 break; 922 }; 923 } 924 return chop_count + 1; 925 } 926
