OLD | NEW |
(Empty) | |
| 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 = true; |
| 17 |
| 18 static struct CubicLineFailures { |
| 19 SkDCubic c; |
| 20 double t; |
| 21 SkDPoint p; |
| 22 } cubicLineFailures[] = { |
| 23 {{{{-164.3726806640625, 36.826904296875}, {-189.045166015625, -953.222045898
4375}, |
| 24 {926.505859375, -897.36175537109375}, {-139.33489990234375, 204.40771484
375}}}, |
| 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.543029785
15625}, |
| 30 {260.8092041015625, -909.34991455078125}, {-976.2125244140625, -18.46551
513671875}}}, |
| 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.35290245
e-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.091035
99e-005, |
| 40 4.38917976e-005, 0.000112348938, 0.000243149242, 0.000433174114, 0.001708802
32, |
| 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 diffMeasure(const SkDCubic& cubic, double t) { |
| 55 SkDPoint pt = cubic.ptAtT(t); |
| 56 // skip answers with no intersections (although note the bug!) or two, or mo
re |
| 57 // see if the line / cubic has a fun range of roots |
| 58 double A, B, C, D; |
| 59 SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D); |
| 60 D -= pt.fY; |
| 61 double allRoots[3], validRoots[3]; |
| 62 int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots); |
| 63 if (realRoots < 2) { |
| 64 return 0; |
| 65 } |
| 66 int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots); |
| 67 if (valid != 1) { |
| 68 return 0; |
| 69 } |
| 70 return fabs(t - validRoots[0]); |
| 71 } |
| 72 |
| 73 static double binary_search(const SkDCubic& cubic, double step, const SkDPoint&
pt, double t, |
| 74 int* iters) { |
| 75 double firstStep = step; |
| 76 do { |
| 77 *iters += 1; |
| 78 SkDPoint cubicAtT = cubic.ptAtT(t); |
| 79 if (cubicAtT.approximatelyEqual(pt)) { |
| 80 break; |
| 81 } |
| 82 double calcX = cubicAtT.fX - pt.fX; |
| 83 double calcY = cubicAtT.fY - pt.fY; |
| 84 double calcDist = calcX * calcX + calcY * calcY; |
| 85 if (step == 0) { |
| 86 SkDebugf("binary search failed: step=%1.9g cubic=", firstStep); |
| 87 cubic.dump(); |
| 88 SkDebugf(" t=%1.9g ", t); |
| 89 pt.dump(); |
| 90 SkDebugf("\n"); |
| 91 return -1; |
| 92 } |
| 93 double lastStep = step; |
| 94 step /= 2; |
| 95 SkDPoint lessPt = cubic.ptAtT(t - lastStep); |
| 96 double lessX = lessPt.fX - pt.fX; |
| 97 double lessY = lessPt.fY - pt.fY; |
| 98 double lessDist = lessX * lessX + lessY * lessY; |
| 99 // use larger x/y difference to choose step |
| 100 if (calcDist > lessDist) { |
| 101 t -= step; |
| 102 t = SkTMax(0., t); |
| 103 } else { |
| 104 SkDPoint morePt = cubic.ptAtT(t + lastStep); |
| 105 double moreX = morePt.fX - pt.fX; |
| 106 double moreY = morePt.fY - pt.fY; |
| 107 double moreDist = moreX * moreX + moreY * moreY; |
| 108 if (calcDist <= moreDist) { |
| 109 continue; |
| 110 } |
| 111 t += step; |
| 112 t = SkTMin(1., t); |
| 113 } |
| 114 } while (true); |
| 115 return t; |
| 116 } |
| 117 |
| 118 static bool r2check(double A, double B, double C, double D, double* R2MinusQ3Ptr
) { |
| 119 if (approximately_zero(A) |
| 120 && approximately_zero_when_compared_to(A, B) |
| 121 && approximately_zero_when_compared_to(A, C) |
| 122 && approximately_zero_when_compared_to(A, D)) { // we're just a qua
dratic |
| 123 return false; |
| 124 } |
| 125 if (approximately_zero_when_compared_to(D, A) |
| 126 && approximately_zero_when_compared_to(D, B) |
| 127 && approximately_zero_when_compared_to(D, C)) { // 0 is one root |
| 128 return false; |
| 129 } |
| 130 if (approximately_zero(A + B + C + D)) { // 1 is one root |
| 131 return false; |
| 132 } |
| 133 double a, b, c; |
| 134 { |
| 135 double invA = 1 / A; |
| 136 a = B * invA; |
| 137 b = C * invA; |
| 138 c = D * invA; |
| 139 } |
| 140 double a2 = a * a; |
| 141 double Q = (a2 - b * 3) / 9; |
| 142 double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; |
| 143 double R2 = R * R; |
| 144 double Q3 = Q * Q * Q; |
| 145 double R2MinusQ3 = R2 - Q3; |
| 146 *R2MinusQ3Ptr = R2MinusQ3; |
| 147 return true; |
| 148 } |
| 149 |
| 150 /* What is the relationship between the accuracy of the root in range and the ma
gnitude of all |
| 151 roots? To find out, create a bunch of cubics, and measure */ |
| 152 |
| 153 DEF_TEST(PathOpsCubicLineRoots, reporter) { |
| 154 if (!gPathOpsCubicLineIntersectionIdeasVerbose) { // takes a while to run -
- so exclude it by default |
| 155 return; |
| 156 } |
| 157 SkRandom ran; |
| 158 double worstStep[256] = {0}; |
| 159 int errors = 0; |
| 160 int iters = 0; |
| 161 double smallestR2 = 0; |
| 162 double largestR2 = 0; |
| 163 for (int index = 0; index < 1000000000; ++index) { |
| 164 SkDPoint origin = {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 10
00)}; |
| 165 SkDCubic cubic = {{origin, |
| 166 {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}, |
| 167 {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)}, |
| 168 {ran.nextRangeF(-1000, 1000), ran.nextRangeF(-1000, 1000)} |
| 169 }}; |
| 170 // construct a line at a known intersection |
| 171 double t = ran.nextRangeF(0, 1); |
| 172 SkDPoint pt = cubic.ptAtT(t); |
| 173 // skip answers with no intersections (although note the bug!) or two, o
r more |
| 174 // see if the line / cubic has a fun range of roots |
| 175 double A, B, C, D; |
| 176 SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D); |
| 177 D -= pt.fY; |
| 178 double allRoots[3] = {0}, validRoots[3] = {0}; |
| 179 int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots); |
| 180 int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots); |
| 181 if (valid != 1) { |
| 182 continue; |
| 183 } |
| 184 if (realRoots == 1) { |
| 185 continue; |
| 186 } |
| 187 t = validRoots[0]; |
| 188 SkDPoint calcPt = cubic.ptAtT(t); |
| 189 if (calcPt.approximatelyEqual(pt)) { |
| 190 continue; |
| 191 } |
| 192 #if 0 |
| 193 double R2MinusQ3; |
| 194 if (r2check(A, B, C, D, &R2MinusQ3)) { |
| 195 smallestR2 = SkTMin(smallestR2, R2MinusQ3); |
| 196 largestR2 = SkTMax(largestR2, R2MinusQ3); |
| 197 } |
| 198 #endif |
| 199 double largest = SkTMax(fabs(allRoots[0]), fabs(allRoots[1])); |
| 200 if (realRoots == 3) { |
| 201 largest = SkTMax(largest, fabs(allRoots[2])); |
| 202 } |
| 203 int largeBits; |
| 204 if (largest <= 1) { |
| 205 #if 0 |
| 206 SkDebugf("realRoots=%d (%1.9g, %1.9g, %1.9g) valid=%d (%1.9g, %1.9g,
%1.9g)\n", |
| 207 realRoots, allRoots[0], allRoots[1], allRoots[2], valid, validRo
ots[0], |
| 208 validRoots[1], validRoots[2]); |
| 209 #endif |
| 210 double smallest = SkTMin(allRoots[0], allRoots[1]); |
| 211 if (realRoots == 3) { |
| 212 smallest = SkTMin(smallest, allRoots[2]); |
| 213 } |
| 214 SK_ALWAYSBREAK(smallest < 0); |
| 215 SK_ALWAYSBREAK(smallest >= -1); |
| 216 largeBits = 0; |
| 217 } else { |
| 218 frexp(largest, &largeBits); |
| 219 SK_ALWAYSBREAK(largeBits >= 0); |
| 220 SK_ALWAYSBREAK(largeBits < 256); |
| 221 } |
| 222 double step = 1e-6; |
| 223 if (largeBits > 21) { |
| 224 step = 1e-1; |
| 225 } else if (largeBits > 18) { |
| 226 step = 1e-2; |
| 227 } else if (largeBits > 15) { |
| 228 step = 1e-3; |
| 229 } else if (largeBits > 12) { |
| 230 step = 1e-4; |
| 231 } else if (largeBits > 9) { |
| 232 step = 1e-5; |
| 233 } |
| 234 double diff; |
| 235 do { |
| 236 double newT = binary_search(cubic, step, pt, t, &iters); |
| 237 if (newT >= 0) { |
| 238 diff = fabs(t - newT); |
| 239 break; |
| 240 } |
| 241 step *= 1.5; |
| 242 SK_ALWAYSBREAK(step < 1); |
| 243 } while (true); |
| 244 worstStep[largeBits] = SkTMax(worstStep[largeBits], diff); |
| 245 #if 0 |
| 246 { |
| 247 cubic.dump(); |
| 248 SkDebugf("\n"); |
| 249 SkDLine line = {{{pt.fX - 1, pt.fY}, {pt.fX + 1, pt.fY}}}; |
| 250 line.dump(); |
| 251 SkDebugf("\n"); |
| 252 } |
| 253 #endif |
| 254 ++errors; |
| 255 } |
| 256 SkDebugf("errors=%d avgIter=%1.9g", errors, (double) iters / errors); |
| 257 SkDebugf(" steps: "); |
| 258 int worstLimit = SK_ARRAY_COUNT(worstStep); |
| 259 while (worstStep[--worstLimit] == 0) ; |
| 260 for (int idx2 = 0; idx2 <= worstLimit; ++idx2) { |
| 261 SkDebugf("%1.9g ", worstStep[idx2]); |
| 262 } |
| 263 SkDebugf("\n"); |
| 264 SkDebugf("smallestR2=%1.9g largestR2=%1.9g\n", smallestR2, largestR2); |
| 265 } |
| 266 |
| 267 static double testOneFailure(const CubicLineFailures& failure) { |
| 268 const SkDCubic& cubic = failure.c; |
| 269 const SkDPoint& pt = failure.p; |
| 270 double A, B, C, D; |
| 271 SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D); |
| 272 D -= pt.fY; |
| 273 double allRoots[3] = {0}, validRoots[3] = {0}; |
| 274 int realRoots = SkDCubic::RootsReal(A, B, C, D, allRoots); |
| 275 int valid = SkDQuad::AddValidTs(allRoots, realRoots, validRoots); |
| 276 SK_ALWAYSBREAK(valid == 1); |
| 277 SK_ALWAYSBREAK(realRoots != 1); |
| 278 double t = validRoots[0]; |
| 279 SkDPoint calcPt = cubic.ptAtT(t); |
| 280 SK_ALWAYSBREAK(!calcPt.approximatelyEqual(pt)); |
| 281 int iters = 0; |
| 282 double newT = binary_search(cubic, 0.1, pt, t, &iters); |
| 283 return newT; |
| 284 } |
| 285 |
| 286 DEF_TEST(PathOpsCubicLineFailures, reporter) { |
| 287 for (int index = 0; index < cubicLineFailuresCount; ++index) { |
| 288 const CubicLineFailures& failure = cubicLineFailures[index]; |
| 289 double newT = testOneFailure(failure); |
| 290 SK_ALWAYSBREAK(newT >= 0); |
| 291 } |
| 292 } |
| 293 |
| 294 DEF_TEST(PathOpsCubicLineOneFailure, reporter) { |
| 295 const CubicLineFailures& failure = cubicLineFailures[1]; |
| 296 double newT = testOneFailure(failure); |
| 297 SK_ALWAYSBREAK(newT >= 0); |
| 298 } |
OLD | NEW |