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    // see parallel routine in line quadratic intersections
97    int intersectRay(double roots[3]) {
98        double adj = fLine[1].fX - fLine[0].fX;
99        double opp = fLine[1].fY - fLine[0].fY;
100        SkDCubic r;
101        for (int n = 0; n < 4; ++n) {
102            r[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp;
103        }
104        double A, B, C, D;
105        SkDCubic::Coefficients(&r[0].fX, &A, &B, &C, &D);
106        return SkDCubic::RootsValidT(A, B, C, D, roots);
107    }
108
109    int intersect() {
110        addExactEndPoints();
111        if (fAllowNear) {
112            addNearEndPoints();
113        }
114        double rootVals[3];
115        int roots = intersectRay(rootVals);
116        for (int index = 0; index < roots; ++index) {
117            double cubicT = rootVals[index];
118            double lineT = findLineT(cubicT);
119            SkDPoint pt;
120            if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized)) {
121    #if ONE_OFF_DEBUG
122                SkDPoint cPt = fCubic.ptAtT(cubicT);
123                SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
124                        cPt.fX, cPt.fY);
125    #endif
126                for (int inner = 0; inner < fIntersections->used(); ++inner) {
127                    if (fIntersections->pt(inner) != pt) {
128                        continue;
129                    }
130                    double existingCubicT = (*fIntersections)[0][inner];
131                    if (cubicT == existingCubicT) {
132                        goto skipInsert;
133                    }
134                    // check if midway on cubic is also same point. If so, discard this
135                    double cubicMidT = (existingCubicT + cubicT) / 2;
136                    SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
137                    if (cubicMidPt.approximatelyEqual(pt)) {
138                        goto skipInsert;
139                    }
140                }
141                fIntersections->insert(cubicT, lineT, pt);
142        skipInsert:
143                ;
144            }
145        }
146        return fIntersections->used();
147    }
148
149    int horizontalIntersect(double axisIntercept, double roots[3]) {
150        double A, B, C, D;
151        SkDCubic::Coefficients(&fCubic[0].fY, &A, &B, &C, &D);
152        D -= axisIntercept;
153        return SkDCubic::RootsValidT(A, B, C, D, roots);
154    }
155
156    int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
157        addExactHorizontalEndPoints(left, right, axisIntercept);
158        if (fAllowNear) {
159            addNearHorizontalEndPoints(left, right, axisIntercept);
160        }
161        double rootVals[3];
162        int roots = horizontalIntersect(axisIntercept, rootVals);
163        for (int index = 0; index < roots; ++index) {
164            double cubicT = rootVals[index];
165            SkDPoint pt = fCubic.ptAtT(cubicT);
166            double lineT = (pt.fX - left) / (right - left);
167            if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
168                fIntersections->insert(cubicT, lineT, pt);
169            }
170        }
171        if (flipped) {
172            fIntersections->flip();
173        }
174        return fIntersections->used();
175    }
176
177    int verticalIntersect(double axisIntercept, double roots[3]) {
178        double A, B, C, D;
179        SkDCubic::Coefficients(&fCubic[0].fX, &A, &B, &C, &D);
180        D -= axisIntercept;
181        return SkDCubic::RootsValidT(A, B, C, D, roots);
182    }
183
184    int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
185        addExactVerticalEndPoints(top, bottom, axisIntercept);
186        if (fAllowNear) {
187            addNearVerticalEndPoints(top, bottom, axisIntercept);
188        }
189        double rootVals[3];
190        int roots = verticalIntersect(axisIntercept, rootVals);
191        for (int index = 0; index < roots; ++index) {
192            double cubicT = rootVals[index];
193            SkDPoint pt = fCubic.ptAtT(cubicT);
194            double lineT = (pt.fY - top) / (bottom - top);
195            if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
196                fIntersections->insert(cubicT, lineT, pt);
197            }
198        }
199        if (flipped) {
200            fIntersections->flip();
201        }
202        return fIntersections->used();
203    }
204
205    protected:
206
207    void addExactEndPoints() {
208        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
209            double lineT = fLine.exactPoint(fCubic[cIndex]);
210            if (lineT < 0) {
211                continue;
212            }
213            double cubicT = (double) (cIndex >> 1);
214            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
215        }
216    }
217
218    /* Note that this does not look for endpoints of the line that are near the cubic.
219       These points are found later when check ends looks for missing points */
220    void addNearEndPoints() {
221        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
222            double cubicT = (double) (cIndex >> 1);
223            if (fIntersections->hasT(cubicT)) {
224                continue;
225            }
226            double lineT = fLine.nearPoint(fCubic[cIndex]);
227            if (lineT < 0) {
228                continue;
229            }
230            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
231        }
232    }
233
234    void addExactHorizontalEndPoints(double left, double right, double y) {
235        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
236            double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
237            if (lineT < 0) {
238                continue;
239            }
240            double cubicT = (double) (cIndex >> 1);
241            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
242        }
243    }
244
245    void addNearHorizontalEndPoints(double left, double right, double y) {
246        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
247            double cubicT = (double) (cIndex >> 1);
248            if (fIntersections->hasT(cubicT)) {
249                continue;
250            }
251            double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
252            if (lineT < 0) {
253                continue;
254            }
255            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
256        }
257        // FIXME: see if line end is nearly on cubic
258    }
259
260    void addExactVerticalEndPoints(double top, double bottom, double x) {
261        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
262            double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
263            if (lineT < 0) {
264                continue;
265            }
266            double cubicT = (double) (cIndex >> 1);
267            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
268        }
269    }
270
271    void addNearVerticalEndPoints(double top, double bottom, double x) {
272        for (int cIndex = 0; cIndex < 4; cIndex += 3) {
273            double cubicT = (double) (cIndex >> 1);
274            if (fIntersections->hasT(cubicT)) {
275                continue;
276            }
277            double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
278            if (lineT < 0) {
279                continue;
280            }
281            fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
282        }
283        // FIXME: see if line end is nearly on cubic
284    }
285
286    double findLineT(double t) {
287        SkDPoint xy = fCubic.ptAtT(t);
288        double dx = fLine[1].fX - fLine[0].fX;
289        double dy = fLine[1].fY - fLine[0].fY;
290        if (fabs(dx) > fabs(dy)) {
291            return (xy.fX - fLine[0].fX) / dx;
292        }
293        return (xy.fY - fLine[0].fY) / dy;
294    }
295
296    bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
297        if (!approximately_one_or_less(*lineT)) {
298            return false;
299        }
300        if (!approximately_zero_or_more(*lineT)) {
301            return false;
302        }
303        double cT = *cubicT = SkPinT(*cubicT);
304        double lT = *lineT = SkPinT(*lineT);
305        if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
306            *pt = fLine.ptAtT(lT);
307        } else if (ptSet == kPointUninitialized) {
308            *pt = fCubic.ptAtT(cT);
309        }
310        SkPoint gridPt = pt->asSkPoint();
311        if (gridPt == fLine[0].asSkPoint()) {
312            *lineT = 0;
313        } else if (gridPt == fLine[1].asSkPoint()) {
314            *lineT = 1;
315        }
316        if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
317            *cubicT = 0;
318        } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
319            *cubicT = 1;
320        }
321        return true;
322    }
323
324private:
325    const SkDCubic& fCubic;
326    const SkDLine& fLine;
327    SkIntersections* fIntersections;
328    bool fAllowNear;
329};
330
331int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
332        bool flipped) {
333    SkDLine line = {{{ left, y }, { right, y }}};
334    LineCubicIntersections c(cubic, line, this);
335    return c.horizontalIntersect(y, left, right, flipped);
336}
337
338int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
339        bool flipped) {
340    SkDLine line = {{{ x, top }, { x, bottom }}};
341    LineCubicIntersections c(cubic, line, this);
342    return c.verticalIntersect(x, top, bottom, flipped);
343}
344
345int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
346    LineCubicIntersections c(cubic, line, this);
347    c.allowNear(fAllowNear);
348    return c.intersect();
349}
350
351int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
352    LineCubicIntersections c(cubic, line, this);
353    fUsed = c.intersectRay(fT[0]);
354    for (int index = 0; index < fUsed; ++index) {
355        fPt[index] = cubic.ptAtT(fT[0][index]);
356    }
357    return fUsed;
358}
359