PathOpsCubicLineIntersectionIdeas.cpp 10 KB

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