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 "SkPathOpsLine.h"
9
10/* Determine the intersection point of two lines. This assumes the lines are not parallel,
11   and that that the lines are infinite.
12   From http://en.wikipedia.org/wiki/Line-line_intersection
13 */
14SkDPoint SkIntersections::Line(const SkDLine& a, const SkDLine& b) {
15    double axLen = a[1].fX - a[0].fX;
16    double ayLen = a[1].fY - a[0].fY;
17    double bxLen = b[1].fX - b[0].fX;
18    double byLen = b[1].fY - b[0].fY;
19    double denom = byLen * axLen - ayLen * bxLen;
20    SkASSERT(denom);
21    double term1 = a[1].fX * a[0].fY - a[1].fY * a[0].fX;
22    double term2 = b[1].fX * b[0].fY - b[1].fY * b[0].fX;
23    SkDPoint p;
24    p.fX = (term1 * bxLen - axLen * term2) / denom;
25    p.fY = (term1 * byLen - ayLen * term2) / denom;
26    return p;
27}
28
29int SkIntersections::computePoints(const SkDLine& line, int used) {
30    fPt[0] = line.ptAtT(fT[0][0]);
31    if ((fUsed = used) == 2) {
32        fPt[1] = line.ptAtT(fT[0][1]);
33    }
34    return fUsed;
35}
36
37int SkIntersections::intersectRay(const SkDLine& a, const SkDLine& b) {
38    double axLen = a[1].fX - a[0].fX;
39    double ayLen = a[1].fY - a[0].fY;
40    double bxLen = b[1].fX - b[0].fX;
41    double byLen = b[1].fY - b[0].fY;
42    /* Slopes match when denom goes to zero:
43                      axLen / ayLen ==                   bxLen / byLen
44    (ayLen * byLen) * axLen / ayLen == (ayLen * byLen) * bxLen / byLen
45             byLen  * axLen         ==  ayLen          * bxLen
46             byLen  * axLen         -   ayLen          * bxLen == 0 ( == denom )
47     */
48    double denom = byLen * axLen - ayLen * bxLen;
49    double ab0y = a[0].fY - b[0].fY;
50    double ab0x = a[0].fX - b[0].fX;
51    double numerA = ab0y * bxLen - byLen * ab0x;
52    double numerB = ab0y * axLen - ayLen * ab0x;
53    numerA /= denom;
54    numerB /= denom;
55    int used;
56    if (!approximately_zero(denom)) {
57        fT[0][0] = numerA;
58        fT[1][0] = numerB;
59        used = 1;
60    } else {
61       /* See if the axis intercepts match:
62                  ay - ax * ayLen / axLen  ==          by - bx * ayLen / axLen
63         axLen * (ay - ax * ayLen / axLen) == axLen * (by - bx * ayLen / axLen)
64         axLen *  ay - ax * ayLen          == axLen *  by - bx * ayLen
65        */
66        if (!AlmostEqualUlps(axLen * a[0].fY - ayLen * a[0].fX,
67                axLen * b[0].fY - ayLen * b[0].fX)) {
68            return fUsed = 0;
69        }
70        // there's no great answer for intersection points for coincident rays, but return something
71        fT[0][0] = fT[1][0] = 0;
72        fT[1][0] = fT[1][1] = 1;
73        used = 2;
74    }
75    return computePoints(a, used);
76}
77
78// note that this only works if both lines are neither horizontal nor vertical
79int SkIntersections::intersect(const SkDLine& a, const SkDLine& b) {
80    // see if end points intersect the opposite line
81    double t;
82    for (int iA = 0; iA < 2; ++iA) {
83        if ((t = b.exactPoint(a[iA])) >= 0) {
84            insert(iA, t, a[iA]);
85        }
86    }
87    for (int iB = 0; iB < 2; ++iB) {
88        if ((t = a.exactPoint(b[iB])) >= 0) {
89            insert(t, iB, b[iB]);
90        }
91    }
92    /* Determine the intersection point of two line segments
93       Return FALSE if the lines don't intersect
94       from: http://paulbourke.net/geometry/lineline2d/ */
95    double axLen = a[1].fX - a[0].fX;
96    double ayLen = a[1].fY - a[0].fY;
97    double bxLen = b[1].fX - b[0].fX;
98    double byLen = b[1].fY - b[0].fY;
99    /* Slopes match when denom goes to zero:
100                      axLen / ayLen ==                   bxLen / byLen
101    (ayLen * byLen) * axLen / ayLen == (ayLen * byLen) * bxLen / byLen
102             byLen  * axLen         ==  ayLen          * bxLen
103             byLen  * axLen         -   ayLen          * bxLen == 0 ( == denom )
104     */
105    double axByLen = axLen * byLen;
106    double ayBxLen = ayLen * bxLen;
107    // detect parallel lines the same way here and in SkOpAngle operator <
108    // so that non-parallel means they are also sortable
109    bool parallel = AlmostEqualUlps(axByLen, ayBxLen);
110    if (!parallel) {
111        double ab0y = a[0].fY - b[0].fY;
112        double ab0x = a[0].fX - b[0].fX;
113        double numerA = ab0y * bxLen - byLen * ab0x;
114        double numerB = ab0y * axLen - ayLen * ab0x;
115        double denom = axByLen - ayBxLen;
116        if (between(0, numerA, denom) && between(0, numerB, denom)) {
117            fT[0][0] = numerA / denom;
118            fT[1][0] = numerB / denom;
119            computePoints(a, 1);
120        }
121    }
122    if (fAllowNear || parallel) {
123        for (int iA = 0; iA < 2; ++iA) {
124            if ((t = b.nearPoint(a[iA])) >= 0) {
125                insert(iA, t, a[iA]);
126            }
127        }
128        for (int iB = 0; iB < 2; ++iB) {
129            if ((t = a.nearPoint(b[iB])) >= 0) {
130                insert(t, iB, b[iB]);
131            }
132        }
133    }
134    return fUsed;
135}
136
137static int horizontal_coincident(const SkDLine& line, double y) {
138    double min = line[0].fY;
139    double max = line[1].fY;
140    if (min > max) {
141        SkTSwap(min, max);
142    }
143    if (min > y || max < y) {
144        return 0;
145    }
146    if (AlmostEqualUlps(min, max) && max - min < fabs(line[0].fX - line[1].fX)) {
147        return 2;
148    }
149    return 1;
150}
151
152static double horizontal_intercept(const SkDLine& line, double y) {
153     return (y - line[0].fY) / (line[1].fY - line[0].fY);
154}
155
156int SkIntersections::horizontal(const SkDLine& line, double y) {
157    int horizontalType = horizontal_coincident(line, y);
158    if (horizontalType == 1) {
159        fT[0][0] = horizontal_intercept(line, y);
160    } else if (horizontalType == 2) {
161        fT[0][0] = 0;
162        fT[0][1] = 1;
163    }
164    return fUsed = horizontalType;
165}
166
167int SkIntersections::horizontal(const SkDLine& line, double left, double right,
168                                double y, bool flipped) {
169    // see if end points intersect the opposite line
170    double t;
171    const SkDPoint leftPt = { left, y };
172    if ((t = line.exactPoint(leftPt)) >= 0) {
173        insert(t, (double) flipped, leftPt);
174    }
175    if (left != right) {
176        const SkDPoint rightPt = { right, y };
177        if ((t = line.exactPoint(rightPt)) >= 0) {
178            insert(t, (double) !flipped, rightPt);
179        }
180        for (int index = 0; index < 2; ++index) {
181            if ((t = SkDLine::ExactPointH(line[index], left, right, y)) >= 0) {
182                insert((double) index, flipped ? 1 - t : t, line[index]);
183            }
184        }
185    }
186    int result = horizontal_coincident(line, y);
187    if (result == 1 && fUsed == 0) {
188        fT[0][0] = horizontal_intercept(line, y);
189        double xIntercept = line[0].fX + fT[0][0] * (line[1].fX - line[0].fX);
190        if (between(left, xIntercept, right)) {
191            fT[1][0] = (xIntercept - left) / (right - left);
192            if (flipped) {
193                // OPTIMIZATION: ? instead of swapping, pass original line, use [1].fX - [0].fX
194                for (int index = 0; index < result; ++index) {
195                    fT[1][index] = 1 - fT[1][index];
196                }
197            }
198            return computePoints(line, result);
199        }
200    }
201    if (!fAllowNear && result != 2) {
202        return fUsed;
203    }
204    if ((t = line.nearPoint(leftPt)) >= 0) {
205        insert(t, (double) flipped, leftPt);
206    }
207    if (left != right) {
208        const SkDPoint rightPt = { right, y };
209        if ((t = line.nearPoint(rightPt)) >= 0) {
210            insert(t, (double) !flipped, rightPt);
211        }
212        for (int index = 0; index < 2; ++index) {
213            if ((t = SkDLine::NearPointH(line[index], left, right, y)) >= 0) {
214                insert((double) index, flipped ? 1 - t : t, line[index]);
215            }
216        }
217    }
218    return fUsed;
219}
220
221static int vertical_coincident(const SkDLine& line, double x) {
222    double min = line[0].fX;
223    double max = line[1].fX;
224    if (min > max) {
225        SkTSwap(min, max);
226    }
227    if (!precisely_between(min, x, max)) {
228        return 0;
229    }
230    if (AlmostEqualUlps(min, max)) {
231        return 2;
232    }
233    return 1;
234}
235
236static double vertical_intercept(const SkDLine& line, double x) {
237    return (x - line[0].fX) / (line[1].fX - line[0].fX);
238}
239
240int SkIntersections::vertical(const SkDLine& line, double x) {
241    int verticalType = vertical_coincident(line, x);
242    if (verticalType == 1) {
243        fT[0][0] = vertical_intercept(line, x);
244    } else if (verticalType == 2) {
245        fT[0][0] = 0;
246        fT[0][1] = 1;
247    }
248    return fUsed = verticalType;
249}
250
251int SkIntersections::vertical(const SkDLine& line, double top, double bottom,
252                              double x, bool flipped) {
253    // see if end points intersect the opposite line
254    double t;
255    SkDPoint topPt = { x, top };
256    if ((t = line.exactPoint(topPt)) >= 0) {
257        insert(t, (double) flipped, topPt);
258    }
259    if (top != bottom) {
260        SkDPoint bottomPt = { x, bottom };
261        if ((t = line.exactPoint(bottomPt)) >= 0) {
262            insert(t, (double) !flipped, bottomPt);
263        }
264        for (int index = 0; index < 2; ++index) {
265            if ((t = SkDLine::ExactPointV(line[index], top, bottom, x)) >= 0) {
266                insert((double) index, flipped ? 1 - t : t, line[index]);
267            }
268        }
269    }
270    int result = vertical_coincident(line, x);
271    if (result == 1 && fUsed == 0) {
272        fT[0][0] = vertical_intercept(line, x);
273        double yIntercept = line[0].fY + fT[0][0] * (line[1].fY - line[0].fY);
274        if (between(top, yIntercept, bottom)) {
275            fT[1][0] = (yIntercept - top) / (bottom - top);
276            if (flipped) {
277                // OPTIMIZATION: instead of swapping, pass original line, use [1].fY - [0].fY
278                for (int index = 0; index < result; ++index) {
279                    fT[1][index] = 1 - fT[1][index];
280                }
281            }
282            return computePoints(line, result);
283        }
284    }
285    if (!fAllowNear && result != 2) {
286        return fUsed;
287    }
288    if ((t = line.nearPoint(topPt)) >= 0) {
289        insert(t, (double) flipped, topPt);
290    }
291    if (top != bottom) {
292        SkDPoint bottomPt = { x, bottom };
293        if ((t = line.nearPoint(bottomPt)) >= 0) {
294            insert(t, (double) !flipped, bottomPt);
295        }
296        for (int index = 0; index < 2; ++index) {
297            if ((t = SkDLine::NearPointV(line[index], top, bottom, x)) >= 0) {
298                insert((double) index, flipped ? 1 - t : t, line[index]);
299            }
300        }
301    }
302    return fUsed;
303}
304
305// from http://www.bryceboe.com/wordpress/wp-content/uploads/2006/10/intersect.py
306// 4 subs, 2 muls, 1 cmp
307static bool ccw(const SkDPoint& A, const SkDPoint& B, const SkDPoint& C) {
308    return (C.fY - A.fY) * (B.fX - A.fX) > (B.fY - A.fY) * (C.fX - A.fX);
309}
310
311// 16 subs, 8 muls, 6 cmps
312bool SkIntersections::Test(const SkDLine& a, const SkDLine& b) {
313    return ccw(a[0], b[0], b[1]) != ccw(a[1], b[0], b[1])
314            && ccw(a[0], a[1], b[0]) != ccw(a[0], a[1], b[1]);
315}
316