1/*
2 * Copyright 2012 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 "CurveIntersection.h"
8#include "CurveUtilities.h"
9#include "CubicIntersection_TestData.h"
10#include "Intersection_Tests.h"
11#include "Intersections.h"
12#include "TestUtilities.h"
13
14#define SHOW_ORIGINAL 1
15
16const int firstCubicIntersectionTest = 9;
17
18static void standardTestCases() {
19    for (size_t index = firstCubicIntersectionTest; index < tests_count; ++index) {
20        const Cubic& cubic1 = tests[index][0];
21        const Cubic& cubic2 = tests[index][1];
22        Cubic reduce1, reduce2;
23        int order1 = reduceOrder(cubic1, reduce1, kReduceOrder_NoQuadraticsAllowed,
24            kReduceOrder_TreatAsFill);
25        int order2 = reduceOrder(cubic2, reduce2, kReduceOrder_NoQuadraticsAllowed,
26            kReduceOrder_TreatAsFill);
27        if (order1 < 4) {
28            printf("%s [%d] cubic1 order=%d\n", __FUNCTION__, (int) index, order1);
29            continue;
30        }
31        if (order2 < 4) {
32            printf("%s [%d] cubic2 order=%d\n", __FUNCTION__, (int) index, order2);
33            continue;
34        }
35        if (implicit_matches(reduce1, reduce2)) {
36            printf("%s [%d] coincident\n", __FUNCTION__, (int) index);
37            continue;
38        }
39        Intersections tIntersections;
40        intersect(reduce1, reduce2, tIntersections);
41        if (!tIntersections.intersected()) {
42            printf("%s [%d] no intersection\n", __FUNCTION__, (int) index);
43            continue;
44        }
45        for (int pt = 0; pt < tIntersections.used(); ++pt) {
46            double tt1 = tIntersections.fT[0][pt];
47            double tx1, ty1;
48            xy_at_t(cubic1, tt1, tx1, ty1);
49            double tt2 = tIntersections.fT[1][pt];
50            double tx2, ty2;
51            xy_at_t(cubic2, tt2, tx2, ty2);
52            if (!AlmostEqualUlps(tx1, tx2)) {
53                printf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
54                    __FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
55            }
56            if (!AlmostEqualUlps(ty1, ty2)) {
57                printf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
58                    __FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
59            }
60        }
61    }
62}
63
64static const Cubic testSet[] = {
65{{67.426548091427676, 37.993772624988935}, {23.483695892376684, 90.476863174921306}, {35.597065061143162, 79.872482633158796}, {75.38634169631932, 18.244890038969412}},
66{{67.4265481, 37.9937726}, {23.4836959, 90.4768632}, {35.5970651, 79.8724826}, {75.3863417, 18.24489}},
67
68{{0, 0}, {0, 1}, {1, 1}, {1, 0}},
69{{1, 0}, {0, 0}, {0, 1}, {1, 1}},
70
71{{0,1}, {4,5}, {1,0}, {5,3}},
72{{0,1}, {3,5}, {1,0}, {5,4}},
73
74{{0, 1}, {1, 6}, {1, 0}, {1, 0}},
75{{0, 1}, {0, 1}, {1, 0}, {6, 1}},
76
77{{0,1}, {3,4}, {1,0}, {5,1}},
78{{0,1}, {1,5}, {1,0}, {4,3}},
79
80{{0,1}, {1,2}, {1,0}, {6,1}},
81{{0,1}, {1,6}, {1,0}, {2,1}},
82
83{{0,1}, {0,5}, {1,0}, {4,0}},
84{{0,1}, {0,4}, {1,0}, {5,0}},
85
86{{0,1}, {3,4}, {1,0}, {3,0}},
87{{0,1}, {0,3}, {1,0}, {4,3}},
88
89{{0, 0}, {1, 2}, {3, 4}, {4, 4}},
90{{0, 0}, {1, 2}, {3, 4}, {4, 4}},
91{{4, 4}, {3, 4}, {1, 2}, {0, 0}},
92
93{{0,1}, {2,3}, {1,0}, {1,0}},
94{{0,1}, {0,1}, {1,0}, {3,2}},
95
96{{0,2}, {0,1}, {1,0}, {1,0}},
97{{0,1}, {0,1}, {2,0}, {1,0}},
98
99{{0, 1}, {0, 2}, {1, 0}, {1, 0}},
100{{0, 1}, {0, 1}, {1, 0}, {2, 0}},
101
102{{0, 1}, {1, 6}, {1, 0}, {2, 0}},
103{{0, 1}, {0, 2}, {1, 0}, {6, 1}},
104
105{{0, 1}, {5, 6}, {1, 0}, {1, 0}},
106{{0, 1}, {0, 1}, {1, 0}, {6, 5}},
107
108{{95.837747722788592, 45.025976907939643}, {16.564570095652982, 0.72959763963222402}, {63.209855865319199, 68.047528419665767}, {57.640240647662544, 59.524565264361243}},
109{{51.593891741518817, 38.53849970667553}, {62.34752929878772, 74.924924725166022}, {74.810149322641152, 34.17966562983564}, {29.368398119401373, 94.66719277886078}},
110
111{{39.765160968417838, 33.060396198677083}, {5.1922921581157908, 66.854301452103215}, {31.619281802149157, 25.269248720849514}, {81.541621071073038, 70.025341524754353}},
112{{46.078911165743556, 48.259962651999651}, {20.24450549867214, 49.403916182650214}, {0.26325131778756683, 24.46489805563581}, {15.915006546264051, 83.515023059917155}},
113
114{{65.454505973241524, 93.881892270353575}, {45.867360264932437, 92.723972719499827}, {2.1464054482739447, 74.636369140183717}, {33.774068594804994, 40.770872887582925}},
115{{72.963387832494163, 95.659300729473728}, {11.809496633619768, 82.209921247423594}, {13.456139067865974, 57.329313623406605}, {36.060621606214262, 70.867335643091849}},
116
117{{32.484981432782945, 75.082940782924624}, {42.467313093350882, 48.131159948246157}, {3.5963115764764657, 43.208665839959245}, {79.442476890721579, 89.709102357602262}},
118{{18.98573861410177, 93.308887208490106}, {40.405250173250792, 91.039661826118675}, {8.0467721950480584, 42.100282172719147}, {40.883324221187891, 26.030185504830527}},
119
120{{7.5374809128872498, 82.441702896003477}, {22.444346930107265, 22.138854312775123}, {66.76091829629658, 50.753805856571446}, {78.193478508942519, 97.7932997968948}},
121{{97.700573130371311, 53.53260215070685}, {87.72443481149358, 84.575876772671876}, {19.215031396232092, 47.032676472809484}, {11.989686410869325, 10.659507480757082}},
122
123{{26.192053931854691, 9.8504326817814416}, {10.174241480498686, 98.476562741434464}, {21.177712558385782, 33.814968789841501}, {75.329030899018534, 55.02231980442177}},
124{{56.222082700683771, 24.54395039218662}, {95.589995289030483, 81.050822735322086}, {28.180450866082897, 28.837706255185282}, {60.128952916771617, 87.311672180570511}},
125
126{{42.449716172390481, 52.379709366885805}, {27.896043159019225, 48.797373636065686}, {92.770268299044233, 89.899302036454571}, {12.102066544863426, 99.43241951960718}},
127{{45.77532924980639, 45.958701495993274}, {37.458701356062065, 68.393691335056758}, {37.569326692060258, 27.673713456687381}, {60.674866037757539, 62.47349659096146}},
128
129{{67.426548091427676, 37.993772624988935}, {23.483695892376684, 90.476863174921306}, {35.597065061143162, 79.872482633158796}, {75.38634169631932, 18.244890038969412}},
130{{61.336508189019057, 82.693132843213675}, {44.639380902349664, 54.074825790745592}, {16.815615499771951, 20.049704667203923}, {41.866884958868326, 56.735503699973002}},
131
132{{67.4265481, 37.9937726}, {23.4836959, 90.4768632}, {35.5970651, 79.8724826}, {75.3863417, 18.24489}},
133{{61.3365082, 82.6931328}, {44.6393809, 54.0748258}, {16.8156155, 20.0497047}, {41.866885, 56.7355037}},
134
135{{18.1312339, 31.6473732}, {95.5711034, 63.5350219}, {92.3283165, 62.0158945}, {18.5656052, 32.1268808}},
136{{97.402018, 35.7169972}, {33.1127443, 25.8935163}, {1.13970027, 54.9424981}, {56.4860195, 60.529264}},
137};
138
139const size_t testSetCount = sizeof(testSet) / sizeof(testSet[0]);
140
141static const Cubic newTestSet[] = {
142{{1,3}, {5,6}, {5,3}, {5,4}},
143{{3,5}, {4,5}, {3,1}, {6,5}},
144
145{{0,5}, {0,5}, {5,4}, {6,4}},
146{{4,5}, {4,6}, {5,0}, {5,0}},
147
148{{0,4}, {1,3}, {5,4}, {4,2}},
149{{4,5}, {2,4}, {4,0}, {3,1}},
150
151{{0,2}, {1,5}, {3,2}, {4,1}},
152{{2,3}, {1,4}, {2,0}, {5,1}},
153
154{{0,2}, {2,3}, {5,1}, {3,2}},
155{{1,5}, {2,3}, {2,0}, {3,2}},
156
157{{2,6}, {4,5}, {1,0}, {6,1}},
158{{0,1}, {1,6}, {6,2}, {5,4}},
159
160{{0,1}, {1,2}, {6,5}, {5,4}},
161{{5,6}, {4,5}, {1,0}, {2,1}},
162
163{{2.5119999999999996, 1.5710000000000002}, {2.6399999999999983, 1.6599999999999997}, {2.8000000000000007, 1.8000000000000003}, {3, 2}},
164{{2.4181876227114887, 1.9849772580462195}, {2.8269904869227211, 2.009330650246834}, {3.2004679292461624, 1.9942047174679169}, {3.4986199496818058, 2.0035994597094731}},
165
166{{2,3}, {1,4}, {1,0}, {6,0}},
167{{0,1}, {0,6}, {3,2}, {4,1}},
168
169{{0,2}, {1,5}, {1,0}, {6,1}},
170{{0,1}, {1,6}, {2,0}, {5,1}},
171
172{{0,1}, {1,5}, {2,1}, {4,0}},
173{{1,2}, {0,4}, {1,0}, {5,1}},
174
175{{0,1}, {3,5}, {2,1}, {3,1}},
176{{1,2}, {1,3}, {1,0}, {5,3}},
177
178{{0,1}, {2,5}, {6,0}, {5,3}},
179{{0,6}, {3,5}, {1,0}, {5,2}},
180
181{{0,1}, {3,6}, {1,0}, {5,2}},
182{{0,1}, {2,5}, {1,0}, {6,3}},
183
184{{1,2},{5,6},{1,0},{1,0}},
185{{0,1},{0,1},{2,1},{6,5}},
186
187{{0,6},{1,2},{1,0},{1,0}},
188{{0,1},{0,1},{6,0},{2,1}},
189
190{{0,2},{0,1},{3,0},{1,0}},
191{{0,3},{0,1},{2,0},{1,0}},
192};
193
194const size_t newTestSetCount = sizeof(newTestSet) / sizeof(newTestSet[0]);
195
196#if 0
197static void oneOff(const Cubic& cubic1, const Cubic& cubic2) {
198    SkTDArray<Quadratic> quads1;
199    cubic_to_quadratics(cubic1, calcPrecision(cubic1), quads1);
200#if SHOW_ORIGINAL
201    SkDebugf("computed quadratics given\n");
202    SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}}, {%1.9g,%1.9g}},\n",
203        cubic1[0].x, cubic1[0].y, cubic1[1].x, cubic1[1].y,
204        cubic1[2].x, cubic1[2].y, cubic1[3].x, cubic1[3].y));
205    SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}}, {%1.9g,%1.9g}},\n",
206        cubic2[0].x, cubic2[0].y, cubic2[1].x, cubic2[1].y,
207        cubic2[2].x, cubic2[2].y, cubic2[3].x, cubic2[3].y));
208#endif
209#if ONE_OFF_DEBUG
210    SkDebugf("computed quadratics set 1\n");
211    for (int index = 0; index < quads1.count(); ++index) {
212        const Quadratic& q = quads1[index];
213        SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
214                 q[1].x, q[1].y,  q[2].x, q[2].y);
215    }
216#endif
217    SkTDArray<Quadratic> quads2;
218    cubic_to_quadratics(cubic2, calcPrecision(cubic2), quads2);
219#if ONE_OFF_DEBUG
220    SkDebugf("computed quadratics set 2\n");
221    for (int index = 0; index < quads2.count(); ++index) {
222        const Quadratic& q = quads2[index];
223        SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
224                 q[1].x, q[1].y,  q[2].x, q[2].y);
225    }
226#endif
227    Intersections intersections2, intersections3;
228    intersect2(cubic1, cubic2, intersections2);
229    intersect3(cubic1, cubic2, intersections3);
230    int pt1, pt2, pt3;
231    bool found;
232    double tt1, tt2, last = -1;
233    _Point xy1, xy2;
234    for (pt1 = 0; pt1 < intersections2.used(); ++pt1) {
235        tt1 = intersections2.fT[0][pt1];
236        SkASSERT(!approximately_equal(last, tt1));
237        last = tt1;
238        xy_at_t(cubic1, tt1, xy1.x, xy1.y);
239        pt2 = intersections2.fFlip ? intersections2.used() - pt1 - 1 : pt1;
240        tt2 = intersections2.fT[1][pt2];
241        xy_at_t(cubic2, tt2, xy2.x, xy2.y);
242#if ONE_OFF_DEBUG
243        SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n",
244                __FUNCTION__, tt1, xy1.x, xy1.y, intersections2.fPt[pt1].x,
245                intersections2.fPt[pt1].y, xy2.x, xy2.y, tt2);
246#endif
247        SkASSERT(xy1.approximatelyEqual(xy2));
248#ifdef SK_DEBUG
249        found = false;
250        for (pt3 = 0; pt3 < intersections3.used(); ++pt3) {
251            if (roughly_equal(tt1, intersections3.fT[0][pt3])) {
252                found = true;
253                break;
254            }
255        }
256        SkASSERT(found);
257#endif
258    }
259    last = -1;
260    for (pt3 = 0; pt3 < intersections3.used(); ++pt3) {
261        found = false;
262        double tt3 = intersections3.fT[0][pt3];
263        SkASSERT(!approximately_equal(last, tt3));
264        last = tt3;
265        for (pt1 = 0; pt1 < intersections2.used(); ++pt1) {
266            if (approximately_equal(tt3, intersections2.fT[0][pt1])) {
267                found = true;
268                break;
269            }
270        }
271        if (!found) {
272            tt1 = intersections3.fT[0][pt3];
273            xy_at_t(cubic1, tt1, xy1.x, xy1.y);
274            pt2 = intersections3.fFlip ? intersections3.used() - pt3 - 1 : pt3;
275            tt2 = intersections3.fT[1][pt2];
276            xy_at_t(cubic2, tt2, xy2.x, xy2.y);
277    #if ONE_OFF_DEBUG
278            SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n",
279                    __FUNCTION__, tt1, xy1.x, xy1.y, intersections3.fPt[pt1].x,
280                    intersections3.fPt[pt1].y, xy2.x, xy2.y, tt2);
281    #endif
282            SkASSERT(xy1.approximatelyEqual(xy2));
283            SkDebugf("%s missing in intersect2\n", __FUNCTION__);
284        }
285    }
286}
287#endif
288
289static void oneOff3(const Cubic& cubic1, const Cubic& cubic2) {
290#if ONE_OFF_DEBUG
291    SkDebugf("computed quadratics given\n");
292    SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n",
293        cubic1[0].x, cubic1[0].y, cubic1[1].x, cubic1[1].y,
294        cubic1[2].x, cubic1[2].y, cubic1[3].x, cubic1[3].y);
295    SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n",
296        cubic2[0].x, cubic2[0].y, cubic2[1].x, cubic2[1].y,
297        cubic2[2].x, cubic2[2].y, cubic2[3].x, cubic2[3].y);
298#endif
299    SkTDArray<Quadratic> quads1;
300    cubic_to_quadratics(cubic1, calcPrecision(cubic1), quads1);
301#if ONE_OFF_DEBUG
302    SkDebugf("computed quadratics set 1\n");
303    for (int index = 0; index < quads1.count(); ++index) {
304        const Quadratic& q = quads1[index];
305        SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
306                 q[1].x, q[1].y,  q[2].x, q[2].y);
307    }
308#endif
309    SkTDArray<Quadratic> quads2;
310    cubic_to_quadratics(cubic2, calcPrecision(cubic2), quads2);
311#if ONE_OFF_DEBUG
312    SkDebugf("computed quadratics set 2\n");
313    for (int index = 0; index < quads2.count(); ++index) {
314        const Quadratic& q = quads2[index];
315        SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
316                 q[1].x, q[1].y,  q[2].x, q[2].y);
317    }
318#endif
319    Intersections intersections3;
320    intersect3(cubic1, cubic2, intersections3);
321    int pt2, pt3;
322    double tt1, tt2, last = -1;
323    _Point xy1, xy2;
324    for (pt3 = 0; pt3 < intersections3.used(); ++pt3) {
325        double tt3 = intersections3.fT[0][pt3];
326     //   SkASSERT(!approximately_equal(last, tt3));
327        last = tt3;
328        tt1 = intersections3.fT[0][pt3];
329        xy_at_t(cubic1, tt1, xy1.x, xy1.y);
330        pt2 = intersections3.fFlip ? intersections3.used() - pt3 - 1 : pt3;
331        tt2 = intersections3.fT[1][pt2];
332        xy_at_t(cubic2, tt2, xy2.x, xy2.y);
333#if ONE_OFF_DEBUG
334        SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n",
335                __FUNCTION__, tt1, xy1.x, xy1.y, intersections3.fPt[pt3].x,
336                intersections3.fPt[pt3].y, xy2.x, xy2.y, tt2);
337#endif
338        SkASSERT(xy1.approximatelyEqual(xy2));
339    }
340}
341
342#if 0
343static int fails[][2] = {   {0, 23}, // fails in intersect2 recursing
344                            {2, 7},  // answers differ, but neither is correct ('3' is closer)
345                            {3, 26}, // fails in intersect2 recursing
346                            {4, 9},  // fails in intersect2 recursing
347                            {4, 10}, // fails in intersect2 recursing
348                            {10, 17}, // fails in intersect2 recursing
349                            {12, 14}, // loops indefinitely
350                            {12, 21}, // fails in intersect2 recursing
351                            {13, 21}, // fails in intersect2 recursing
352                            {14, 21}, // fails in intersect2 recursing
353                            {17, 25}, // fails in intersect2 recursing
354                            {23, 25}, // fails in intersect2 recursing
355};
356
357static int failCount = sizeof(fails) / sizeof(fails[0]);
358#endif
359
360static void oneOff(int outer, int inner) {
361    const Cubic& cubic1 = testSet[outer];
362    const Cubic& cubic2 = testSet[inner];
363#if 0
364    bool failing = false;
365    for (int i = 0; i < failCount; ++i) {
366        if ((fails[i][0] == outer && fails[i][1] == inner)
367                || (fails[i][1] == outer && fails[i][0] == inner)) {
368            failing = true;
369            break;
370        }
371    }
372    if (!failing) {
373        oneOff(cubic1, cubic2);
374    } else {
375#endif
376        oneOff3(cubic1, cubic2);
377//    }
378}
379
380void CubicIntersection_OneOffTest() {
381    oneOff(0, 1);
382}
383
384static void newOneOff(int outer, int inner) {
385    const Cubic& cubic1 = newTestSet[outer];
386    const Cubic& cubic2 = newTestSet[inner];
387    oneOff3(cubic1, cubic2);
388}
389
390void CubicIntersection_NewOneOffTest() {
391    newOneOff(0, 1);
392}
393
394static void oneOffTests() {
395    for (size_t outer = 0; outer < testSetCount - 1; ++outer) {
396        for (size_t inner = outer + 1; inner < testSetCount; ++inner) {
397            oneOff(outer, inner);
398        }
399    }
400}
401
402void CubicIntersection_OneOffTests() {
403    oneOffTests();
404}
405
406#define DEBUG_CRASH 0
407
408class CubicChopper {
409public:
410
411// only finds one intersection
412CubicChopper(const Cubic& c1, const Cubic& c2)
413    : cubic1(c1)
414    , cubic2(c2)
415    , depth(0) {
416}
417
418bool intersect(double minT1, double maxT1, double minT2, double maxT2) {
419    Cubic sub1, sub2;
420    // FIXME: carry last subdivide and reduceOrder result with cubic
421    sub_divide(cubic1, minT1, maxT1, sub1);
422    sub_divide(cubic2, minT2, maxT2, sub2);
423    Intersections i;
424    intersect3(sub1, sub2, i);
425    if (i.used() == 0) {
426        return false;
427    }
428    double x1, y1, x2, y2;
429    t1 = minT1 + i.fT[0][0] * (maxT1 - minT1);
430    t2 = minT2 + i.fT[1][0] * (maxT2 - minT2);
431    xy_at_t(cubic1, t1, x1, y1);
432    xy_at_t(cubic2, t2, x2, y2);
433    if (AlmostEqualUlps(x1, x2) && AlmostEqualUlps(y1, y2)) {
434        return true;
435    }
436    double half1 = (minT1 + maxT1) / 2;
437    double half2 = (minT2 + maxT2) / 2;
438    ++depth;
439    bool result;
440    if (depth & 1) {
441        result = intersect(minT1, half1, minT2, maxT2) || intersect(half1, maxT1, minT2, maxT2)
442            || intersect(minT1, maxT1, minT2, half2) || intersect(minT1, maxT1, half2, maxT2);
443    } else {
444        result = intersect(minT1, maxT1, minT2, half2) || intersect(minT1, maxT1, half2, maxT2)
445            || intersect(minT1, half1, minT2, maxT2) || intersect(half1, maxT1, minT2, maxT2);
446    }
447    --depth;
448    return result;
449}
450
451const Cubic& cubic1;
452const Cubic& cubic2;
453double t1;
454double t2;
455int depth;
456};
457
458#define TRY_OLD 0 // old way fails on test == 1
459
460void CubicIntersection_RandTestOld() {
461    srand(0);
462    const int tests = 1000000; // 10000000;
463    double largestFactor = DBL_MAX;
464    for (int test = 0; test < tests; ++test) {
465        Cubic cubic1, cubic2;
466        for (int i = 0; i < 4; ++i) {
467            cubic1[i].x = (double) rand() / RAND_MAX * 100;
468            cubic1[i].y = (double) rand() / RAND_MAX * 100;
469            cubic2[i].x = (double) rand() / RAND_MAX * 100;
470            cubic2[i].y = (double) rand() / RAND_MAX * 100;
471        }
472        if (test == 2513) { // the pair crosses three times, but the quadratic approximation
473            continue; // only sees one -- should be OK to ignore the other two?
474        }
475        if (test == 12932) { // this exposes a weakness when one cubic touches the other but
476            continue; // does not touch the quad approximation. Captured in qc.htm as cubic15
477        }
478    #if DEBUG_CRASH
479        char str[1024];
480        sprintf(str, "{{%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}},\n"
481            "{{%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}},\n",
482                cubic1[0].x, cubic1[0].y,  cubic1[1].x, cubic1[1].y, cubic1[2].x, cubic1[2].y,
483                cubic1[3].x, cubic1[3].y,
484                cubic2[0].x, cubic2[0].y,  cubic2[1].x, cubic2[1].y, cubic2[2].x, cubic2[2].y,
485                cubic2[3].x, cubic2[3].y);
486    #endif
487        _Rect rect1, rect2;
488        rect1.setBounds(cubic1);
489        rect2.setBounds(cubic2);
490        bool boundsIntersect = rect1.left <= rect2.right && rect2.left <= rect2.right
491                && rect1.top <= rect2.bottom && rect2.top <= rect1.bottom;
492        Intersections i1, i2;
493    #if TRY_OLD
494        bool oldIntersects = intersect(cubic1, cubic2, i1);
495    #else
496        bool oldIntersects = false;
497    #endif
498        if (test == -1) {
499            SkDebugf("ready...\n");
500        }
501        bool newIntersects = intersect3(cubic1, cubic2, i2);
502        if (!boundsIntersect && (oldIntersects || newIntersects)) {
503    #if DEBUG_CRASH
504            SkDebugf("%s %d unexpected intersection boundsIntersect=%d oldIntersects=%d"
505                    " newIntersects=%d\n%s %s\n", __FUNCTION__, test, boundsIntersect,
506                    oldIntersects, newIntersects, __FUNCTION__, str);
507    #endif
508            SkASSERT(0);
509        }
510        if (oldIntersects && !newIntersects) {
511    #if DEBUG_CRASH
512            SkDebugf("%s %d missing intersection oldIntersects=%d newIntersects=%d\n%s %s\n",
513                    __FUNCTION__, test, oldIntersects, newIntersects, __FUNCTION__, str);
514    #endif
515            SkASSERT(0);
516        }
517        if (!oldIntersects && !newIntersects) {
518            continue;
519        }
520        if (i2.used() > 1) {
521            continue;
522            // just look at single intercepts for simplicity
523        }
524        Intersections self1, self2; // self-intersect checks
525        if (intersect(cubic1, self1)) {
526            continue;
527        }
528        if (intersect(cubic2, self2)) {
529            continue;
530        }
531        // binary search for range necessary to enclose real intersection
532        CubicChopper c(cubic1, cubic2);
533        bool result = c.intersect(0, 1, 0, 1);
534        if (!result) {
535            // FIXME: a failure here probably means that a core routine used by CubicChopper is failing
536            continue;
537        }
538        double delta1 = fabs(c.t1 - i2.fT[0][0]);
539        double delta2 = fabs(c.t2 - i2.fT[1][0]);
540        double calc1 = calcPrecision(cubic1);
541        double calc2 = calcPrecision(cubic2);
542        double factor1 = calc1 / delta1;
543        double factor2 = calc2 / delta2;
544        SkDebugf("%s %d calc1=%1.9g delta1=%1.9g factor1=%1.9g calc2=%1.9g delta2=%1.9g"
545                " factor2=%1.9g\n", __FUNCTION__, test,
546                calc1, delta1, factor1, calc2, delta2, factor2);
547        if (factor1 < largestFactor) {
548            SkDebugf("WE HAVE A WINNER! %1.9g\n", factor1);
549    #if DEBUG_CRASH
550            SkDebugf("%s\n", str);
551    #endif
552            oneOff3(cubic1, cubic2);
553            largestFactor = factor1;
554        }
555        if (factor2 < largestFactor) {
556            SkDebugf("WE HAVE A WINNER! %1.9g\n", factor2);
557    #if DEBUG_CRASH
558            SkDebugf("%s\n", str);
559    #endif
560            oneOff3(cubic1, cubic2);
561            largestFactor = factor2;
562        }
563    }
564}
565
566void CubicIntersection_RandTest() {
567    srand(0);
568    const int tests = 10000000;
569    for (int test = 0; test < tests; ++test) {
570        Cubic cubic1, cubic2;
571        for (int i = 0; i < 4; ++i) {
572            cubic1[i].x = (double) rand() / RAND_MAX * 100;
573            cubic1[i].y = (double) rand() / RAND_MAX * 100;
574            cubic2[i].x = (double) rand() / RAND_MAX * 100;
575            cubic2[i].y = (double) rand() / RAND_MAX * 100;
576        }
577    #if DEBUG_CRASH
578        char str[1024];
579        sprintf(str, "{{%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}},\n"
580            "{{%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}},\n",
581                cubic1[0].x, cubic1[0].y,  cubic1[1].x, cubic1[1].y, cubic1[2].x, cubic1[2].y,
582                cubic1[3].x, cubic1[3].y,
583                cubic2[0].x, cubic2[0].y,  cubic2[1].x, cubic2[1].y, cubic2[2].x, cubic2[2].y,
584                cubic2[3].x, cubic2[3].y);
585    #endif
586        _Rect rect1, rect2;
587        rect1.setBounds(cubic1);
588        rect2.setBounds(cubic2);
589        bool boundsIntersect = rect1.left <= rect2.right && rect2.left <= rect2.right
590                && rect1.top <= rect2.bottom && rect2.top <= rect1.bottom;
591        if (test == -1) {
592            SkDebugf("ready...\n");
593        }
594        Intersections intersections2;
595        bool newIntersects = intersect3(cubic1, cubic2, intersections2);
596        if (!boundsIntersect && newIntersects) {
597    #if DEBUG_CRASH
598            SkDebugf("%s %d unexpected intersection boundsIntersect=%d "
599                    " newIntersects=%d\n%s %s\n", __FUNCTION__, test, boundsIntersect,
600                    newIntersects, __FUNCTION__, str);
601    #endif
602            SkASSERT(0);
603        }
604        for (int pt = 0; pt < intersections2.used(); ++pt) {
605            double tt1 = intersections2.fT[0][pt];
606            _Point xy1, xy2;
607            xy_at_t(cubic1, tt1, xy1.x, xy1.y);
608            int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
609            double tt2 = intersections2.fT[1][pt2];
610            xy_at_t(cubic2, tt2, xy2.x, xy2.y);
611        #if 0
612            SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
613                tt1, xy1.x, xy1.y, xy2.x, xy2.y, tt2);
614        #endif
615            SkASSERT(xy1.approximatelyEqual(xy2));
616        }
617    }
618}
619
620static void intersectionFinder(int index0, int index1, double t1Seed, double t2Seed,
621        double t1Step, double t2Step) {
622    const Cubic& cubic1 = newTestSet[index0];
623    const Cubic& cubic2 = newTestSet[index1];
624    _Point t1[3], t2[3];
625    bool toggle = true;
626    do {
627        xy_at_t(cubic1, t1Seed - t1Step, t1[0].x, t1[0].y);
628        xy_at_t(cubic1, t1Seed,          t1[1].x, t1[1].y);
629        xy_at_t(cubic1, t1Seed + t1Step, t1[2].x, t1[2].y);
630        xy_at_t(cubic2, t2Seed - t2Step, t2[0].x, t2[0].y);
631        xy_at_t(cubic2, t2Seed,          t2[1].x, t2[1].y);
632        xy_at_t(cubic2, t2Seed + t2Step, t2[2].x, t2[2].y);
633        double dist[3][3];
634        dist[1][1] = t1[1].distance(t2[1]);
635        int best_i = 1, best_j = 1;
636        for (int i = 0; i < 3; ++i) {
637            for (int j = 0; j < 3; ++j) {
638                if (i == 1 && j == 1) {
639                    continue;
640                }
641                dist[i][j] = t1[i].distance(t2[j]);
642                if (dist[best_i][best_j] > dist[i][j]) {
643                    best_i = i;
644                    best_j = j;
645                }
646            }
647        }
648        if (best_i == 0) {
649            t1Seed -= t1Step;
650        } else if (best_i == 2) {
651            t1Seed += t1Step;
652        }
653        if (best_j == 0) {
654            t2Seed -= t2Step;
655        } else if (best_j == 2) {
656            t2Seed += t2Step;
657        }
658        if (best_i == 1 && best_j == 1) {
659            if ((toggle ^= true)) {
660                t1Step /= 2;
661            } else {
662                t2Step /= 2;
663            }
664        }
665    } while (!t1[1].approximatelyEqual(t2[1]));
666    t1Step = t2Step = 0.1;
667    double t10 = t1Seed - t1Step * 2;
668    double t12 = t1Seed + t1Step * 2;
669    double t20 = t2Seed - t2Step * 2;
670    double t22 = t2Seed + t2Step * 2;
671    _Point test;
672    while (!approximately_zero(t1Step)) {
673        xy_at_t(cubic1, t10, test.x, test.y);
674        t10 += t1[1].approximatelyEqual(test) ? -t1Step : t1Step;
675        t1Step /= 2;
676    }
677    t1Step = 0.1;
678    while (!approximately_zero(t1Step)) {
679        xy_at_t(cubic1, t12, test.x, test.y);
680        t12 -= t1[1].approximatelyEqual(test) ? -t1Step : t1Step;
681        t1Step /= 2;
682    }
683    while (!approximately_zero(t2Step)) {
684        xy_at_t(cubic2, t20, test.x, test.y);
685        t20 += t2[1].approximatelyEqual(test) ? -t2Step : t2Step;
686        t2Step /= 2;
687    }
688    t2Step = 0.1;
689    while (!approximately_zero(t2Step)) {
690        xy_at_t(cubic2, t22, test.x, test.y);
691        t22 -= t2[1].approximatelyEqual(test) ? -t2Step : t2Step;
692        t2Step /= 2;
693    }
694#if ONE_OFF_DEBUG
695    SkDebugf("%s t1=(%1.9g<%1.9g<%1.9g) t2=(%1.9g<%1.9g<%1.9g)\n", __FUNCTION__,
696        t10, t1Seed, t12, t20, t2Seed, t22);
697    _Point p10 = xy_at_t(cubic1, t10);
698    _Point p1Seed = xy_at_t(cubic1, t1Seed);
699    _Point p12 = xy_at_t(cubic1, t12);
700    SkDebugf("%s p1=(%1.9g,%1.9g)<(%1.9g,%1.9g)<(%1.9g,%1.9g)\n", __FUNCTION__,
701        p10.x, p10.y, p1Seed.x, p1Seed.y, p12.x, p12.y);
702    _Point p20 = xy_at_t(cubic2, t20);
703    _Point p2Seed = xy_at_t(cubic2, t2Seed);
704    _Point p22 = xy_at_t(cubic2, t22);
705    SkDebugf("%s p2=(%1.9g,%1.9g)<(%1.9g,%1.9g)<(%1.9g,%1.9g)\n", __FUNCTION__,
706        p20.x, p20.y, p2Seed.x, p2Seed.y, p22.x, p22.y);
707#endif
708}
709
710void CubicIntersection_IntersectionFinder() {
711
712 //   double t1Seed = 0.87;
713 //   double t2Seed = 0.87;
714    double t1Step = 0.000001;
715    double t2Step = 0.000001;
716    intersectionFinder(0, 1, 0.855895664, 0.864850875, t1Step, t2Step);
717    intersectionFinder(0, 1, 0.865207906, 0.865207887, t1Step, t2Step);
718    intersectionFinder(0, 1, 0.865213351, 0.865208087, t1Step, t2Step);
719}
720
721static void coincidentTest() {
722#if 0
723    Cubic cubic1 = {{0, 1}, {0, 2}, {1, 0}, {1, 0}};
724    Cubic cubic2 = {{0, 1}, {0, 2}, {1, 0}, {6, 1}};
725#endif
726}
727
728void CubicIntersection_SelfTest() {
729    const Cubic selfSet[] = {
730        {{0,2}, {2,3}, {5,1}, {3,2}},
731        {{0,2}, {3,5}, {5,0}, {4,2}},
732        {{3.34,8.98}, {1.95,10.27}, {3.76,7.65}, {4.96,10.64}},
733        {{3.13,2.74}, {1.08,4.62}, {3.71,0.94}, {2.01,3.81}},
734        {{6.71,3.14}, {7.99,2.75}, {8.27,1.96}, {6.35,3.57}},
735        {{12.81,7.27}, {7.22,6.98}, {12.49,8.97}, {11.42,6.18}},
736    };
737    size_t selfSetCount = sizeof(selfSet) / sizeof(selfSet[0]);
738    size_t firstFail = 1;
739    for (size_t index = firstFail; index < selfSetCount; ++index) {
740        const Cubic& cubic = selfSet[index];
741    #if ONE_OFF_DEBUG
742        int idx2;
743        double max[3];
744        int ts = find_cubic_max_curvature(cubic, max);
745        for (idx2 = 0; idx2 < ts; ++idx2) {
746            SkDebugf("%s max[%d]=%1.9g (%1.9g, %1.9g)\n", __FUNCTION__, idx2,
747                    max[idx2], xy_at_t(cubic, max[idx2]).x, xy_at_t(cubic, max[idx2]).y);
748        }
749        SkTDArray<double> ts1;
750        SkTDArray<Quadratic> quads1;
751        cubic_to_quadratics(cubic, calcPrecision(cubic), ts1);
752        for (idx2 = 0; idx2 < ts1.count(); ++idx2) {
753            SkDebugf("%s t[%d]=%1.9g\n", __FUNCTION__, idx2, ts1[idx2]);
754        }
755        cubic_to_quadratics(cubic, calcPrecision(cubic), quads1);
756        for (idx2 = 0; idx2 < quads1.count(); ++idx2) {
757            const Quadratic& q = quads1[idx2];
758            SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n",
759                    q[0].x, q[0].y,  q[1].x, q[1].y,  q[2].x, q[2].y);
760        }
761        SkDebugf("\n");
762    #endif
763        Intersections i;
764        SkDEBUGCODE(int result = ) intersect(cubic, i);
765        SkASSERT(result == 1);
766        SkASSERT(i.used() == 1);
767        SkASSERT(!approximately_equal(i.fT[0][0], i.fT[1][0]));
768        _Point pt1 = xy_at_t(cubic, i.fT[0][0]);
769        _Point pt2 = xy_at_t(cubic, i.fT[1][0]);
770        SkASSERT(pt1.approximatelyEqual(pt2));
771    }
772}
773
774void CubicIntersection_Test() {
775    oneOffTests();
776    coincidentTest();
777    standardTestCases();
778}
779