Avoid catastrophic cancellation when computing det for QuadUVMatrix

With data with large values we can get into situations where clearly
degenerate quads are not giving a zero determinant, because of
errors when computing the additions and subtractions in sequence.
The solution is to compute the subtractions first, and then add them
together.

Bug: oss-fuzz:41445
Change-Id: I7febead4b020e11a137f1703b761af70c14e4e47
Reviewed-on: https://skia-review.googlesource.com/c/skia/+/511796
Reviewed-by: Michael Ludwig <michaelludwig@google.com>
Commit-Queue: Jim Van Verth <jvanverth@google.com>
This commit is contained in:
Jim Van Verth 2022-02-23 12:47:43 -05:00 committed by SkCQ
parent a8689e3d4d
commit 8b53da2b37

View File

@ -129,7 +129,6 @@ uint32_t GrPathUtils::generateCubicPoints(const SkPoint& p0,
}
void GrPathUtils::QuadUVMatrix::set(const SkPoint qPts[3]) {
SkMatrix m;
// We want M such that M * xy_pt = uv_pt
// We know M * control_pts = [0 1/2 1]
// [0 0 1]
@ -147,7 +146,12 @@ void GrPathUtils::QuadUVMatrix::set(const SkPoint qPts[3]) {
double y1 = qPts[1].fY;
double x2 = qPts[2].fX;
double y2 = qPts[2].fY;
double det = x0*y1 - y0*x1 + x2*y0 - y2*x0 + x1*y2 - y1*x2;
// pre-calculate some adjugate matrix factors for determinant
double a2 = x1*y2-x2*y1;
double a5 = x2*y0-x0*y2;
double a8 = x0*y1-x1*y0;
double det = a2 + a5 + a8;
if (!sk_float_isfinite(det)
|| SkScalarNearlyZero((float)det, SK_ScalarNearlyZero * SK_ScalarNearlyZero)) {
@ -191,50 +195,21 @@ void GrPathUtils::QuadUVMatrix::set(const SkPoint qPts[3]) {
double scale = 1.0/det;
// compute adjugate matrix
double a2, a3, a4, a5, a6, a7, a8;
a2 = x1*y2-x2*y1;
double a3, a4, a6, a7;
a3 = y2-y0;
a4 = x0-x2;
a5 = x2*y0-x0*y2;
a6 = y0-y1;
a7 = x1-x0;
a8 = x0*y1-x1*y0;
// this performs the uv_pts*adjugate(control_pts) multiply,
// then does the scale by 1/det afterwards to improve precision
m[SkMatrix::kMScaleX] = (float)((0.5*a3 + a6)*scale);
m[SkMatrix::kMSkewX] = (float)((0.5*a4 + a7)*scale);
m[SkMatrix::kMTransX] = (float)((0.5*a5 + a8)*scale);
m[SkMatrix::kMSkewY] = (float)(a6*scale);
m[SkMatrix::kMScaleY] = (float)(a7*scale);
m[SkMatrix::kMTransY] = (float)(a8*scale);
// kMPersp0 & kMPersp1 should algebraically be zero
m[SkMatrix::kMPersp0] = 0.0f;
m[SkMatrix::kMPersp1] = 0.0f;
m[SkMatrix::kMPersp2] = (float)((a2 + a5 + a8)*scale);
// It may not be normalized to have 1.0 in the bottom right
float m33 = m.get(SkMatrix::kMPersp2);
if (1.f != m33) {
m33 = 1.f / m33;
fM[0] = m33 * m.get(SkMatrix::kMScaleX);
fM[1] = m33 * m.get(SkMatrix::kMSkewX);
fM[2] = m33 * m.get(SkMatrix::kMTransX);
fM[3] = m33 * m.get(SkMatrix::kMSkewY);
fM[4] = m33 * m.get(SkMatrix::kMScaleY);
fM[5] = m33 * m.get(SkMatrix::kMTransY);
} else {
fM[0] = m.get(SkMatrix::kMScaleX);
fM[1] = m.get(SkMatrix::kMSkewX);
fM[2] = m.get(SkMatrix::kMTransX);
fM[3] = m.get(SkMatrix::kMSkewY);
fM[4] = m.get(SkMatrix::kMScaleY);
fM[5] = m.get(SkMatrix::kMTransY);
}
fM[0] = (float)((0.5*a3 + a6)*scale);
fM[1] = (float)((0.5*a4 + a7)*scale);
fM[2] = (float)((0.5*a5 + a8)*scale);
fM[3] = (float)(a6*scale);
fM[4] = (float)(a7*scale);
fM[5] = (float)(a8*scale);
}
}