1// Another approach is to start with the implicit form of one curve and solve
2// (seek implicit coefficients in QuadraticParameter.cpp
3// by substituting in the parametric form of the other.
4// The downside of this approach is that early rejects are difficult to come by.
5// http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step
6
7#include "SkDQuadImplicit.h"
8#include "SkIntersections.h"
9#include "SkPathOpsLine.h"
10#include "SkQuarticRoot.h"
11#include "SkTArray.h"
12#include "SkTSort.h"
13
14/* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F
15 * and given x = at^2 + bt + c  (the parameterized form)
16 *           y = dt^2 + et + f
17 * then
18 * 0 = A(at^2+bt+c)(at^2+bt+c)+B(at^2+bt+c)(dt^2+et+f)+C(dt^2+et+f)(dt^2+et+f)+D(at^2+bt+c)+E(dt^2+et+f)+F
19 */
20
21static int findRoots(const SkDQuadImplicit& i, const SkDQuad& quad, double roots[4],
22        bool oneHint, bool flip, int firstCubicRoot) {
23    SkDQuad flipped;
24    const SkDQuad& q = flip ? (flipped = quad.flip()) : quad;
25    double a, b, c;
26    SkDQuad::SetABC(&q[0].fX, &a, &b, &c);
27    double d, e, f;
28    SkDQuad::SetABC(&q[0].fY, &d, &e, &f);
29    const double t4 =     i.x2() *  a * a
30                    +     i.xy() *  a * d
31                    +     i.y2() *  d * d;
32    const double t3 = 2 * i.x2() *  a * b
33                    +     i.xy() * (a * e +     b * d)
34                    + 2 * i.y2() *  d * e;
35    const double t2 =     i.x2() * (b * b + 2 * a * c)
36                    +     i.xy() * (c * d +     b * e + a * f)
37                    +     i.y2() * (e * e + 2 * d * f)
38                    +     i.x()  *  a
39                    +     i.y()  *  d;
40    const double t1 = 2 * i.x2() *  b * c
41                    +     i.xy() * (c * e + b * f)
42                    + 2 * i.y2() *  e * f
43                    +     i.x()  *  b
44                    +     i.y()  *  e;
45    const double t0 =     i.x2() *  c * c
46                    +     i.xy() *  c * f
47                    +     i.y2() *  f * f
48                    +     i.x()  *  c
49                    +     i.y()  *  f
50                    +     i.c();
51    int rootCount = SkReducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots);
52    if (rootCount < 0) {
53        rootCount = SkQuarticRootsReal(firstCubicRoot, t4, t3, t2, t1, t0, roots);
54    }
55    if (flip) {
56        for (int index = 0; index < rootCount; ++index) {
57            roots[index] = 1 - roots[index];
58        }
59    }
60    return rootCount;
61}
62
63static int addValidRoots(const double roots[4], const int count, double valid[4]) {
64    int result = 0;
65    int index;
66    for (index = 0; index < count; ++index) {
67        if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) {
68            continue;
69        }
70        double t = 1 - roots[index];
71        if (approximately_less_than_zero(t)) {
72            t = 0;
73        } else if (approximately_greater_than_one(t)) {
74            t = 1;
75        }
76        valid[result++] = t;
77    }
78    return result;
79}
80
81static bool only_end_pts_in_common(const SkDQuad& q1, const SkDQuad& q2) {
82// the idea here is to see at minimum do a quick reject by rotating all points
83// to either side of the line formed by connecting the endpoints
84// if the opposite curves points are on the line or on the other side, the
85// curves at most intersect at the endpoints
86    for (int oddMan = 0; oddMan < 3; ++oddMan) {
87        const SkDPoint* endPt[2];
88        for (int opp = 1; opp < 3; ++opp) {
89            int end = oddMan ^ opp;  // choose a value not equal to oddMan
90            if (3 == end) {  // and correct so that largest value is 1 or 2
91                end = opp;
92            }
93            endPt[opp - 1] = &q1[end];
94        }
95        double origX = endPt[0]->fX;
96        double origY = endPt[0]->fY;
97        double adj = endPt[1]->fX - origX;
98        double opp = endPt[1]->fY - origY;
99        double sign = (q1[oddMan].fY - origY) * adj - (q1[oddMan].fX - origX) * opp;
100        if (approximately_zero(sign)) {
101            goto tryNextHalfPlane;
102        }
103        for (int n = 0; n < 3; ++n) {
104            double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp;
105            if (test * sign > 0 && !precisely_zero(test)) {
106                goto tryNextHalfPlane;
107            }
108        }
109        return true;
110tryNextHalfPlane:
111        ;
112    }
113    return false;
114}
115
116// returns false if there's more than one intercept or the intercept doesn't match the point
117// returns true if the intercept was successfully added or if the
118// original quads need to be subdivided
119static bool add_intercept(const SkDQuad& q1, const SkDQuad& q2, double tMin, double tMax,
120                          SkIntersections* i, bool* subDivide) {
121    double tMid = (tMin + tMax) / 2;
122    SkDPoint mid = q2.ptAtT(tMid);
123    SkDLine line;
124    line[0] = line[1] = mid;
125    SkDVector dxdy = q2.dxdyAtT(tMid);
126    line[0] -= dxdy;
127    line[1] += dxdy;
128    SkIntersections rootTs;
129    rootTs.allowNear(false);
130    int roots = rootTs.intersect(q1, line);
131    if (roots == 0) {
132        if (subDivide) {
133            *subDivide = true;
134        }
135        return true;
136    }
137    if (roots == 2) {
138        return false;
139    }
140    SkDPoint pt2 = q1.ptAtT(rootTs[0][0]);
141    if (!pt2.approximatelyEqual(mid)) {
142        return false;
143    }
144    i->insertSwap(rootTs[0][0], tMid, pt2);
145    return true;
146}
147
148static bool is_linear_inner(const SkDQuad& q1, double t1s, double t1e, const SkDQuad& q2,
149                            double t2s, double t2e, SkIntersections* i, bool* subDivide) {
150    SkDQuad hull = q1.subDivide(t1s, t1e);
151    SkDLine line = {{hull[2], hull[0]}};
152    const SkDLine* testLines[] = { &line, (const SkDLine*) &hull[0], (const SkDLine*) &hull[1] };
153    const size_t kTestCount = SK_ARRAY_COUNT(testLines);
154    SkSTArray<kTestCount * 2, double, true> tsFound;
155    for (size_t index = 0; index < kTestCount; ++index) {
156        SkIntersections rootTs;
157        rootTs.allowNear(false);
158        int roots = rootTs.intersect(q2, *testLines[index]);
159        for (int idx2 = 0; idx2 < roots; ++idx2) {
160            double t = rootTs[0][idx2];
161#if 0 // def SK_DEBUG   // FIXME : accurate for error = 16, error of 17.5 seen
162// {{{136.08723965397621, 1648.2814535211637}, {593.49031197259478, 1190.8784277439891}, {593.49031197259478, 544.0128173828125}}}
163// {{{-968.181396484375, 544.0128173828125}, {592.2825927734375, 870.552490234375}, {593.435302734375, 557.8828125}}}
164
165            SkDPoint qPt = q2.ptAtT(t);
166            SkDPoint lPt = testLines[index]->ptAtT(rootTs[1][idx2]);
167            SkASSERT(qPt.approximatelyDEqual(lPt));
168#endif
169            if (approximately_negative(t - t2s) || approximately_positive(t - t2e)) {
170                continue;
171            }
172            tsFound.push_back(rootTs[0][idx2]);
173        }
174    }
175    int tCount = tsFound.count();
176    if (tCount <= 0) {
177        return true;
178    }
179    double tMin, tMax;
180    if (tCount == 1) {
181        tMin = tMax = tsFound[0];
182    } else {
183        SkASSERT(tCount > 1);
184        SkTQSort<double>(tsFound.begin(), tsFound.end() - 1);
185        tMin = tsFound[0];
186        tMax = tsFound[tsFound.count() - 1];
187    }
188    SkDPoint end = q2.ptAtT(t2s);
189    bool startInTriangle = hull.pointInHull(end);
190    if (startInTriangle) {
191        tMin = t2s;
192    }
193    end = q2.ptAtT(t2e);
194    bool endInTriangle = hull.pointInHull(end);
195    if (endInTriangle) {
196        tMax = t2e;
197    }
198    int split = 0;
199    SkDVector dxy1, dxy2;
200    if (tMin != tMax || tCount > 2) {
201        dxy2 = q2.dxdyAtT(tMin);
202        for (int index = 1; index < tCount; ++index) {
203            dxy1 = dxy2;
204            dxy2 = q2.dxdyAtT(tsFound[index]);
205            double dot = dxy1.dot(dxy2);
206            if (dot < 0) {
207                split = index - 1;
208                break;
209            }
210        }
211    }
212    if (split == 0) {  // there's one point
213        if (add_intercept(q1, q2, tMin, tMax, i, subDivide)) {
214            return true;
215        }
216        i->swap();
217        return is_linear_inner(q2, tMin, tMax, q1, t1s, t1e, i, subDivide);
218    }
219    // At this point, we have two ranges of t values -- treat each separately at the split
220    bool result;
221    if (add_intercept(q1, q2, tMin, tsFound[split - 1], i, subDivide)) {
222        result = true;
223    } else {
224        i->swap();
225        result = is_linear_inner(q2, tMin, tsFound[split - 1], q1, t1s, t1e, i, subDivide);
226    }
227    if (add_intercept(q1, q2, tsFound[split], tMax, i, subDivide)) {
228        result = true;
229    } else {
230        i->swap();
231        result |= is_linear_inner(q2, tsFound[split], tMax, q1, t1s, t1e, i, subDivide);
232    }
233    return result;
234}
235
236static double flat_measure(const SkDQuad& q) {
237    SkDVector mid = q[1] - q[0];
238    SkDVector dxy = q[2] - q[0];
239    double length = dxy.length();  // OPTIMIZE: get rid of sqrt
240    return fabs(mid.cross(dxy) / length);
241}
242
243// FIXME ? should this measure both and then use the quad that is the flattest as the line?
244static bool is_linear(const SkDQuad& q1, const SkDQuad& q2, SkIntersections* i) {
245    double measure = flat_measure(q1);
246    // OPTIMIZE: (get rid of sqrt) use approximately_zero
247    if (!approximately_zero_sqrt(measure)) {
248        return false;
249    }
250    return is_linear_inner(q1, 0, 1, q2, 0, 1, i, NULL);
251}
252
253// FIXME: if flat measure is sufficiently large, then probably the quartic solution failed
254// avoid imprecision incurred with chopAt
255static void relaxed_is_linear(const SkDQuad* q1, double s1, double e1, const SkDQuad* q2,
256        double s2, double e2, SkIntersections* i) {
257    double m1 = flat_measure(*q1);
258    double m2 = flat_measure(*q2);
259    i->reset();
260    const SkDQuad* rounder, *flatter;
261    double sf, midf, ef, sr, er;
262    if (m2 < m1) {
263        rounder = q1;
264        sr = s1;
265        er = e1;
266        flatter = q2;
267        sf = s2;
268        midf = (s2 + e2) / 2;
269        ef = e2;
270    } else {
271        rounder = q2;
272        sr = s2;
273        er = e2;
274        flatter = q1;
275        sf = s1;
276        midf = (s1 + e1) / 2;
277        ef = e1;
278    }
279    bool subDivide = false;
280    is_linear_inner(*flatter, sf, ef, *rounder, sr, er, i, &subDivide);
281    if (subDivide) {
282        relaxed_is_linear(flatter, sf, midf, rounder, sr, er, i);
283        relaxed_is_linear(flatter, midf, ef, rounder, sr, er, i);
284    }
285    if (m2 < m1) {
286        i->swapPts();
287    }
288}
289
290// each time through the loop, this computes values it had from the last loop
291// if i == j == 1, the center values are still good
292// otherwise, for i != 1 or j != 1, four of the values are still good
293// and if i == 1 ^ j == 1, an additional value is good
294static bool binary_search(const SkDQuad& quad1, const SkDQuad& quad2, double* t1Seed,
295                          double* t2Seed, SkDPoint* pt) {
296    double tStep = ROUGH_EPSILON;
297    SkDPoint t1[3], t2[3];
298    int calcMask = ~0;
299    do {
300        if (calcMask & (1 << 1)) t1[1] = quad1.ptAtT(*t1Seed);
301        if (calcMask & (1 << 4)) t2[1] = quad2.ptAtT(*t2Seed);
302        if (t1[1].approximatelyEqual(t2[1])) {
303            *pt = t1[1];
304    #if ONE_OFF_DEBUG
305            SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __FUNCTION__,
306                    t1Seed, t2Seed, t1[1].fX, t1[1].fY, t2[1].fX, t2[1].fY);
307    #endif
308            return true;
309        }
310        if (calcMask & (1 << 0)) t1[0] = quad1.ptAtT(SkTMax(0., *t1Seed - tStep));
311        if (calcMask & (1 << 2)) t1[2] = quad1.ptAtT(SkTMin(1., *t1Seed + tStep));
312        if (calcMask & (1 << 3)) t2[0] = quad2.ptAtT(SkTMax(0., *t2Seed - tStep));
313        if (calcMask & (1 << 5)) t2[2] = quad2.ptAtT(SkTMin(1., *t2Seed + tStep));
314        double dist[3][3];
315        // OPTIMIZE: using calcMask value permits skipping some distance calcuations
316        //   if prior loop's results are moved to correct slot for reuse
317        dist[1][1] = t1[1].distanceSquared(t2[1]);
318        int best_i = 1, best_j = 1;
319        for (int i = 0; i < 3; ++i) {
320            for (int j = 0; j < 3; ++j) {
321                if (i == 1 && j == 1) {
322                    continue;
323                }
324                dist[i][j] = t1[i].distanceSquared(t2[j]);
325                if (dist[best_i][best_j] > dist[i][j]) {
326                    best_i = i;
327                    best_j = j;
328                }
329            }
330        }
331        if (best_i == 1 && best_j == 1) {
332            tStep /= 2;
333            if (tStep < FLT_EPSILON_HALF) {
334                break;
335            }
336            calcMask = (1 << 0) | (1 << 2) | (1 << 3) | (1 << 5);
337            continue;
338        }
339        if (best_i == 0) {
340            *t1Seed -= tStep;
341            t1[2] = t1[1];
342            t1[1] = t1[0];
343            calcMask = 1 << 0;
344        } else if (best_i == 2) {
345            *t1Seed += tStep;
346            t1[0] = t1[1];
347            t1[1] = t1[2];
348            calcMask = 1 << 2;
349        } else {
350            calcMask = 0;
351        }
352        if (best_j == 0) {
353            *t2Seed -= tStep;
354            t2[2] = t2[1];
355            t2[1] = t2[0];
356            calcMask |= 1 << 3;
357        } else if (best_j == 2) {
358            *t2Seed += tStep;
359            t2[0] = t2[1];
360            t2[1] = t2[2];
361            calcMask |= 1 << 5;
362        }
363    } while (true);
364#if ONE_OFF_DEBUG
365    SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCTION__,
366        t1Seed, t2Seed, t1[1].fX, t1[1].fY, t1[2].fX, t1[2].fY);
367#endif
368    return false;
369}
370
371static void lookNearEnd(const SkDQuad& q1, const SkDQuad& q2, int testT,
372        const SkIntersections& orig, bool swap, SkIntersections* i) {
373    if (orig.used() == 1 && orig[!swap][0] == testT) {
374        return;
375    }
376    if (orig.used() == 2 && orig[!swap][1] == testT) {
377        return;
378    }
379    SkDLine tmpLine;
380    int testTIndex = testT << 1;
381    tmpLine[0] = tmpLine[1] = q2[testTIndex];
382    tmpLine[1].fX += q2[1].fY - q2[testTIndex].fY;
383    tmpLine[1].fY -= q2[1].fX - q2[testTIndex].fX;
384    SkIntersections impTs;
385    impTs.intersectRay(q1, tmpLine);
386    for (int index = 0; index < impTs.used(); ++index) {
387        SkDPoint realPt = impTs.pt(index);
388        if (!tmpLine[0].approximatelyPEqual(realPt)) {
389            continue;
390        }
391        if (swap) {
392            i->insert(testT, impTs[0][index], tmpLine[0]);
393        } else {
394            i->insert(impTs[0][index], testT, tmpLine[0]);
395        }
396    }
397}
398
399int SkIntersections::intersect(const SkDQuad& q1, const SkDQuad& q2) {
400    fMax = 4;
401    // if the quads share an end point, check to see if they overlap
402    for (int i1 = 0; i1 < 3; i1 += 2) {
403        for (int i2 = 0; i2 < 3; i2 += 2) {
404            if (q1[i1].asSkPoint() == q2[i2].asSkPoint()) {
405                insert(i1 >> 1, i2 >> 1, q1[i1]);
406            }
407        }
408    }
409    SkASSERT(fUsed < 3);
410    if (only_end_pts_in_common(q1, q2)) {
411        return fUsed;
412    }
413    if (only_end_pts_in_common(q2, q1)) {
414        return fUsed;
415    }
416    // see if either quad is really a line
417    // FIXME: figure out why reduce step didn't find this earlier
418    if (is_linear(q1, q2, this)) {
419        return fUsed;
420    }
421    SkIntersections swapped;
422    swapped.setMax(fMax);
423    if (is_linear(q2, q1, &swapped)) {
424        swapped.swapPts();
425        *this = swapped;
426        return fUsed;
427    }
428    SkIntersections copyI(*this);
429    lookNearEnd(q1, q2, 0, *this, false, &copyI);
430    lookNearEnd(q1, q2, 1, *this, false, &copyI);
431    lookNearEnd(q2, q1, 0, *this, true, &copyI);
432    lookNearEnd(q2, q1, 1, *this, true, &copyI);
433    int innerEqual = 0;
434    if (copyI.fUsed >= 2) {
435        SkASSERT(copyI.fUsed <= 4);
436        double width = copyI[0][1] - copyI[0][0];
437        int midEnd = 1;
438        for (int index = 2; index < copyI.fUsed; ++index) {
439            double testWidth = copyI[0][index] - copyI[0][index - 1];
440            if (testWidth <= width) {
441                continue;
442            }
443            midEnd = index;
444        }
445        for (int index = 0; index < 2; ++index) {
446            double testT = (copyI[0][midEnd] * (index + 1)
447                    + copyI[0][midEnd - 1] * (2 - index)) / 3;
448            SkDPoint testPt1 = q1.ptAtT(testT);
449            testT = (copyI[1][midEnd] * (index + 1) + copyI[1][midEnd - 1] * (2 - index)) / 3;
450            SkDPoint testPt2 = q2.ptAtT(testT);
451            innerEqual += testPt1.approximatelyEqual(testPt2);
452        }
453    }
454    bool expectCoincident = copyI.fUsed >= 2 && innerEqual == 2;
455    if (expectCoincident) {
456        reset();
457        insertCoincident(copyI[0][0], copyI[1][0], copyI.fPt[0]);
458        int last = copyI.fUsed - 1;
459        insertCoincident(copyI[0][last], copyI[1][last], copyI.fPt[last]);
460        return fUsed;
461    }
462    SkDQuadImplicit i1(q1);
463    SkDQuadImplicit i2(q2);
464    int index;
465    bool flip1 = q1[2] == q2[0];
466    bool flip2 = q1[0] == q2[2];
467    bool useCubic = q1[0] == q2[0];
468    double roots1[4];
469    int rootCount = findRoots(i2, q1, roots1, useCubic, flip1, 0);
470    // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
471    double roots1Copy[4];
472    int r1Count = addValidRoots(roots1, rootCount, roots1Copy);
473    SkDPoint pts1[4];
474    for (index = 0; index < r1Count; ++index) {
475        pts1[index] = q1.ptAtT(roots1Copy[index]);
476    }
477    double roots2[4];
478    int rootCount2 = findRoots(i1, q2, roots2, useCubic, flip2, 0);
479    double roots2Copy[4];
480    int r2Count = addValidRoots(roots2, rootCount2, roots2Copy);
481    SkDPoint pts2[4];
482    for (index = 0; index < r2Count; ++index) {
483        pts2[index] = q2.ptAtT(roots2Copy[index]);
484    }
485    if (r1Count == r2Count && r1Count <= 1) {
486        if (r1Count == 1 && used() == 0) {
487            if (pts1[0].approximatelyEqual(pts2[0])) {
488                insert(roots1Copy[0], roots2Copy[0], pts1[0]);
489            } else {
490                // find intersection by chasing t
491                if (binary_search(q1, q2, roots1Copy, roots2Copy, pts1)) {
492                    insert(roots1Copy[0], roots2Copy[0], pts1[0]);
493                }
494            }
495        }
496        return fUsed;
497    }
498    int closest[4];
499    double dist[4];
500    bool foundSomething = false;
501    for (index = 0; index < r1Count; ++index) {
502        dist[index] = DBL_MAX;
503        closest[index] = -1;
504        for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) {
505            if (!pts2[ndex2].approximatelyEqual(pts1[index])) {
506                continue;
507            }
508            double dx = pts2[ndex2].fX - pts1[index].fX;
509            double dy = pts2[ndex2].fY - pts1[index].fY;
510            double distance = dx * dx + dy * dy;
511            if (dist[index] <= distance) {
512                continue;
513            }
514            for (int outer = 0; outer < index; ++outer) {
515                if (closest[outer] != ndex2) {
516                    continue;
517                }
518                if (dist[outer] < distance) {
519                    goto next;
520                }
521                closest[outer] = -1;
522            }
523            dist[index] = distance;
524            closest[index] = ndex2;
525            foundSomething = true;
526        next:
527            ;
528        }
529    }
530    if (r1Count && r2Count && !foundSomething) {
531        relaxed_is_linear(&q1, 0, 1, &q2, 0, 1, this);
532        return fUsed;
533    }
534    int used = 0;
535    do {
536        double lowest = DBL_MAX;
537        int lowestIndex = -1;
538        for (index = 0; index < r1Count; ++index) {
539            if (closest[index] < 0) {
540                continue;
541            }
542            if (roots1Copy[index] < lowest) {
543                lowestIndex = index;
544                lowest = roots1Copy[index];
545            }
546        }
547        if (lowestIndex < 0) {
548            break;
549        }
550        insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]],
551                pts1[lowestIndex]);
552        closest[lowestIndex] = -1;
553    } while (++used < r1Count);
554    return fUsed;
555}
556