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
16static bool gPathOpsCubicLineIntersectionIdeasVerbose = false;
17
18static 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
34int cubicLineFailuresCount = (int) SK_ARRAY_COUNT(cubicLineFailures);
35
36double 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
54static 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
100static 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
136DEF_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
250static 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
269DEF_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
278DEF_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}
284