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 "CurveIntersection.h"
8#include "Intersections.h"
9#include "IntersectionUtilities.h"
10#include "LineIntersection.h"
11#include "LineUtilities.h"
12#include "QuadraticLineSegments.h"
13#include "QuadraticUtilities.h"
14#include <algorithm> // for swap
15
16static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
17
18class QuadraticIntersections {
19public:
20
21QuadraticIntersections(const Quadratic& q1, const Quadratic& q2, Intersections& i)
22    : quad1(q1)
23    , quad2(q2)
24    , intersections(i)
25    , depth(0)
26    , splits(0)
27    , coinMinT1(-1) {
28}
29
30bool intersect() {
31    double minT1, minT2, maxT1, maxT2;
32    if (!bezier_clip(quad2, quad1, minT1, maxT1)) {
33        return false;
34    }
35    if (!bezier_clip(quad1, quad2, minT2, maxT2)) {
36        return false;
37    }
38    quad1Divisions = 1 / subDivisions(quad1);
39    quad2Divisions = 1 / subDivisions(quad2);
40    int split;
41    if (maxT1 - minT1 < maxT2 - minT2) {
42        intersections.swap();
43        minT2 = 0;
44        maxT2 = 1;
45        split = maxT1 - minT1 > tClipLimit;
46    } else {
47        minT1 = 0;
48        maxT1 = 1;
49        split = (maxT2 - minT2 > tClipLimit) << 1;
50    }
51    return chop(minT1, maxT1, minT2, maxT2, split);
52}
53
54protected:
55
56bool intersect(double minT1, double maxT1, double minT2, double maxT2) {
57    bool t1IsLine = maxT1 - minT1 <= quad1Divisions;
58    bool t2IsLine = maxT2 - minT2 <= quad2Divisions;
59    if (t1IsLine | t2IsLine) {
60        return intersectAsLine(minT1, maxT1, minT2, maxT2, t1IsLine, t2IsLine);
61    }
62    Quadratic smaller, larger;
63    // FIXME: carry last subdivide and reduceOrder result with quad
64    sub_divide(quad1, minT1, maxT1, intersections.swapped() ? larger : smaller);
65    sub_divide(quad2, minT2, maxT2, intersections.swapped() ? smaller : larger);
66    double minT, maxT;
67    if (!bezier_clip(smaller, larger, minT, maxT)) {
68        if (approximately_equal(minT, maxT)) {
69            double smallT, largeT;
70            _Point q2pt, q1pt;
71            if (intersections.swapped()) {
72                largeT = interp(minT2, maxT2, minT);
73                xy_at_t(quad2, largeT, q2pt.x, q2pt.y);
74                xy_at_t(quad1, minT1, q1pt.x, q1pt.y);
75                if (AlmostEqualUlps(q2pt.x, q1pt.x) && AlmostEqualUlps(q2pt.y, q1pt.y)) {
76                    smallT = minT1;
77                } else {
78                    xy_at_t(quad1, maxT1, q1pt.x, q1pt.y); // FIXME: debug code
79                    SkASSERT(AlmostEqualUlps(q2pt.x, q1pt.x) && AlmostEqualUlps(q2pt.y, q1pt.y));
80                    smallT = maxT1;
81                }
82            } else {
83                smallT = interp(minT1, maxT1, minT);
84                xy_at_t(quad1, smallT, q1pt.x, q1pt.y);
85                xy_at_t(quad2, minT2, q2pt.x, q2pt.y);
86                if (AlmostEqualUlps(q2pt.x, q1pt.x) && AlmostEqualUlps(q2pt.y, q1pt.y)) {
87                    largeT = minT2;
88                } else {
89                    xy_at_t(quad2, maxT2, q2pt.x, q2pt.y); // FIXME: debug code
90                    SkASSERT(AlmostEqualUlps(q2pt.x, q1pt.x) && AlmostEqualUlps(q2pt.y, q1pt.y));
91                    largeT = maxT2;
92                }
93            }
94            intersections.add(smallT, largeT);
95            return true;
96        }
97        return false;
98    }
99    int split;
100    if (intersections.swapped()) {
101        double newMinT1 = interp(minT1, maxT1, minT);
102        double newMaxT1 = interp(minT1, maxT1, maxT);
103        split = (newMaxT1 - newMinT1 > (maxT1 - minT1) * tClipLimit) << 1;
104#define VERBOSE 0
105#if VERBOSE
106        printf("%s d=%d s=%d new1=(%g,%g) old1=(%g,%g) split=%d\n", __FUNCTION__, depth,
107            splits, newMinT1, newMaxT1, minT1, maxT1, split);
108#endif
109        minT1 = newMinT1;
110        maxT1 = newMaxT1;
111    } else {
112        double newMinT2 = interp(minT2, maxT2, minT);
113        double newMaxT2 = interp(minT2, maxT2, maxT);
114        split = newMaxT2 - newMinT2 > (maxT2 - minT2) * tClipLimit;
115#if VERBOSE
116        printf("%s d=%d s=%d new2=(%g,%g) old2=(%g,%g) split=%d\n", __FUNCTION__, depth,
117            splits, newMinT2, newMaxT2, minT2, maxT2, split);
118#endif
119        minT2 = newMinT2;
120        maxT2 = newMaxT2;
121    }
122    return chop(minT1, maxT1, minT2, maxT2, split);
123}
124
125bool intersectAsLine(double minT1, double maxT1, double minT2, double maxT2,
126       bool treat1AsLine, bool treat2AsLine)
127{
128    _Line line1, line2;
129    if (intersections.swapped()) {
130        SkTSwap(treat1AsLine, treat2AsLine);
131        SkTSwap(minT1, minT2);
132        SkTSwap(maxT1, maxT2);
133    }
134    if (coinMinT1 >= 0) {
135        bool earlyExit;
136        if ((earlyExit = coinMaxT1 == minT1)) {
137            coinMaxT1 = maxT1;
138        }
139        if (coinMaxT2 == minT2) {
140            coinMaxT2 = maxT2;
141            return true;
142        }
143        if (earlyExit) {
144            return true;
145        }
146        coinMinT1 = -1;
147    }
148    // do line/quadratic or even line/line intersection instead
149    if (treat1AsLine) {
150        xy_at_t(quad1, minT1, line1[0].x, line1[0].y);
151        xy_at_t(quad1, maxT1, line1[1].x, line1[1].y);
152    }
153    if (treat2AsLine) {
154        xy_at_t(quad2, minT2, line2[0].x, line2[0].y);
155        xy_at_t(quad2, maxT2, line2[1].x, line2[1].y);
156    }
157    int pts;
158    double smallT1, largeT1, smallT2, largeT2;
159    if (treat1AsLine & treat2AsLine) {
160        double t1[2], t2[2];
161        pts = ::intersect(line1, line2, t1, t2);
162        if (pts == 2) {
163            smallT1 = interp(minT1, maxT1, t1[0]);
164            largeT1 = interp(minT2, maxT2, t2[0]);
165            smallT2 = interp(minT1, maxT1, t1[1]);
166            largeT2 = interp(minT2, maxT2, t2[1]);
167            intersections.addCoincident(smallT1, smallT2, largeT1, largeT2);
168        } else {
169            smallT1 = interp(minT1, maxT1, t1[0]);
170            largeT1 = interp(minT2, maxT2, t2[0]);
171            intersections.add(smallT1, largeT1);
172        }
173    } else {
174        Intersections lq;
175        pts = ::intersect(treat1AsLine ? quad2 : quad1,
176                treat1AsLine ? line1 : line2, lq);
177        if (pts == 2) { // if the line and edge are coincident treat differently
178            _Point midQuad, midLine;
179            double midQuadT = (lq.fT[0][0] + lq.fT[0][1]) / 2;
180            xy_at_t(treat1AsLine ? quad2 : quad1, midQuadT, midQuad.x, midQuad.y);
181            double lineT = t_at(treat1AsLine ? line1 : line2, midQuad);
182            xy_at_t(treat1AsLine ? line1 : line2, lineT, midLine.x, midLine.y);
183            if (AlmostEqualUlps(midQuad.x, midLine.x)
184                    && AlmostEqualUlps(midQuad.y, midLine.y)) {
185                smallT1 = lq.fT[0][0];
186                largeT1 = lq.fT[1][0];
187                smallT2 = lq.fT[0][1];
188                largeT2 = lq.fT[1][1];
189                if (treat2AsLine) {
190                    smallT1 = interp(minT1, maxT1, smallT1);
191                    smallT2 = interp(minT1, maxT1, smallT2);
192                } else {
193                    largeT1 = interp(minT2, maxT2, largeT1);
194                    largeT2 = interp(minT2, maxT2, largeT2);
195                }
196                intersections.addCoincident(smallT1, smallT2, largeT1, largeT2);
197                goto setCoinMinMax;
198            }
199        }
200        for (int index = 0; index < pts; ++index) {
201            smallT1 = lq.fT[0][index];
202            largeT1 = lq.fT[1][index];
203            if (treat2AsLine) {
204                smallT1 = interp(minT1, maxT1, smallT1);
205            } else {
206                largeT1 = interp(minT2, maxT2, largeT1);
207            }
208            intersections.add(smallT1, largeT1);
209        }
210    }
211    if (pts > 0) {
212setCoinMinMax:
213        coinMinT1 = minT1;
214        coinMaxT1 = maxT1;
215        coinMinT2 = minT2;
216        coinMaxT2 = maxT2;
217    }
218    return pts > 0;
219}
220
221bool chop(double minT1, double maxT1, double minT2, double maxT2, int split) {
222    ++depth;
223    intersections.swap();
224    if (split) {
225        ++splits;
226        if (split & 2) {
227            double middle1 = (maxT1 + minT1) / 2;
228            intersect(minT1, middle1, minT2, maxT2);
229            intersect(middle1, maxT1, minT2, maxT2);
230        } else {
231            double middle2 = (maxT2 + minT2) / 2;
232            intersect(minT1, maxT1, minT2, middle2);
233            intersect(minT1, maxT1, middle2, maxT2);
234        }
235        --splits;
236        intersections.swap();
237        --depth;
238        return intersections.intersected();
239    }
240    bool result = intersect(minT1, maxT1, minT2, maxT2);
241    intersections.swap();
242    --depth;
243    return result;
244}
245
246private:
247
248const Quadratic& quad1;
249const Quadratic& quad2;
250Intersections& intersections;
251int depth;
252int splits;
253double quad1Divisions; // line segments to approximate original within error
254double quad2Divisions;
255double coinMinT1; // range of Ts where approximate line intersected curve
256double coinMaxT1;
257double coinMinT2;
258double coinMaxT2;
259};
260
261#include "LineParameters.h"
262
263static void hackToFixPartialCoincidence(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
264    // look to see if non-coincident data basically has unsortable tangents
265
266    // look to see if a point between non-coincident data is on the curve
267    int cIndex;
268    for (int uIndex = 0; uIndex < i.fUsed; ) {
269        double bestDist1 = 1;
270        double bestDist2 = 1;
271        int closest1 = -1;
272        int closest2 = -1;
273        for (cIndex = 0; cIndex < i.fCoincidentUsed; ++cIndex) {
274            double dist = fabs(i.fT[0][uIndex] - i.fCoincidentT[0][cIndex]);
275            if (bestDist1 > dist) {
276                bestDist1 = dist;
277                closest1 = cIndex;
278            }
279            dist = fabs(i.fT[1][uIndex] - i.fCoincidentT[1][cIndex]);
280            if (bestDist2 > dist) {
281                bestDist2 = dist;
282                closest2 = cIndex;
283            }
284        }
285        _Line ends;
286        _Point mid;
287        double t1 = i.fT[0][uIndex];
288        xy_at_t(q1, t1, ends[0].x, ends[0].y);
289        xy_at_t(q1, i.fCoincidentT[0][closest1], ends[1].x, ends[1].y);
290        double midT = (t1 + i.fCoincidentT[0][closest1]) / 2;
291        xy_at_t(q1, midT, mid.x, mid.y);
292        LineParameters params;
293        params.lineEndPoints(ends);
294        double midDist = params.pointDistance(mid);
295        // Note that we prefer to always measure t error, which does not scale,
296        // instead of point error, which is scale dependent. FIXME
297        if (!approximately_zero(midDist)) {
298            ++uIndex;
299            continue;
300        }
301        double t2 = i.fT[1][uIndex];
302        xy_at_t(q2, t2, ends[0].x, ends[0].y);
303        xy_at_t(q2, i.fCoincidentT[1][closest2], ends[1].x, ends[1].y);
304        midT = (t2 + i.fCoincidentT[1][closest2]) / 2;
305        xy_at_t(q2, midT, mid.x, mid.y);
306        params.lineEndPoints(ends);
307        midDist = params.pointDistance(mid);
308        if (!approximately_zero(midDist)) {
309            ++uIndex;
310            continue;
311        }
312        // if both midpoints are close to the line, lengthen coincident span
313        int cEnd = closest1 ^ 1; // assume coincidence always travels in pairs
314        if (!between(i.fCoincidentT[0][cEnd], t1, i.fCoincidentT[0][closest1])) {
315            i.fCoincidentT[0][closest1] = t1;
316        }
317        cEnd = closest2 ^ 1;
318        if (!between(i.fCoincidentT[0][cEnd], t2, i.fCoincidentT[0][closest2])) {
319            i.fCoincidentT[0][closest2] = t2;
320        }
321        int remaining = --i.fUsed - uIndex;
322        if (remaining > 0) {
323            memmove(&i.fT[0][uIndex], &i.fT[0][uIndex + 1], sizeof(i.fT[0][0]) * remaining);
324            memmove(&i.fT[1][uIndex], &i.fT[1][uIndex + 1], sizeof(i.fT[1][0]) * remaining);
325        }
326    }
327    // if coincident data is subjectively a tiny span, replace it with a single point
328    for (cIndex = 0; cIndex < i.fCoincidentUsed; ) {
329        double start1 = i.fCoincidentT[0][cIndex];
330        double end1 = i.fCoincidentT[0][cIndex + 1];
331        _Line ends1;
332        xy_at_t(q1, start1, ends1[0].x, ends1[0].y);
333        xy_at_t(q1, end1, ends1[1].x, ends1[1].y);
334        if (!AlmostEqualUlps(ends1[0].x, ends1[1].x) || AlmostEqualUlps(ends1[0].y, ends1[1].y)) {
335            cIndex += 2;
336            continue;
337        }
338        double start2 = i.fCoincidentT[1][cIndex];
339        double end2 = i.fCoincidentT[1][cIndex + 1];
340        _Line ends2;
341        xy_at_t(q2, start2, ends2[0].x, ends2[0].y);
342        xy_at_t(q2, end2, ends2[1].x, ends2[1].y);
343        // again, approximately should be used with T values, not points FIXME
344        if (!AlmostEqualUlps(ends2[0].x, ends2[1].x) || AlmostEqualUlps(ends2[0].y, ends2[1].y)) {
345            cIndex += 2;
346            continue;
347        }
348        if (approximately_less_than_zero(start1) || approximately_less_than_zero(end1)) {
349            start1 = 0;
350        } else if (approximately_greater_than_one(start1) || approximately_greater_than_one(end1)) {
351            start1 = 1;
352        } else {
353            start1 = (start1 + end1) / 2;
354        }
355        if (approximately_less_than_zero(start2) || approximately_less_than_zero(end2)) {
356            start2 = 0;
357        } else if (approximately_greater_than_one(start2) || approximately_greater_than_one(end2)) {
358            start2 = 1;
359        } else {
360            start2 = (start2 + end2) / 2;
361        }
362        i.insert(start1, start2);
363        i.fCoincidentUsed -= 2;
364        int remaining = i.fCoincidentUsed - cIndex;
365        if (remaining > 0) {
366            memmove(&i.fCoincidentT[0][cIndex], &i.fCoincidentT[0][cIndex + 2], sizeof(i.fCoincidentT[0][0]) * remaining);
367            memmove(&i.fCoincidentT[1][cIndex], &i.fCoincidentT[1][cIndex + 2], sizeof(i.fCoincidentT[1][0]) * remaining);
368        }
369    }
370}
371
372bool intersect(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
373    if (implicit_matches(q1, q2)) {
374        // FIXME: compute T values
375        // compute the intersections of the ends to find the coincident span
376        bool useVertical = fabs(q1[0].x - q1[2].x) < fabs(q1[0].y - q1[2].y);
377        double t;
378        if ((t = axialIntersect(q1, q2[0], useVertical)) >= 0) {
379            i.addCoincident(t, 0);
380        }
381        if ((t = axialIntersect(q1, q2[2], useVertical)) >= 0) {
382            i.addCoincident(t, 1);
383        }
384        useVertical = fabs(q2[0].x - q2[2].x) < fabs(q2[0].y - q2[2].y);
385        if ((t = axialIntersect(q2, q1[0], useVertical)) >= 0) {
386            i.addCoincident(0, t);
387        }
388        if ((t = axialIntersect(q2, q1[2], useVertical)) >= 0) {
389            i.addCoincident(1, t);
390        }
391        SkASSERT(i.fCoincidentUsed <= 2);
392        return i.fCoincidentUsed > 0;
393    }
394    QuadraticIntersections q(q1, q2, i);
395    bool result = q.intersect();
396    // FIXME: partial coincidence detection is currently poor. For now, try
397    // to fix up the data after the fact. In the future, revisit the error
398    // term to try to avoid this kind of result in the first place.
399    if (i.fUsed && i.fCoincidentUsed) {
400        hackToFixPartialCoincidence(q1, q2, i);
401    }
402    return result;
403}
404