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