When solving the cubic line intersection directly fails, use binary search as a fallback.
The cubic line intersection math empirically works 99.99% of the time (fails 3100 out of 1B random tests) but when it fails, an intersection may be missed altogether.
The binary search is may not find a solution if the cubic line failed to find any solutions at all, but so far that case hasn't arisen.
BUG=skia:2504
TBR=reed@google.com
Author: caryclark@google.com
Review URL: https://codereview.chromium.org/266063003
git-svn-id: http://skia.googlecode.com/svn/trunk@14614 2bbb7eff-a529-9590-31e7-b0007b416f81
diff --git a/tests/PathOpsCubicLineIntersectionIdeas.cpp b/tests/PathOpsCubicLineIntersectionIdeas.cpp
new file mode 100644
index 0000000..2887a2c
--- /dev/null
+++ b/tests/PathOpsCubicLineIntersectionIdeas.cpp
@@ -0,0 +1,283 @@
+/*
+ * Copyright 2014 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+#include "PathOpsTestCommon.h"
+#include "SkIntersections.h"
+#include "SkPathOpsCubic.h"
+#include "SkPathOpsLine.h"
+#include "SkPathOpsQuad.h"
+#include "SkRandom.h"
+#include "SkReduceOrder.h"
+#include "Test.h"
+
+static bool gPathOpsCubicLineIntersectionIdeasVerbose = false;
+
+static struct CubicLineFailures {
+ SkDCubic c;
+ double t;
+ SkDPoint p;
+} cubicLineFailures[] = {
+ {{{{-164.3726806640625, 36.826904296875}, {-189.045166015625, -953.2220458984375},
+ {926.505859375, -897.36175537109375}, {-139.33489990234375, 204.40771484375}}},
+ 0.37329583, {107.54935269006289, -632.13736293162208}},
+ {{{{784.056884765625, -554.8350830078125}, {67.5489501953125, 509.0224609375},
+ {-447.713134765625, 751.375}, {415.7784423828125, 172.22265625}}},
+ 0.660005242, {-32.973148967736151, 478.01341797403569}},
+ {{{{-580.6834716796875, -127.044921875}, {-872.8983154296875, -945.54302978515625},
+ {260.8092041015625, -909.34991455078125}, {-976.2125244140625, -18.46551513671875}}},
+ 0.578826774, {-390.17910153915489, -687.21144412296007}},
+};
+
+int cubicLineFailuresCount = (int) SK_ARRAY_COUNT(cubicLineFailures);
+
+double measuredSteps[] = {
+ 9.15910731e-007, 8.6600277e-007, 7.4122059e-007, 6.92087618e-007, 8.35290245e-007,
+ 3.29763199e-007, 5.07547773e-007, 4.41294224e-007, 0, 0,
+ 3.76879167e-006, 1.06126249e-006, 2.36873967e-006, 1.62421134e-005, 3.09103599e-005,
+ 4.38917976e-005, 0.000112348938, 0.000243149242, 0.000433174114, 0.00170880232,
+ 0.00272619724, 0.00518844604, 0.000352621078, 0.00175960064, 0.027875185,
+ 0.0351329803, 0.103964925,
+};
+
+/* last output : errors=3121
+ 9.1796875e-007 8.59375e-007 7.5e-007 6.875e-007 8.4375e-007
+ 3.125e-007 5e-007 4.375e-007 0 0
+ 3.75e-006 1.09375e-006 2.1875e-006 1.640625e-005 3.0859375e-005
+ 4.38964844e-005 0.000112304687 0.000243164063 0.000433181763 0.00170898437
+ 0.00272619247 0.00518844604 0.000352621078 0.00175960064 0.027875185
+ 0.0351329803 0.103964925
+*/
+
+static double binary_search(const SkDCubic& cubic, double step, const SkDPoint& pt, double t,
+ int* iters) {
+ double firstStep = step;
+ do {
+ *iters += 1;
+ SkDPoint cubicAtT = cubic.ptAtT(t);
+ if (cubicAtT.approximatelyEqual(pt)) {
+ break;
+ }
+ double calcX = cubicAtT.fX - pt.fX;
+ double calcY = cubicAtT.fY - pt.fY;
+ double calcDist = calcX * calcX + calcY * calcY;
+ if (step == 0) {
+ SkDebugf("binary search failed: step=%1.9g cubic=", firstStep);
+ cubic.dump();
+ SkDebugf(" t=%1.9g ", t);
+ pt.dump();
+ SkDebugf("\n");
+ return -1;
+ }
+ double lastStep = step;
+ step /= 2;
+ SkDPoint lessPt = cubic.ptAtT(t - lastStep);
+ double lessX = lessPt.fX - pt.fX;
+ double lessY = lessPt.fY - pt.fY;
+ double lessDist = lessX * lessX + lessY * lessY;
+ // use larger x/y difference to choose step
+ if (calcDist > lessDist) {
+ t -= step;
+ t = SkTMax(0., t);
+ } else {
+ SkDPoint morePt = cubic.ptAtT(t + lastStep);
+ double moreX = morePt.fX - pt.fX;
+ double moreY = morePt.fY - pt.fY;
+ double moreDist = moreX * moreX + moreY * moreY;
+ if (calcDist <= moreDist) {
+ continue;
+ }
+ t += step;
+ t = SkTMin(1., t);
+ }
+ } while (true);
+ return t;
+}
+
+#if 0
+static bool r2check(double A, double B, double C, double D, double* R2MinusQ3Ptr) {
+ if (approximately_zero(A)
+ && approximately_zero_when_compared_to(A, B)
+ && approximately_zero_when_compared_to(A, C)
+ && approximately_zero_when_compared_to(A, D)) { // we're just a quadratic
+ return false;
+ }
+ if (approximately_zero_when_compared_to(D, A)
+ && approximately_zero_when_compared_to(D, B)
+ && approximately_zero_when_compared_to(D, C)) { // 0 is one root
+ return false;
+ }
+ if (approximately_zero(A + B + C + D)) { // 1 is one root
+ return false;
+ }
+ double a, b, c;
+ {
+ double invA = 1 / A;
+ a = B * invA;
+ b = C * invA;
+ c = D * invA;
+ }
+ double a2 = a * a;
+ double Q = (a2 - b * 3) / 9;
+ double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
+ double R2 = R * R;
+ double Q3 = Q * Q * Q;
+ double R2MinusQ3 = R2 - Q3;
+ *R2MinusQ3Ptr = R2MinusQ3;
+ return true;
+}
+#endif
+
+/* What is the relationship between the accuracy of the root in range and the magnitude of all
+ roots? To find out, create a bunch of cubics, and measure */
+
+DEF_TEST(PathOpsCubicLineRoots, reporter) {
+ if (!gPathOpsCubicLineIntersectionIdeasVerbose) { // slow; exclude it by default
+ return;
+ }
+ SkRandom ran;
+ double worstStep[256] = {0};
+ int errors = 0;
+ int iters = 0;
+ double smallestR2 = 0;
+ double largestR2 = 0;
+ for (int index = 0; index < 1000000000; ++index) {
+ SkDPoint origin = {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)};
+ SkDCubic cubic = {{origin,
+ {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
+ {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
+ {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}
+ }};
+ // construct a line at a known intersection
+ double t = ran.nextRangeF(0, 1);
+ SkDPoint pt = cubic.ptAtT(t);
+ // skip answers with no intersections (although note the bug!) or two, or more
+ // see if the line / cubic has a fun range of roots
+ double A, B, C, D;
+ SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D);
+ D -= pt.fY;
+ double allRoots[3] = {0}, validRoots[3] = {0};
+ int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots);
+ int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots);
+ if (valid != 1) {
+ continue;
+ }
+ if (realRoots == 1) {
+ continue;
+ }
+ t = validRoots[0];
+ SkDPoint calcPt = cubic.ptAtT(t);
+ if (calcPt.approximatelyEqual(pt)) {
+ continue;
+ }
+#if 0
+ double R2MinusQ3;
+ if (r2check(A, B, C, D, &R2MinusQ3)) {
+ smallestR2 = SkTMin(smallestR2, R2MinusQ3);
+ largestR2 = SkTMax(largestR2, R2MinusQ3);
+ }
+#endif
+ double largest = SkTMax(fabs(allRoots[0]), fabs(allRoots[1]));
+ if (realRoots == 3) {
+ largest = SkTMax(largest, fabs(allRoots[2]));
+ }
+ int largeBits;
+ if (largest <= 1) {
+#if 0
+ SkDebugf("realRoots=%d (%1.9g, %1.9g, %1.9g) valid=%d (%1.9g, %1.9g, %1.9g)\n",
+ realRoots, allRoots[0], allRoots[1], allRoots[2], valid, validRoots[0],
+ validRoots[1], validRoots[2]);
+#endif
+ double smallest = SkTMin(allRoots[0], allRoots[1]);
+ if (realRoots == 3) {
+ smallest = SkTMin(smallest, allRoots[2]);
+ }
+ SK_ALWAYSBREAK(smallest < 0);
+ SK_ALWAYSBREAK(smallest >= -1);
+ largeBits = 0;
+ } else {
+ frexp(largest, &largeBits);
+ SK_ALWAYSBREAK(largeBits >= 0);
+ SK_ALWAYSBREAK(largeBits < 256);
+ }
+ double step = 1e-6;
+ if (largeBits > 21) {
+ step = 1e-1;
+ } else if (largeBits > 18) {
+ step = 1e-2;
+ } else if (largeBits > 15) {
+ step = 1e-3;
+ } else if (largeBits > 12) {
+ step = 1e-4;
+ } else if (largeBits > 9) {
+ step = 1e-5;
+ }
+ double diff;
+ do {
+ double newT = binary_search(cubic, step, pt, t, &iters);
+ if (newT >= 0) {
+ diff = fabs(t - newT);
+ break;
+ }
+ step *= 1.5;
+ SK_ALWAYSBREAK(step < 1);
+ } while (true);
+ worstStep[largeBits] = SkTMax(worstStep[largeBits], diff);
+#if 0
+ {
+ cubic.dump();
+ SkDebugf("\n");
+ SkDLine line = {{{pt.fX - 1, pt.fY}, {pt.fX + 1, pt.fY}}};
+ line.dump();
+ SkDebugf("\n");
+ }
+#endif
+ ++errors;
+ }
+ SkDebugf("errors=%d avgIter=%1.9g", errors, (double) iters / errors);
+ SkDebugf(" steps: ");
+ int worstLimit = SK_ARRAY_COUNT(worstStep);
+ while (worstStep[--worstLimit] == 0) ;
+ for (int idx2 = 0; idx2 <= worstLimit; ++idx2) {
+ SkDebugf("%1.9g ", worstStep[idx2]);
+ }
+ SkDebugf("\n");
+ SkDebugf("smallestR2=%1.9g largestR2=%1.9g\n", smallestR2, largestR2);
+}
+
+static double testOneFailure(const CubicLineFailures& failure) {
+ const SkDCubic& cubic = failure.c;
+ const SkDPoint& pt = failure.p;
+ double A, B, C, D;
+ SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D);
+ D -= pt.fY;
+ double allRoots[3] = {0}, validRoots[3] = {0};
+ int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots);
+ int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots);
+ SK_ALWAYSBREAK(valid == 1);
+ SK_ALWAYSBREAK(realRoots != 1);
+ double t = validRoots[0];
+ SkDPoint calcPt = cubic.ptAtT(t);
+ SK_ALWAYSBREAK(!calcPt.approximatelyEqual(pt));
+ int iters = 0;
+ double newT = binary_search(cubic, 0.1, pt, t, &iters);
+ return newT;
+}
+
+DEF_TEST(PathOpsCubicLineFailures, reporter) {
+ return; // disable for now
+ for (int index = 0; index < cubicLineFailuresCount; ++index) {
+ const CubicLineFailures& failure = cubicLineFailures[index];
+ double newT = testOneFailure(failure);
+ SK_ALWAYSBREAK(newT >= 0);
+ }
+}
+
+DEF_TEST(PathOpsCubicLineOneFailure, reporter) {
+ return; // disable for now
+ const CubicLineFailures& failure = cubicLineFailures[1];
+ double newT = testOneFailure(failure);
+ SK_ALWAYSBREAK(newT >= 0);
+}
diff --git a/tests/PathOpsCubicLineIntersectionTest.cpp b/tests/PathOpsCubicLineIntersectionTest.cpp
index 1a2e188..234a538 100644
--- a/tests/PathOpsCubicLineIntersectionTest.cpp
+++ b/tests/PathOpsCubicLineIntersectionTest.cpp
@@ -49,6 +49,13 @@
}
static lineCubic lineCubicTests[] = {
+ {{{{-634.60540771484375, -481.262939453125}, {266.2696533203125, -752.70867919921875},
+ {-751.8370361328125, -317.37921142578125}, {-969.7427978515625, 824.7255859375}}},
+ {{{-287.9506133720805678, -557.1376476615772617},
+ {-285.9506133720805678, -557.1376476615772617}}}},
+
+ {{{{36.7184372,0.888650894}, {36.7184372,0.888650894}, {35.1233864,0.554015458},
+ {34.5114098,-0.115255356}}}, {{{35.4531212,0}, {31.9375,0}}}},
{{{{421, 378}, {421, 380.209137f}, {418.761414f, 382}, {416, 382}}},
{{{320, 378}, {421, 378.000031f}}}},
@@ -83,6 +90,32 @@
static const size_t lineCubicTests_count = SK_ARRAY_COUNT(lineCubicTests);
+static int doIntersect(SkIntersections& intersections, const SkDCubic& cubic, const SkDLine& line) {
+ int result;
+ bool flipped = false;
+ if (line[0].fX == line[1].fX) {
+ double top = line[0].fY;
+ double bottom = line[1].fY;
+ flipped = top > bottom;
+ if (flipped) {
+ SkTSwap<double>(top, bottom);
+ }
+ result = intersections.vertical(cubic, top, bottom, line[0].fX, flipped);
+ } else if (line[0].fY == line[1].fY) {
+ double left = line[0].fX;
+ double right = line[1].fX;
+ flipped = left > right;
+ if (flipped) {
+ SkTSwap<double>(left, right);
+ }
+ result = intersections.horizontal(cubic, left, right, line[0].fY, flipped);
+ } else {
+ intersections.intersect(cubic, line);
+ result = intersections.used();
+ }
+ return result;
+}
+
static void testOne(skiatest::Reporter* reporter, int iIndex) {
const SkDCubic& cubic = lineCubicTests[iIndex].cubic;
SkASSERT(ValidCubic(cubic));
@@ -102,7 +135,7 @@
}
if (order1 == 4 && order2 == 2) {
SkIntersections i;
- int roots = i.intersect(cubic, line);
+ int roots = doIntersect(i, cubic, line);
for (int pt = 0; pt < roots; ++pt) {
double tt1 = i[0][pt];
SkDPoint xy1 = cubic.ptAtT(tt1);
diff --git a/tests/PathOpsOpTest.cpp b/tests/PathOpsOpTest.cpp
index 86baea4..75b6030 100644
--- a/tests/PathOpsOpTest.cpp
+++ b/tests/PathOpsOpTest.cpp
@@ -3329,10 +3329,30 @@
testPathOp(reporter, path1, path2, kDifference_PathOp, filename);
}
+static void issue2504(skiatest::Reporter* reporter, const char* filename) {
+ SkPath path1;
+ path1.moveTo(34.2421875, -5.976562976837158203125);
+ path1.lineTo(35.453121185302734375, 0);
+ path1.lineTo(31.9375, 0);
+ path1.close();
+
+ SkPath path2;
+ path2.moveTo(36.71843719482421875, 0.8886508941650390625);
+ path2.cubicTo(36.71843719482421875, 0.8886508941650390625,
+ 35.123386383056640625, 0.554015457630157470703125,
+ 34.511409759521484375, -0.1152553558349609375);
+ path2.cubicTo(33.899425506591796875, -0.7845261096954345703125,
+ 34.53484344482421875, -5.6777553558349609375,
+ 34.53484344482421875, -5.6777553558349609375);
+ path2.close();
+ testPathOp(reporter, path1, path2, kUnion_PathOp, filename);
+}
+
static void (*firstTest)(skiatest::Reporter* , const char* filename) = 0;
static void (*stopTest)(skiatest::Reporter* , const char* filename) = 0;
static struct TestDesc tests[] = {
+ TEST(issue2504),
TEST(kari1),
TEST(quadOp10i),
#if 0 // FIXME: serpentine curve is ordered the wrong way