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 "SkIntersections.h"
8#include "SkPathOpsCubic.h"
9#include "SkPathOpsCurve.h"
10#include "SkPathOpsLine.h"
11
12/*
13Find the interection of a line and cubic by solving for valid t values.
14
15Analogous to line-quadratic intersection, solve line-cubic intersection by
16representing the cubic as:
17  x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
18  y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
19and the line as:
20  y = i*x + j  (if the line is more horizontal)
21or:
22  x = i*y + j  (if the line is more vertical)
23
24Then using Mathematica, solve for the values of t where the cubic intersects the
25line:
26
27  (in) Resultant[
28        a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
29        e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
30  (out) -e     +   j     +
31       3 e t   - 3 f t   -
32       3 e t^2 + 6 f t^2 - 3 g t^2 +
33         e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
34     i ( a     -
35       3 a t + 3 b t +
36       3 a t^2 - 6 b t^2 + 3 c t^2 -
37         a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
38
39if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
40
41  (in) Resultant[
42        a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
43        e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y,       y]
44  (out)  a     -   j     -
45       3 a t   + 3 b t   +
46       3 a t^2 - 6 b t^2 + 3 c t^2 -
47         a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
48     i ( e     -
49       3 e t   + 3 f t   +
50       3 e t^2 - 6 f t^2 + 3 g t^2 -
51         e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
52
53Solving this with Mathematica produces an expression with hundreds of terms;
54instead, use Numeric Solutions recipe to solve the cubic.
55
56The near-horizontal case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
57    A =   (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d)     )
58    B = 3*(-( e - 2*f +   g    ) + i*( a - 2*b +   c    )     )
59    C = 3*(-(-e +   f          ) + i*(-a +   b          )     )
60    D =   (-( e                ) + i*( a                ) + j )
61
62The near-vertical case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
63    A =   ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h)     )
64    B = 3*( ( a - 2*b +   c    ) - i*( e - 2*f +   g    )     )
65    C = 3*( (-a +   b          ) - i*(-e +   f          )     )
66    D =   ( ( a                ) - i*( e                ) - j )
67
68For horizontal lines:
69(in) Resultant[
70      a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
71      e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
72(out)  e     -   j     -
73     3 e t   + 3 f t   +
74     3 e t^2 - 6 f t^2 + 3 g t^2 -
75       e t^3 + 3 f t^3 - 3 g t^3 + h t^3
76 */
77
78class LineCubicIntersections {
79public:
80    enum PinTPoint {
81        kPointUninitialized,
82        kPointInitialized
83    };
84
85    LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections* i)
86        : fCubic(c)
87        , fLine(l)
88        , fIntersections(i)
89        , fAllowNear(true) {
90        i->setMax(4);
91    }
92
93    void allowNear(bool allow) {
94        fAllowNear = allow;
95    }
96
97    void checkCoincident() {
98        int last = fIntersections->used() - 1;
99        for (int index = 0; index < last; ) {
100            double cubicMidT = ((*fIntersections)[0][index] + (*fIntersections)[0][index + 1]) / 2;
101            SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
102            double t = fLine.nearPoint(cubicMidPt, nullptr);
103            if (t < 0) {
104                ++index;
105                continue;
106            }
107            if (fIntersections->isCoincident(index)) {
108                fIntersections->removeOne(index);
109                --last;
110            } else if (fIntersections->isCoincident(index + 1)) {
111                fIntersections->removeOne(index + 1);
112                --last;
113            } else {
114                fIntersections->setCoincident(index++);
115            }
116            fIntersections->setCoincident(index);
117        }
118    }
119
120    // see parallel routine in line quadratic intersections
121    int intersectRay(double roots[3]) {
122        double adj = fLine[1].fX - fLine[0].fX;
123        double opp = fLine[1].fY - fLine[0].fY;
124        SkDCubic c;
125        SkDEBUGCODE(c.fDebugGlobalState = fIntersections->globalState());
126        for (int n = 0; n < 4; ++n) {
127            c[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp;
128        }
129        double A, B, C, D;
130        SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D);
131        int count = SkDCubic::RootsValidT(A, B, C, D, roots);
132        for (int index = 0; index < count; ++index) {
133            SkDPoint calcPt = c.ptAtT(roots[index]);
134            if (!approximately_zero(calcPt.fX)) {
135                for (int n = 0; n < 4; ++n) {
136                    c[n].fY = (fCubic[n].fY - fLine[0].fY) * opp
137                            + (fCubic[n].fX - fLine[0].fX) * adj;
138                }
139                double extremeTs[6];
140                int extrema = SkDCubic::FindExtrema(&c[0].fX, extremeTs);
141                count = c.searchRoots(extremeTs, extrema, 0, SkDCubic::kXAxis, roots);
142                break;
143            }
144        }
145        return count;
146    }
147
148    int intersect() {
149        addExactEndPoints();
150        if (fAllowNear) {
151            addNearEndPoints();
152        }
153        double rootVals[3];
154        int roots = intersectRay(rootVals);
155        for (int index = 0; index < roots; ++index) {
156            double cubicT = rootVals[index];
157            double lineT = findLineT(cubicT);
158            SkDPoint pt;
159            if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized) && uniqueAnswer(cubicT, pt)) {
160                fIntersections->insert(cubicT, lineT, pt);
161            }
162        }
163        checkCoincident();
164        return fIntersections->used();
165    }
166
167    static int HorizontalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) {
168        double A, B, C, D;
169        SkDCubic::Coefficients(&c[0].fY, &A, &B, &C, &D);
170        D -= axisIntercept;
171        int count = SkDCubic::RootsValidT(A, B, C, D, roots);
172        for (int index = 0; index < count; ++index) {
173            SkDPoint calcPt = c.ptAtT(roots[index]);
174            if (!approximately_equal(calcPt.fY, axisIntercept)) {
175                double extremeTs[6];
176                int extrema = SkDCubic::FindExtrema(&c[0].fY, extremeTs);
177                count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kYAxis, roots);
178                break;
179            }
180        }
181        return count;
182    }
183
184    int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
185        addExactHorizontalEndPoints(left, right, axisIntercept);
186        if (fAllowNear) {
187            addNearHorizontalEndPoints(left, right, axisIntercept);
188        }
189        double roots[3];
190        int count = HorizontalIntersect(fCubic, axisIntercept, roots);
191        for (int index = 0; index < count; ++index) {
192            double cubicT = roots[index];
193            SkDPoint pt = { fCubic.ptAtT(cubicT).fX,  axisIntercept };
194            double lineT = (pt.fX - left) / (right - left);
195            if (pinTs(&cubicT, &lineT, &pt, kPointInitialized) && uniqueAnswer(cubicT, pt)) {
196                fIntersections->insert(cubicT, lineT, pt);
197            }
198        }
199        if (flipped) {
200            fIntersections->flip();
201        }
202        checkCoincident();
203        return fIntersections->used();
204    }
205
206        bool uniqueAnswer(double cubicT, const SkDPoint& pt) {
207            for (int inner = 0; inner < fIntersections->used(); ++inner) {
208                if (fIntersections->pt(inner) != pt) {
209                    continue;
210                }
211                double existingCubicT = (*fIntersections)[0][inner];
212                if (cubicT == existingCubicT) {
213                    return false;
214                }
215                // check if midway on cubic is also same point. If so, discard this
216                double cubicMidT = (existingCubicT + cubicT) / 2;
217                SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
218                if (cubicMidPt.approximatelyEqual(pt)) {
219                    return false;
220                }
221            }
222#if ONE_OFF_DEBUG
223            SkDPoint cPt = fCubic.ptAtT(cubicT);
224            SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
225                    cPt.fX, cPt.fY);
226#endif
227            return true;
228        }
229
230    static int VerticalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) {
231        double A, B, C, D;
232        SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D);
233        D -= axisIntercept;
234        int count = SkDCubic::RootsValidT(A, B, C, D, roots);
235        for (int index = 0; index < count; ++index) {
236            SkDPoint calcPt = c.ptAtT(roots[index]);
237            if (!approximately_equal(calcPt.fX, axisIntercept)) {
238                double extremeTs[6];
239                int extrema = SkDCubic::FindExtrema(&c[0].fX, extremeTs);
240                count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kXAxis, roots);
241                break;
242            }
243        }
244        return count;
245    }
246
247    int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
248        addExactVerticalEndPoints(top, bottom, axisIntercept);
249        if (fAllowNear) {
250            addNearVerticalEndPoints(top, bottom, axisIntercept);
251        }
252        double roots[3];
253        int count = VerticalIntersect(fCubic, axisIntercept, roots);
254        for (int index = 0; index < count; ++index) {
255            double cubicT = roots[index];
256            SkDPoint pt = { axisIntercept, fCubic.ptAtT(cubicT).fY };
257            double lineT = (pt.fY - top) / (bottom - top);
258            if (pinTs(&cubicT, &lineT, &pt, kPointInitialized) && uniqueAnswer(cubicT, pt)) {
259                fIntersections->insert(cubicT, lineT, pt);
260            }
261        }
262        if (flipped) {
263            fIntersections->flip();
264        }
265        checkCoincident();
266        return fIntersections->used();
267    }
268
269    protected:
270
271    void addExactEndPoints() {
272        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
273            double lineT = fLine.exactPoint(fCubic[cIndex]);
274            if (lineT < 0) {
275                continue;
276            }
277            double cubicT = (double) (cIndex >> 1);
278            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
279        }
280    }
281
282    /* Note that this does not look for endpoints of the line that are near the cubic.
283       These points are found later when check ends looks for missing points */
284    void addNearEndPoints() {
285        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
286            double cubicT = (double) (cIndex >> 1);
287            if (fIntersections->hasT(cubicT)) {
288                continue;
289            }
290            double lineT = fLine.nearPoint(fCubic[cIndex], nullptr);
291            if (lineT < 0) {
292                continue;
293            }
294            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
295        }
296        this->addLineNearEndPoints();
297    }
298
299    void addLineNearEndPoints() {
300        for (int lIndex = 0; lIndex < 2; ++lIndex) {
301            double lineT = (double) lIndex;
302            if (fIntersections->hasOppT(lineT)) {
303                continue;
304            }
305            double cubicT = ((SkDCurve*) &fCubic)->nearPoint(SkPath::kCubic_Verb,
306                fLine[lIndex], fLine[!lIndex]);
307            if (cubicT < 0) {
308                continue;
309            }
310            fIntersections->insert(cubicT, lineT, fLine[lIndex]);
311        }
312    }
313
314    void addExactHorizontalEndPoints(double left, double right, double y) {
315        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
316            double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
317            if (lineT < 0) {
318                continue;
319            }
320            double cubicT = (double) (cIndex >> 1);
321            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
322        }
323    }
324
325    void addNearHorizontalEndPoints(double left, double right, double y) {
326        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
327            double cubicT = (double) (cIndex >> 1);
328            if (fIntersections->hasT(cubicT)) {
329                continue;
330            }
331            double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
332            if (lineT < 0) {
333                continue;
334            }
335            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
336        }
337        this->addLineNearEndPoints();
338    }
339
340    void addExactVerticalEndPoints(double top, double bottom, double x) {
341        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
342            double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
343            if (lineT < 0) {
344                continue;
345            }
346            double cubicT = (double) (cIndex >> 1);
347            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
348        }
349    }
350
351    void addNearVerticalEndPoints(double top, double bottom, double x) {
352        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
353            double cubicT = (double) (cIndex >> 1);
354            if (fIntersections->hasT(cubicT)) {
355                continue;
356            }
357            double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
358            if (lineT < 0) {
359                continue;
360            }
361            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
362        }
363        this->addLineNearEndPoints();
364    }
365
366    double findLineT(double t) {
367        SkDPoint xy = fCubic.ptAtT(t);
368        double dx = fLine[1].fX - fLine[0].fX;
369        double dy = fLine[1].fY - fLine[0].fY;
370        if (fabs(dx) > fabs(dy)) {
371            return (xy.fX - fLine[0].fX) / dx;
372        }
373        return (xy.fY - fLine[0].fY) / dy;
374    }
375
376    bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
377        if (!approximately_one_or_less(*lineT)) {
378            return false;
379        }
380        if (!approximately_zero_or_more(*lineT)) {
381            return false;
382        }
383        double cT = *cubicT = SkPinT(*cubicT);
384        double lT = *lineT = SkPinT(*lineT);
385        SkDPoint lPt = fLine.ptAtT(lT);
386        SkDPoint cPt = fCubic.ptAtT(cT);
387        if (!lPt.roughlyEqual(cPt)) {
388            return false;
389        }
390        // FIXME: if points are roughly equal but not approximately equal, need to do
391        // a binary search like quad/quad intersection to find more precise t values
392        if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
393            *pt = lPt;
394        } else if (ptSet == kPointUninitialized) {
395            *pt = cPt;
396        }
397        SkPoint gridPt = pt->asSkPoint();
398        if (gridPt == fLine[0].asSkPoint()) {
399            *lineT = 0;
400        } else if (gridPt == fLine[1].asSkPoint()) {
401            *lineT = 1;
402        }
403        if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
404            *cubicT = 0;
405        } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
406            *cubicT = 1;
407        }
408        return true;
409    }
410
411private:
412    const SkDCubic& fCubic;
413    const SkDLine& fLine;
414    SkIntersections* fIntersections;
415    bool fAllowNear;
416};
417
418int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
419        bool flipped) {
420    SkDLine line = {{{ left, y }, { right, y }}};
421    LineCubicIntersections c(cubic, line, this);
422    return c.horizontalIntersect(y, left, right, flipped);
423}
424
425int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
426        bool flipped) {
427    SkDLine line = {{{ x, top }, { x, bottom }}};
428    LineCubicIntersections c(cubic, line, this);
429    return c.verticalIntersect(x, top, bottom, flipped);
430}
431
432int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
433    LineCubicIntersections c(cubic, line, this);
434    c.allowNear(fAllowNear);
435    return c.intersect();
436}
437
438int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
439    LineCubicIntersections c(cubic, line, this);
440    fUsed = c.intersectRay(fT[0]);
441    for (int index = 0; index < fUsed; ++index) {
442        fPt[index] = cubic.ptAtT(fT[0][index]);
443    }
444    return fUsed;
445}
446
447// SkDCubic accessors to Intersection utilities
448
449int SkDCubic::horizontalIntersect(double yIntercept, double roots[3]) const {
450    return LineCubicIntersections::HorizontalIntersect(*this, yIntercept, roots);
451}
452
453int SkDCubic::verticalIntersect(double xIntercept, double roots[3]) const {
454    return LineCubicIntersections::VerticalIntersect(*this, xIntercept, roots);
455}
456