From 2db7fe7d3b7ee875e1099a22f0af17520696f5d7 Mon Sep 17 00:00:00 2001 From: "commit-bot@chromium.org" Date: Wed, 7 May 2014 15:31:40 +0000 Subject: 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 --- tests/PathOpsCubicLineIntersectionIdeas.cpp | 283 ++++++++++++++++++++++++++++ 1 file changed, 283 insertions(+) create mode 100644 tests/PathOpsCubicLineIntersectionIdeas.cpp (limited to 'tests/PathOpsCubicLineIntersectionIdeas.cpp') diff --git a/tests/PathOpsCubicLineIntersectionIdeas.cpp b/tests/PathOpsCubicLineIntersectionIdeas.cpp new file mode 100644 index 0000000000..2887a2ccfc --- /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); +} -- cgit v1.2.3