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