From 2db7fe7d3b7ee875e1099a22f0af17520696f5d7 Mon Sep 17 00:00:00 2001 From: "commit-bot@chromium.org" Date: Wed, 7 May 2014 15:31:40 +0000 Subject: [PATCH] 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 --- gyp/pathops_unittest.gyp | 1 + src/pathops/SkDCubicLineIntersection.cpp | 78 ++++-- src/pathops/SkPathOpsCubic.cpp | 72 ++++- src/pathops/SkPathOpsCubic.h | 34 ++- src/pathops/SkPathOpsDebug.h | 2 + src/pathops/SkPathOpsTypes.cpp | 7 + src/pathops/SkPathOpsTypes.h | 4 +- tests/PathOpsCubicLineIntersectionIdeas.cpp | 283 ++++++++++++++++++++ tests/PathOpsCubicLineIntersectionTest.cpp | 35 ++- tests/PathOpsOpTest.cpp | 20 ++ tools/pathops_sorter.htm | 39 +++ 11 files changed, 538 insertions(+), 37 deletions(-) create mode 100644 tests/PathOpsCubicLineIntersectionIdeas.cpp diff --git a/gyp/pathops_unittest.gyp b/gyp/pathops_unittest.gyp index 51af1e9ac5..35eeabd209 100644 --- a/gyp/pathops_unittest.gyp +++ b/gyp/pathops_unittest.gyp @@ -23,6 +23,7 @@ ], 'sources': [ '../tests/PathOpsAngleIdeas.cpp', + '../tests/PathOpsCubicLineIntersectionIdeas.cpp', '../tests/PathOpsDebug.cpp', '../tests/PathOpsOpLoopThreadedTest.cpp', '../tests/PathOpsSkpClipTest.cpp', diff --git a/src/pathops/SkDCubicLineIntersection.cpp b/src/pathops/SkDCubicLineIntersection.cpp index da4b983d45..41828bb714 100644 --- a/src/pathops/SkDCubicLineIntersection.cpp +++ b/src/pathops/SkDCubicLineIntersection.cpp @@ -97,13 +97,27 @@ public: int intersectRay(double roots[3]) { double adj = fLine[1].fX - fLine[0].fX; double opp = fLine[1].fY - fLine[0].fY; - SkDCubic r; + SkDCubic c; for (int n = 0; n < 4; ++n) { - r[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp; + c[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp; } double A, B, C, D; - SkDCubic::Coefficients(&r[0].fX, &A, &B, &C, &D); - return SkDCubic::RootsValidT(A, B, C, D, roots); + SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D); + int count = SkDCubic::RootsValidT(A, B, C, D, roots); + for (int index = 0; index < count; ++index) { + SkDPoint calcPt = c.ptAtT(roots[index]); + if (!approximately_zero(calcPt.fX)) { + for (int n = 0; n < 4; ++n) { + c[n].fY = (fCubic[n].fY - fLine[0].fY) * opp + + (fCubic[n].fX - fLine[0].fX) * adj; + } + double extremeTs[6]; + int extrema = SkDCubic::FindExtrema(c[0].fX, c[1].fX, c[2].fX, c[3].fX, extremeTs); + count = c.searchRoots(extremeTs, extrema, 0, SkDCubic::kXAxis, roots); + break; + } + } + return count; } int intersect() { @@ -146,11 +160,21 @@ public: return fIntersections->used(); } - int horizontalIntersect(double axisIntercept, double roots[3]) { + static int HorizontalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) { double A, B, C, D; - SkDCubic::Coefficients(&fCubic[0].fY, &A, &B, &C, &D); + SkDCubic::Coefficients(&c[0].fY, &A, &B, &C, &D); D -= axisIntercept; - return SkDCubic::RootsValidT(A, B, C, D, roots); + int count = SkDCubic::RootsValidT(A, B, C, D, roots); + for (int index = 0; index < count; ++index) { + SkDPoint calcPt = c.ptAtT(roots[index]); + if (!approximately_equal(calcPt.fY, axisIntercept)) { + double extremeTs[6]; + int extrema = SkDCubic::FindExtrema(c[0].fY, c[1].fY, c[2].fY, c[3].fY, extremeTs); + count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kYAxis, roots); + break; + } + } + return count; } int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) { @@ -158,11 +182,13 @@ public: if (fAllowNear) { addNearHorizontalEndPoints(left, right, axisIntercept); } - double rootVals[3]; - int roots = horizontalIntersect(axisIntercept, rootVals); - for (int index = 0; index < roots; ++index) { - double cubicT = rootVals[index]; - SkDPoint pt = fCubic.ptAtT(cubicT); + double roots[3]; + int count = HorizontalIntersect(fCubic, axisIntercept, roots); + for (int index = 0; index < count; ++index) { + double cubicT = roots[index]; + SkDPoint pt; + pt.fX = fCubic.ptAtT(cubicT).fX; + pt.fY = axisIntercept; double lineT = (pt.fX - left) / (right - left); if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) { fIntersections->insert(cubicT, lineT, pt); @@ -174,11 +200,21 @@ public: return fIntersections->used(); } - int verticalIntersect(double axisIntercept, double roots[3]) { + static int VerticalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) { double A, B, C, D; - SkDCubic::Coefficients(&fCubic[0].fX, &A, &B, &C, &D); + SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D); D -= axisIntercept; - return SkDCubic::RootsValidT(A, B, C, D, roots); + int count = SkDCubic::RootsValidT(A, B, C, D, roots); + for (int index = 0; index < count; ++index) { + SkDPoint calcPt = c.ptAtT(roots[index]); + if (!approximately_equal(calcPt.fX, axisIntercept)) { + double extremeTs[6]; + int extrema = SkDCubic::FindExtrema(c[0].fX, c[1].fX, c[2].fX, c[3].fX, extremeTs); + count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kXAxis, roots); + break; + } + } + return count; } int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) { @@ -186,11 +222,13 @@ public: if (fAllowNear) { addNearVerticalEndPoints(top, bottom, axisIntercept); } - double rootVals[3]; - int roots = verticalIntersect(axisIntercept, rootVals); - for (int index = 0; index < roots; ++index) { - double cubicT = rootVals[index]; - SkDPoint pt = fCubic.ptAtT(cubicT); + double roots[3]; + int count = VerticalIntersect(fCubic, axisIntercept, roots); + for (int index = 0; index < count; ++index) { + double cubicT = roots[index]; + SkDPoint pt; + pt.fX = axisIntercept; + pt.fY = fCubic.ptAtT(cubicT).fY; double lineT = (pt.fY - top) / (bottom - top); if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) { fIntersections->insert(cubicT, lineT, pt); diff --git a/src/pathops/SkPathOpsCubic.cpp b/src/pathops/SkPathOpsCubic.cpp index a89604f94c..9d70d58ec1 100644 --- a/src/pathops/SkPathOpsCubic.cpp +++ b/src/pathops/SkPathOpsCubic.cpp @@ -9,9 +9,58 @@ #include "SkPathOpsLine.h" #include "SkPathOpsQuad.h" #include "SkPathOpsRect.h" +#include "SkTSort.h" const int SkDCubic::gPrecisionUnit = 256; // FIXME: test different values in test framework +// give up when changing t no longer moves point +// also, copy point rather than recompute it when it does change +double SkDCubic::binarySearch(double min, double max, double axisIntercept, + SearchAxis xAxis) const { + double t = (min + max) / 2; + double step = (t - min) / 2; + SkDPoint cubicAtT = ptAtT(t); + double calcPos = (&cubicAtT.fX)[xAxis]; + double calcDist = calcPos - axisIntercept; + do { + double priorT = t - step; + SkASSERT(priorT >= min); + SkDPoint lessPt = ptAtT(priorT); + if (approximately_equal(lessPt.fX, cubicAtT.fX) + && approximately_equal(lessPt.fY, cubicAtT.fY)) { + return -1; // binary search found no point at this axis intercept + } + double lessDist = (&lessPt.fX)[xAxis] - axisIntercept; +#if DEBUG_CUBIC_BINARY_SEARCH + SkDebugf("t=%1.9g calc=%1.9g dist=%1.9g step=%1.9g less=%1.9g\n", t, calcPos, calcDist, + step, lessDist); +#endif + double lastStep = step; + step /= 2; + if (calcDist > 0 ? calcDist > lessDist : calcDist < lessDist) { + t = priorT; + } else { + double nextT = t + lastStep; + SkASSERT(nextT <= max); + SkDPoint morePt = ptAtT(nextT); + if (approximately_equal(morePt.fX, cubicAtT.fX) + && approximately_equal(morePt.fY, cubicAtT.fY)) { + return -1; // binary search found no point at this axis intercept + } + double moreDist = (&morePt.fX)[xAxis] - axisIntercept; + if (calcDist > 0 ? calcDist <= moreDist : calcDist >= moreDist) { + continue; + } + t = nextT; + } + SkDPoint testAtT = ptAtT(t); + cubicAtT = testAtT; + calcPos = (&cubicAtT.fX)[xAxis]; + calcDist = calcPos - axisIntercept; + } while (!approximately_equal(calcPos, axisIntercept)); + return t; +} + // FIXME: cache keep the bounds and/or precision with the caller? double SkDCubic::calcPrecision() const { SkDRect dRect; @@ -93,6 +142,27 @@ bool SkDCubic::monotonicInY() const { && between(fPts[0].fY, fPts[2].fY, fPts[3].fY); } +int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept, + SearchAxis xAxis, double* validRoots) const { + extrema += findInflections(&extremeTs[extrema]); + extremeTs[extrema++] = 0; + extremeTs[extrema] = 1; + SkTQSort(extremeTs, extremeTs + extrema); + int validCount = 0; + for (int index = 0; index < extrema; ) { + double min = extremeTs[index]; + double max = extremeTs[++index]; + if (min == max) { + continue; + } + double newT = binarySearch(min, max, axisIntercept, xAxis); + if (newT >= 0) { + validRoots[validCount++] = newT; + } + } + return validCount; +} + bool SkDCubic::serpentine() const { #if 0 // FIXME: enabling this fixes cubicOp114 but breaks cubicOp58d and cubicOp53d double tValues[2]; @@ -210,7 +280,7 @@ int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) { } r = A - adiv3; *roots++ = r; - if (AlmostDequalUlps(R2, Q3)) { + if (AlmostDequalUlps((double) R2, (double) Q3)) { r = -A / 2 - adiv3; if (!AlmostDequalUlps(s[0], r)) { *roots++ = r; diff --git a/src/pathops/SkPathOpsCubic.h b/src/pathops/SkPathOpsCubic.h index 500319607d..1037cae4f7 100644 --- a/src/pathops/SkPathOpsCubic.h +++ b/src/pathops/SkPathOpsCubic.h @@ -19,21 +19,16 @@ struct SkDCubicPair { }; struct SkDCubic { - SkDPoint fPts[4]; - - void set(const SkPoint pts[4]) { - fPts[0] = pts[0]; - fPts[1] = pts[1]; - fPts[2] = pts[2]; - fPts[3] = pts[3]; - } - - static const int gPrecisionUnit; + enum SearchAxis { + kXAxis, + kYAxis + }; const SkDPoint& operator[](int n) const { SkASSERT(n >= 0 && n < 4); return fPts[n]; } SkDPoint& operator[](int n) { SkASSERT(n >= 0 && n < 4); return fPts[n]; } void align(int endIndex, int ctrlIndex, SkDPoint* dstPt) const; + double binarySearch(double min, double max, double axisIntercept, SearchAxis xAxis) const; double calcPrecision() const; SkDCubicPair chopAt(double t) const; bool clockwise() const; @@ -42,9 +37,9 @@ struct SkDCubic { SkDVector dxdyAtT(double t) const; bool endsAreExtremaInXOrY() const; static int FindExtrema(double a, double b, double c, double d, double tValue[2]); - int findInflections(double tValues[]) const; + int findInflections(double tValues[2]) const; - static int FindInflections(const SkPoint a[4], double tValues[]) { + static int FindInflections(const SkPoint a[4], double tValues[2]) { SkDCubic cubic; cubic.set(a); return cubic.findInflections(tValues); @@ -56,7 +51,18 @@ struct SkDCubic { SkDPoint ptAtT(double t) const; static int RootsReal(double A, double B, double C, double D, double t[3]); static int RootsValidT(const double A, const double B, const double C, double D, double s[3]); + + int searchRoots(double extremes[6], int extrema, double axisIntercept, + SearchAxis xAxis, double* validRoots) const; bool serpentine() const; + + void set(const SkPoint pts[4]) { + fPts[0] = pts[0]; + fPts[1] = pts[1]; + fPts[2] = pts[2]; + fPts[3] = pts[3]; + } + SkDCubic subDivide(double t1, double t2) const; static SkDCubic SubDivide(const SkPoint a[4], double t1, double t2) { @@ -81,6 +87,10 @@ struct SkDCubic { // utilities callable by the user from the debugger when the implementation code is linked in void dump() const; void dumpNumber() const; + + static const int gPrecisionUnit; + + SkDPoint fPts[4]; }; #endif diff --git a/src/pathops/SkPathOpsDebug.h b/src/pathops/SkPathOpsDebug.h index 39d5a6dda8..bb54039244 100644 --- a/src/pathops/SkPathOpsDebug.h +++ b/src/pathops/SkPathOpsDebug.h @@ -50,6 +50,7 @@ #define DEBUG_CHECK_TINY 0 #define DEBUG_CONCIDENT 0 #define DEBUG_CROSS 0 +#define DEBUG_CUBIC_BINARY_SEARCH 0 #define DEBUG_FLAT_QUADS 0 #define DEBUG_FLOW 0 #define DEBUG_LIMIT_WIND_SUM 0 @@ -84,6 +85,7 @@ #define DEBUG_CHECK_TINY 1 #define DEBUG_CONCIDENT 1 #define DEBUG_CROSS 01 +#define DEBUG_CUBIC_BINARY_SEARCH 1 #define DEBUG_FLAT_QUADS 0 #define DEBUG_FLOW 1 #define DEBUG_LIMIT_WIND_SUM 4 diff --git a/src/pathops/SkPathOpsTypes.cpp b/src/pathops/SkPathOpsTypes.cpp index dbed086fbd..2c8d778c69 100644 --- a/src/pathops/SkPathOpsTypes.cpp +++ b/src/pathops/SkPathOpsTypes.cpp @@ -102,6 +102,13 @@ bool AlmostDequalUlps(float a, float b) { return d_equal_ulps(a, b, UlpsEpsilon); } +bool AlmostDequalUlps(double a, double b) { + if (SkScalarIsFinite(a) || SkScalarIsFinite(b)) { + return AlmostDequalUlps(SkDoubleToScalar(a), SkDoubleToScalar(b)); + } + return fabs(a - b) / SkTMax(fabs(a), fabs(b)) < FLT_EPSILON * 16; +} + bool AlmostEqualUlps(float a, float b) { const int UlpsEpsilon = 16; return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon); diff --git a/src/pathops/SkPathOpsTypes.h b/src/pathops/SkPathOpsTypes.h index e8054ad476..4f8bd15344 100644 --- a/src/pathops/SkPathOpsTypes.h +++ b/src/pathops/SkPathOpsTypes.h @@ -30,9 +30,7 @@ inline bool AlmostEqualUlps(double a, double b) { // Use Almost Dequal when comparing should not special case denormalized values. bool AlmostDequalUlps(float a, float b); -inline bool AlmostDequalUlps(double a, double b) { - return AlmostDequalUlps(SkDoubleToScalar(a), SkDoubleToScalar(b)); -} +bool AlmostDequalUlps(double a, double b); bool NotAlmostEqualUlps(float a, float b); inline bool NotAlmostEqualUlps(double a, double b) { 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); +} diff --git a/tests/PathOpsCubicLineIntersectionTest.cpp b/tests/PathOpsCubicLineIntersectionTest.cpp index 1a2e188589..234a53805f 100644 --- a/tests/PathOpsCubicLineIntersectionTest.cpp +++ b/tests/PathOpsCubicLineIntersectionTest.cpp @@ -49,6 +49,13 @@ static void testFail(skiatest::Reporter* reporter, int iIndex) { } 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 lineCubic lineCubicTests[] = { 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(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(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 @@ static void testOne(skiatest::Reporter* reporter, int iIndex) { } 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 86baea423c..75b6030aa5 100644 --- a/tests/PathOpsOpTest.cpp +++ b/tests/PathOpsOpTest.cpp @@ -3329,10 +3329,30 @@ static void kari1(skiatest::Reporter* reporter, const char* filename) { 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 diff --git a/tools/pathops_sorter.htm b/tools/pathops_sorter.htm index 59ae2b6d90..41314f62a1 100644 --- a/tools/pathops_sorter.htm +++ b/tools/pathops_sorter.htm @@ -867,11 +867,50 @@ op intersect {{{900.0235595703125, 551.60284423828125}, {900.06072998046875, 551.29705810546875}, {900.15655517578125, 551.0157470703125}}} +
+{{{-634.60540771484375, -481.262939453125}, {266.2696533203125, -752.70867919921875}, {-751.8370361328125, -317.37921142578125}, {-969.7427978515625, 824.7255859375}}} +{{{-287.9506133720805678, -557.1376476615772617}, {-285.9506133720805678, -557.1376476615772617}}} +
+ +
+{{{-818.4456787109375, 248.218017578125}, {944.18505859375, -252.2330322265625}, {957.3946533203125, -45.43280029296875}, {-591.766357421875, 868.6187744140625}}} +{{{435.1963493079119871, -16.42683763243891093}, {437.1963493079119871, -16.42683763243891093}}} +
+ +
+{{{-818.4456787109375, 248.218017578125}, {944.18505859375, -252.2330322265625}, {957.3946533203125, -45.43280029296875}, {-591.766357421875, 868.6187744140625}}} +{{{397.5007682490800676, -17.35020084021140008}, {399.5007682490800676, -17.35020084021140008}}} +
+ +
+{{{-652.660888671875, -384.6475830078125}, {-551.7723388671875, -925.5025634765625}, {-321.06658935546875, -813.10345458984375}, {142.6982421875, -47.4503173828125}}} +{{{-478.4372049758064236, -717.868282575075682}, {-476.4372049758064236, -717.868282575075682}}} +
+ +
+{{{-954.4322509765625, 827.2216796875}, {-420.24017333984375, -7.80560302734375}, {799.134765625, -971.4295654296875}, {-556.23486328125, 344.400146484375}}} + +{{{58.57411390280688579, -302.8879316712078662}, {60.57411390280688579, -302.8879316712078662}}} +
+ +
+{{{-634.60540771484375, -481.262939453125}, {266.2696533203125, -752.70867919921875}, {-751.8370361328125, -317.37921142578125}, {-969.7427978515625, 824.7255859375}}} +{{{-287.95061337208057, -557.13764766157726}, {-285.95061337208057, -557.13764766157726}}} +{{{-308.65463091760211, -549.4520029924679} -308.65463091760211, -569.4520029924679 +
+ +