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 CubicPts 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 CubicPts cuPts = {{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 SkDCubic cubic; 156 cubic.debugSet(cuPts.fPts); 157 SkDPoint pt = cubic.ptAtT(t); 158 // skip answers with no intersections (although note the bug!) or two, or more 159 // see if the line / cubic has a fun range of roots 160 double A, B, C, D; 161 SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D); 162 D -= pt.fY; 163 double allRoots[3] = {0}, validRoots[3] = {0}; 164 int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots); 165 int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots); 166 if (valid != 1) { 167 continue; 168 } 169 if (realRoots == 1) { 170 continue; 171 } 172 t = validRoots[0]; 173 SkDPoint calcPt = cubic.ptAtT(t); 174 if (calcPt.approximatelyEqual(pt)) { 175 continue; 176 } 177 #if 0 178 double R2MinusQ3; 179 if (r2check(A, B, C, D, &R2MinusQ3)) { 180 smallestR2 = SkTMin(smallestR2, R2MinusQ3); 181 largestR2 = SkTMax(largestR2, R2MinusQ3); 182 } 183 #endif 184 double largest = SkTMax(fabs(allRoots[0]), fabs(allRoots[1])); 185 if (realRoots == 3) { 186 largest = SkTMax(largest, fabs(allRoots[2])); 187 } 188 int largeBits; 189 if (largest <= 1) { 190 #if 0 191 SkDebugf("realRoots=%d (%1.9g, %1.9g, %1.9g) valid=%d (%1.9g, %1.9g, %1.9g)\n", 192 realRoots, allRoots[0], allRoots[1], allRoots[2], valid, validRoots[0], 193 validRoots[1], validRoots[2]); 194 #endif 195 double smallest = SkTMin(allRoots[0], allRoots[1]); 196 if (realRoots == 3) { 197 smallest = SkTMin(smallest, allRoots[2]); 198 } 199 SkASSERT_RELEASE(smallest < 0); 200 SkASSERT_RELEASE(smallest >= -1); 201 largeBits = 0; 202 } else { 203 frexp(largest, &largeBits); 204 SkASSERT_RELEASE(largeBits >= 0); 205 SkASSERT_RELEASE(largeBits < 256); 206 } 207 double step = 1e-6; 208 if (largeBits > 21) { 209 step = 1e-1; 210 } else if (largeBits > 18) { 211 step = 1e-2; 212 } else if (largeBits > 15) { 213 step = 1e-3; 214 } else if (largeBits > 12) { 215 step = 1e-4; 216 } else if (largeBits > 9) { 217 step = 1e-5; 218 } 219 double diff; 220 do { 221 double newT = binary_search(cubic, step, pt, t, &iters); 222 if (newT >= 0) { 223 diff = fabs(t - newT); 224 break; 225 } 226 step *= 1.5; 227 SkASSERT_RELEASE(step < 1); 228 } while (true); 229 worstStep[largeBits] = SkTMax(worstStep[largeBits], diff); 230 #if 0 231 { 232 cubic.dump(); 233 SkDebugf("\n"); 234 SkDLine line = {{{pt.fX - 1, pt.fY}, {pt.fX + 1, pt.fY}}}; 235 line.dump(); 236 SkDebugf("\n"); 237 } 238 #endif 239 ++errors; 240 } 241 SkDebugf("errors=%d avgIter=%1.9g", errors, (double) iters / errors); 242 SkDebugf(" steps: "); 243 int worstLimit = SK_ARRAY_COUNT(worstStep); 244 while (worstStep[--worstLimit] == 0) ; 245 for (int idx2 = 0; idx2 <= worstLimit; ++idx2) { 246 SkDebugf("%1.9g ", worstStep[idx2]); 247 } 248 SkDebugf("\n"); 249 SkDebugf("smallestR2=%1.9g largestR2=%1.9g\n", smallestR2, largestR2); 250 } 251 252 static double testOneFailure(const CubicLineFailures& failure) { 253 const CubicPts& c = failure.c; 254 SkDCubic cubic; 255 cubic.debugSet(c.fPts); 256 const SkDPoint& pt = failure.p; 257 double A, B, C, D; 258 SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D); 259 D -= pt.fY; 260 double allRoots[3] = {0}, validRoots[3] = {0}; 261 int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots); 262 int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots); 263 SkASSERT_RELEASE(valid == 1); 264 SkASSERT_RELEASE(realRoots != 1); 265 double t = validRoots[0]; 266 SkDPoint calcPt = cubic.ptAtT(t); 267 SkASSERT_RELEASE(!calcPt.approximatelyEqual(pt)); 268 int iters = 0; 269 double newT = binary_search(cubic, 0.1, pt, t, &iters); 270 return newT; 271 } 272 273 DEF_TEST(PathOpsCubicLineFailures, reporter) { 274 return; // disable for now 275 for (int index = 0; index < cubicLineFailuresCount; ++index) { 276 const CubicLineFailures& failure = cubicLineFailures[index]; 277 double newT = testOneFailure(failure); 278 SkASSERT_RELEASE(newT >= 0); 279 } 280 } 281 282 DEF_TEST(PathOpsCubicLineOneFailure, reporter) { 283 return; // disable for now 284 const CubicLineFailures& failure = cubicLineFailures[1]; 285 double newT = testOneFailure(failure); 286 SkASSERT_RELEASE(newT >= 0); 287 } 288
