src/third_party/skia/tests/PathOpsCubicLineIntersectionIdeas.cpp

   1 /*
   2  * Copyright 2014 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 #include "PathOpsTestCommon.h"
   8 #include "SkIntersections.h"
   9 #include "SkPathOpsCubic.h"
  10 #include "SkPathOpsLine.h"
  11 #include "SkPathOpsQuad.h"
  12 #include "SkRandom.h"
  13 #include "SkReduceOrder.h"
  14 #include "Test.h"
  15
  16 static bool gPathOpsCubicLineIntersectionIdeasVerbose = false;
  17
  18 static struct CubicLineFailures {
  19     SkDCubic c;
  20     double t;
  21     SkDPoint p;
  22 } cubicLineFailures[] = {
  23     {{{{-164.3726806640625, 36.826904296875}, {-189.045166015625, -953.2220458984375},
  24         {926.505859375, -897.36175537109375}, {-139.33489990234375, 204.40771484375}}},
  25         0.37329583, {107.54935269006289, -632.13736293162208}},
  26     {{{{784.056884765625, -554.8350830078125}, {67.5489501953125, 509.0224609375},
  27         {-447.713134765625, 751.375}, {415.7784423828125, 172.22265625}}},
  28         0.660005242, {-32.973148967736151, 478.01341797403569}},
  29     {{{{-580.6834716796875, -127.044921875}, {-872.8983154296875, -945.54302978515625},
  30         {260.8092041015625, -909.34991455078125}, {-976.2125244140625, -18.46551513671875}}},
  31         0.578826774, {-390.17910153915489, -687.21144412296007}},
  32 };
  33
  34 int cubicLineFailuresCount = (int) SK_ARRAY_COUNT(cubicLineFailures);
  35
  36 double measuredSteps[] = {
  37     9.15910731e-007, 8.6600277e-007, 7.4122059e-007, 6.92087618e-007, 8.35290245e-007,
  38     3.29763199e-007, 5.07547773e-007, 4.41294224e-007, 0, 0,
  39     3.76879167e-006, 1.06126249e-006, 2.36873967e-006, 1.62421134e-005, 3.09103599e-005,
  40     4.38917976e-005, 0.000112348938, 0.000243149242, 0.000433174114, 0.00170880232,
  41     0.00272619724, 0.00518844604, 0.000352621078, 0.00175960064, 0.027875185,
  42     0.0351329803, 0.103964925,
  43 };
  44
  45 /* last output : errors=3121
  46     9.1796875e-007 8.59375e-007 7.5e-007 6.875e-007 8.4375e-007
  47     3.125e-007 5e-007 4.375e-007 0 0
  48     3.75e-006 1.09375e-006 2.1875e-006 1.640625e-005 3.0859375e-005
  49     4.38964844e-005 0.000112304687 0.000243164063 0.000433181763 0.00170898437
  50     0.00272619247 0.00518844604 0.000352621078 0.00175960064 0.027875185
  51     0.0351329803 0.103964925
  52 */
  53
  54 static double binary_search(const SkDCubic& cubic, double step, const SkDPoint& pt, double t,
  55         int* iters) {
  56     double firstStep = step;
  57     do {
  58         *iters += 1;
  59         SkDPoint cubicAtT = cubic.ptAtT(t);
  60         if (cubicAtT.approximatelyEqual(pt)) {
  61             break;
  62         }
  63         double calcX = cubicAtT.fX - pt.fX;
  64         double calcY = cubicAtT.fY - pt.fY;
  65         double calcDist = calcX * calcX + calcY * calcY;
  66         if (step == 0) {
  67             SkDebugf("binary search failed: step=%1.9g cubic=", firstStep);
  68             cubic.dump();
  69             SkDebugf(" t=%1.9g ", t);
  70             pt.dump();
  71             SkDebugf("\n");
  72             return -1;
  73         }
  74         double lastStep = step;
  75         step /= 2;
  76         SkDPoint lessPt = cubic.ptAtT(t - lastStep);
  77         double lessX = lessPt.fX - pt.fX;
  78         double lessY = lessPt.fY - pt.fY;
  79         double lessDist = lessX * lessX + lessY * lessY;
  80         // use larger x/y difference to choose step
  81         if (calcDist > lessDist) {
  82             t -= step;
  83             t = SkTMax(0., t);
  84         } else {
  85             SkDPoint morePt = cubic.ptAtT(t + lastStep);
  86             double moreX = morePt.fX - pt.fX;
  87             double moreY = morePt.fY - pt.fY;
  88             double moreDist = moreX * moreX + moreY * moreY;
  89             if (calcDist <= moreDist) {
  90                 continue;
  91             }
  92             t += step;
  93             t = SkTMin(1., t);
  94         }
  95     } while (true);
  96     return t;
  97 }
  98
  99 #if 0
 100 static bool r2check(double A, double B, double C, double D, double* R2MinusQ3Ptr) {
 101     if (approximately_zero(A)
 102             && approximately_zero_when_compared_to(A, B)
 103             && approximately_zero_when_compared_to(A, C)
 104             && approximately_zero_when_compared_to(A, D)) {  // we're just a quadratic
 105         return false;
 106     }
 107     if (approximately_zero_when_compared_to(D, A)
 108             && approximately_zero_when_compared_to(D, B)
 109             && approximately_zero_when_compared_to(D, C)) {  // 0 is one root
 110         return false;
 111     }
 112     if (approximately_zero(A + B + C + D)) {  // 1 is one root
 113         return false;
 114     }
 115     double a, b, c;
 116     {
 117         double invA = 1 / A;
 118         a = B * invA;
 119         b = C * invA;
 120         c = D * invA;
 121     }
 122     double a2 = a * a;
 123     double Q = (a2 - b * 3) / 9;
 124     double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
 125     double R2 = R * R;
 126     double Q3 = Q * Q * Q;
 127     double R2MinusQ3 = R2 - Q3;
 128     *R2MinusQ3Ptr = R2MinusQ3;
 129     return true;
 130 }
 131 #endif
 132
 133 /* What is the relationship between the accuracy of the root in range and the magnitude of all
 134    roots? To find out, create a bunch of cubics, and measure */
 135
 136 DEF_TEST(PathOpsCubicLineRoots, reporter) {
 137     if (!gPathOpsCubicLineIntersectionIdeasVerbose) {  // slow; exclude it by default
 138         return;
 139     }
 140     SkRandom ran;
 141     double worstStep[256] = {0};
 142     int errors = 0;
 143     int iters = 0;
 144     double smallestR2 = 0;
 145     double largestR2 = 0;
 146     for (int index = 0; index < 1000000000; ++index) {
 147         SkDPoint origin = {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)};
 148         SkDCubic cubic = {{origin,
 149                 {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
 150                 {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)},
 151                 {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}
 152         }};
 153         // construct a line at a known intersection
 154         double t = ran.nextRangeF(0, 1);
 155         SkDPoint pt = cubic.ptAtT(t);
 156         // skip answers with no intersections (although note the bug!) or two, or more
 157         // see if the line / cubic has a fun range of roots
 158         double A, B, C, D;
 159         SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D);
 160         D -= pt.fY;
 161         double allRoots[3] = {0}, validRoots[3] = {0};
 162         int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots);
 163         int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots);
 164         if (valid != 1) {
 165             continue;
 166         }
 167         if (realRoots == 1) {
 168             continue;
 169         }
 170         t = validRoots[0];
 171         SkDPoint calcPt = cubic.ptAtT(t);
 172         if (calcPt.approximatelyEqual(pt)) {
 173             continue;
 174         }
 175 #if 0
 176         double R2MinusQ3;
 177         if (r2check(A, B, C, D, &R2MinusQ3)) {
 178             smallestR2 = SkTMin(smallestR2, R2MinusQ3);
 179             largestR2 = SkTMax(largestR2, R2MinusQ3);
 180         }
 181 #endif
 182         double largest = SkTMax(fabs(allRoots[0]), fabs(allRoots[1]));
 183         if (realRoots == 3) {
 184             largest = SkTMax(largest, fabs(allRoots[2]));
 185         }
 186         int largeBits;
 187         if (largest <= 1) {
 188 #if 0
 189             SkDebugf("realRoots=%d (%1.9g, %1.9g, %1.9g) valid=%d (%1.9g, %1.9g, %1.9g)\n",
 190                 realRoots, allRoots[0], allRoots[1], allRoots[2], valid, validRoots[0],
 191                 validRoots[1], validRoots[2]);
 192 #endif
 193             double smallest = SkTMin(allRoots[0], allRoots[1]);
 194             if (realRoots == 3) {
 195                 smallest = SkTMin(smallest, allRoots[2]);
 196             }
 197             SK_ALWAYSBREAK(smallest < 0);
 198             SK_ALWAYSBREAK(smallest >= -1);
 199             largeBits = 0;
 200         } else {
 201             frexp(largest, &largeBits);
 202             SK_ALWAYSBREAK(largeBits >= 0);
 203             SK_ALWAYSBREAK(largeBits < 256);
 204         }
 205         double step = 1e-6;
 206         if (largeBits > 21) {
 207             step = 1e-1;
 208         } else if (largeBits > 18) {
 209             step = 1e-2;
 210         } else if (largeBits > 15) {
 211             step = 1e-3;
 212         } else if (largeBits > 12) {
 213             step = 1e-4;
 214         } else if (largeBits > 9) {
 215             step = 1e-5;
 216         }
 217         double diff;
 218         do {
 219             double newT = binary_search(cubic, step, pt, t, &iters);
 220             if (newT >= 0) {
 221                 diff = fabs(t - newT);
 222                 break;
 223             }
 224             step *= 1.5;
 225             SK_ALWAYSBREAK(step < 1);
 226         } while (true);
 227         worstStep[largeBits] = SkTMax(worstStep[largeBits], diff);
 228 #if 0
 229         {
 230             cubic.dump();
 231             SkDebugf("\n");
 232             SkDLine line = {{{pt.fX - 1, pt.fY}, {pt.fX + 1, pt.fY}}};
 233             line.dump();
 234             SkDebugf("\n");
 235         }
 236 #endif
 237         ++errors;
 238     }
 239     SkDebugf("errors=%d avgIter=%1.9g", errors, (double) iters / errors);
 240     SkDebugf(" steps: ");
 241     int worstLimit = SK_ARRAY_COUNT(worstStep);
 242     while (worstStep[--worstLimit] == 0) ;
 243     for (int idx2 = 0; idx2 <= worstLimit; ++idx2) {
 244         SkDebugf("%1.9g ", worstStep[idx2]);
 245     }
 246     SkDebugf("\n");
 247     SkDebugf("smallestR2=%1.9g largestR2=%1.9g\n", smallestR2, largestR2);
 248 }
 249
 250 static double testOneFailure(const CubicLineFailures& failure) {
 251     const SkDCubic& cubic = failure.c;
 252     const SkDPoint& pt = failure.p;
 253     double A, B, C, D;
 254     SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D);
 255     D -= pt.fY;
 256     double allRoots[3] = {0}, validRoots[3] = {0};
 257     int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots);
 258     int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots);
 259     SK_ALWAYSBREAK(valid == 1);
 260     SK_ALWAYSBREAK(realRoots != 1);
 261     double t = validRoots[0];
 262     SkDPoint calcPt = cubic.ptAtT(t);
 263     SK_ALWAYSBREAK(!calcPt.approximatelyEqual(pt));
 264     int iters = 0;
 265     double newT = binary_search(cubic, 0.1, pt, t, &iters);
 266     return newT;
 267 }
 268
 269 DEF_TEST(PathOpsCubicLineFailures, reporter) {
 270     return;  // disable for now
 271     for (int index = 0; index < cubicLineFailuresCount; ++index) {
 272         const CubicLineFailures& failure = cubicLineFailures[index];
 273         double newT = testOneFailure(failure);
 274         SK_ALWAYSBREAK(newT >= 0);
 275     }
 276 }
 277
 278 DEF_TEST(PathOpsCubicLineOneFailure, reporter) {
 279     return;  // disable for now
 280     const CubicLineFailures& failure = cubicLineFailures[1];
 281     double newT = testOneFailure(failure);
 282     SK_ALWAYSBREAK(newT >= 0);
 283 }