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
binary_search(const SkDCubic & cubic,double step,const SkDPoint & pt,double t,int * iters)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
DEF_TEST(PathOpsCubicLineRoots,reporter)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
testOneFailure(const CubicLineFailures & failure)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
DEF_TEST(PathOpsCubicLineFailures,reporter)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
DEF_TEST(PathOpsCubicLineOneFailure,reporter)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