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 "SkLineParameters.h"
9#include "SkPathOpsCubic.h"
10#include "SkPathOpsCurve.h"
11#include "SkPathOpsQuad.h"
12
13// from blackpawn.com/texts/pointinpoly
14static bool pointInTriangle(const SkDPoint fPts[3], const SkDPoint& test) {
15    SkDVector v0 = fPts[2] - fPts[0];
16    SkDVector v1 = fPts[1] - fPts[0];
17    SkDVector v2 = test - fPts[0];
18    double dot00 = v0.dot(v0);
19    double dot01 = v0.dot(v1);
20    double dot02 = v0.dot(v2);
21    double dot11 = v1.dot(v1);
22    double dot12 = v1.dot(v2);
23    // Compute barycentric coordinates
24    double invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
25    double u = (dot11 * dot02 - dot01 * dot12) * invDenom;
26    double v = (dot00 * dot12 - dot01 * dot02) * invDenom;
27    // Check if point is in triangle
28    return u >= 0 && v >= 0 && u + v < 1;
29}
30
31static bool matchesEnd(const SkDPoint fPts[3], const SkDPoint& test) {
32    return fPts[0] == test || fPts[2] == test;
33}
34
35/* started with at_most_end_pts_in_common from SkDQuadIntersection.cpp */
36// Do a quick reject by rotating all points relative to a line formed by
37// a pair of one quad's points. If the 2nd quad's points
38// are on the line or on the opposite side from the 1st quad's 'odd man', the
39// curves at most intersect at the endpoints.
40/* if returning true, check contains true if quad's hull collapsed, making the cubic linear
41   if returning false, check contains true if the the quad pair have only the end point in common
42*/
43bool SkDQuad::hullIntersects(const SkDQuad& q2, bool* isLinear) const {
44    bool linear = true;
45    for (int oddMan = 0; oddMan < kPointCount; ++oddMan) {
46        const SkDPoint* endPt[2];
47        this->otherPts(oddMan, endPt);
48        double origX = endPt[0]->fX;
49        double origY = endPt[0]->fY;
50        double adj = endPt[1]->fX - origX;
51        double opp = endPt[1]->fY - origY;
52        double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp;
53        if (approximately_zero(sign)) {
54            continue;
55        }
56        linear = false;
57        bool foundOutlier = false;
58        for (int n = 0; n < kPointCount; ++n) {
59            double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp;
60            if (test * sign > 0 && !precisely_zero(test)) {
61                foundOutlier = true;
62                break;
63            }
64        }
65        if (!foundOutlier) {
66            return false;
67        }
68    }
69    if (linear && !matchesEnd(fPts, q2.fPts[0]) && !matchesEnd(fPts, q2.fPts[2])) {
70        // if the end point of the opposite quad is inside the hull that is nearly a line,
71        // then representing the quad as a line may cause the intersection to be missed.
72        // Check to see if the endpoint is in the triangle.
73        if (pointInTriangle(fPts, q2.fPts[0]) || pointInTriangle(fPts, q2.fPts[2])) {
74            linear = false;
75        }
76    }
77    *isLinear = linear;
78    return true;
79}
80
81bool SkDQuad::hullIntersects(const SkDConic& conic, bool* isLinear) const {
82    return conic.hullIntersects(*this, isLinear);
83}
84
85bool SkDQuad::hullIntersects(const SkDCubic& cubic, bool* isLinear) const {
86    return cubic.hullIntersects(*this, isLinear);
87}
88
89/* bit twiddling for finding the off curve index (x&~m is the pair in [0,1,2] excluding oddMan)
90oddMan    opp   x=oddMan^opp  x=x-oddMan  m=x>>2   x&~m
91    0       1         1            1         0       1
92            2         2            2         0       2
93    1       1         0           -1        -1       0
94            2         3            2         0       2
95    2       1         3            1         0       1
96            2         0           -2        -1       0
97*/
98void SkDQuad::otherPts(int oddMan, const SkDPoint* endPt[2]) const {
99    for (int opp = 1; opp < kPointCount; ++opp) {
100        int end = (oddMan ^ opp) - oddMan;  // choose a value not equal to oddMan
101        end &= ~(end >> 2);  // if the value went negative, set it to zero
102        endPt[opp - 1] = &fPts[end];
103    }
104}
105
106int SkDQuad::AddValidTs(double s[], int realRoots, double* t) {
107    int foundRoots = 0;
108    for (int index = 0; index < realRoots; ++index) {
109        double tValue = s[index];
110        if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
111            if (approximately_less_than_zero(tValue)) {
112                tValue = 0;
113            } else if (approximately_greater_than_one(tValue)) {
114                tValue = 1;
115            }
116            for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
117                if (approximately_equal(t[idx2], tValue)) {
118                    goto nextRoot;
119                }
120            }
121            t[foundRoots++] = tValue;
122        }
123nextRoot:
124        {}
125    }
126    return foundRoots;
127}
128
129// note: caller expects multiple results to be sorted smaller first
130// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
131//  analysis of the quadratic equation, suggesting why the following looks at
132//  the sign of B -- and further suggesting that the greatest loss of precision
133//  is in b squared less two a c
134int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) {
135    double s[2];
136    int realRoots = RootsReal(A, B, C, s);
137    int foundRoots = AddValidTs(s, realRoots, t);
138    return foundRoots;
139}
140
141/*
142Numeric Solutions (5.6) suggests to solve the quadratic by computing
143       Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
144and using the roots
145      t1 = Q / A
146      t2 = C / Q
147*/
148// this does not discard real roots <= 0 or >= 1
149int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) {
150    const double p = B / (2 * A);
151    const double q = C / A;
152    if (!A || (approximately_zero(A) && (approximately_zero_inverse(p)
153            || approximately_zero_inverse(q)))) {
154        if (approximately_zero(B)) {
155            s[0] = 0;
156            return C == 0;
157        }
158        s[0] = -C / B;
159        return 1;
160    }
161    /* normal form: x^2 + px + q = 0 */
162    const double p2 = p * p;
163    if (!AlmostDequalUlps(p2, q) && p2 < q) {
164        return 0;
165    }
166    double sqrt_D = 0;
167    if (p2 > q) {
168        sqrt_D = sqrt(p2 - q);
169    }
170    s[0] = sqrt_D - p;
171    s[1] = -sqrt_D - p;
172    return 1 + !AlmostDequalUlps(s[0], s[1]);
173}
174
175bool SkDQuad::isLinear(int startIndex, int endIndex) const {
176    SkLineParameters lineParameters;
177    lineParameters.quadEndPoints(*this, startIndex, endIndex);
178    // FIXME: maybe it's possible to avoid this and compare non-normalized
179    lineParameters.normalize();
180    double distance = lineParameters.controlPtDistance(*this);
181    double tiniest = SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(fPts[0].fX, fPts[0].fY),
182            fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY);
183    double largest = SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(fPts[0].fX, fPts[0].fY),
184            fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY);
185    largest = SkTMax(largest, -tiniest);
186    return approximately_zero_when_compared_to(distance, largest);
187}
188
189SkDVector SkDQuad::dxdyAtT(double t) const {
190    double a = t - 1;
191    double b = 1 - 2 * t;
192    double c = t;
193    SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
194            a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
195    if (result.fX == 0 && result.fY == 0) {
196        if (zero_or_one(t)) {
197            result = fPts[2] - fPts[0];
198        } else {
199            // incomplete
200            SkDebugf("!q");
201        }
202    }
203    return result;
204}
205
206// OPTIMIZE: assert if caller passes in t == 0 / t == 1 ?
207SkDPoint SkDQuad::ptAtT(double t) const {
208    if (0 == t) {
209        return fPts[0];
210    }
211    if (1 == t) {
212        return fPts[2];
213    }
214    double one_t = 1 - t;
215    double a = one_t * one_t;
216    double b = 2 * one_t * t;
217    double c = t * t;
218    SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
219            a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
220    return result;
221}
222
223static double interp_quad_coords(const double* src, double t) {
224    if (0 == t) {
225        return src[0];
226    }
227    if (1 == t) {
228        return src[4];
229    }
230    double ab = SkDInterp(src[0], src[2], t);
231    double bc = SkDInterp(src[2], src[4], t);
232    double abc = SkDInterp(ab, bc, t);
233    return abc;
234}
235
236bool SkDQuad::monotonicInX() const {
237    return between(fPts[0].fX, fPts[1].fX, fPts[2].fX);
238}
239
240bool SkDQuad::monotonicInY() const {
241    return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
242}
243
244/*
245Given a quadratic q, t1, and t2, find a small quadratic segment.
246
247The new quadratic is defined by A, B, and C, where
248 A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
249 C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
250
251To find B, compute the point halfway between t1 and t2:
252
253q(at (t1 + t2)/2) == D
254
255Next, compute where D must be if we know the value of B:
256
257_12 = A/2 + B/2
25812_ = B/2 + C/2
259123 = A/4 + B/2 + C/4
260    = D
261
262Group the known values on one side:
263
264B   = D*2 - A/2 - C/2
265*/
266
267// OPTIMIZE? : special case  t1 = 1 && t2 = 0
268SkDQuad SkDQuad::subDivide(double t1, double t2) const {
269    if (0 == t1 && 1 == t2) {
270        return *this;
271    }
272    SkDQuad dst;
273    double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
274    double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
275    double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
276    double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
277    double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
278    double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
279    /* bx = */ dst[1].fX = 2 * dx - (ax + cx) / 2;
280    /* by = */ dst[1].fY = 2 * dy - (ay + cy) / 2;
281    return dst;
282}
283
284void SkDQuad::align(int endIndex, SkDPoint* dstPt) const {
285    if (fPts[endIndex].fX == fPts[1].fX) {
286        dstPt->fX = fPts[endIndex].fX;
287    }
288    if (fPts[endIndex].fY == fPts[1].fY) {
289        dstPt->fY = fPts[endIndex].fY;
290    }
291}
292
293SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
294    SkASSERT(t1 != t2);
295    SkDPoint b;
296    SkDQuad sub = subDivide(t1, t2);
297    SkDLine b0 = {{a, sub[1] + (a - sub[0])}};
298    SkDLine b1 = {{c, sub[1] + (c - sub[2])}};
299    SkIntersections i;
300    i.intersectRay(b0, b1);
301    if (i.used() == 1 && i[0][0] >= 0 && i[1][0] >= 0) {
302        b = i.pt(0);
303    } else {
304        SkASSERT(i.used() <= 2);
305        return SkDPoint::Mid(b0[1], b1[1]);
306    }
307    if (t1 == 0 || t2 == 0) {
308        align(0, &b);
309    }
310    if (t1 == 1 || t2 == 1) {
311        align(2, &b);
312    }
313    if (AlmostBequalUlps(b.fX, a.fX)) {
314        b.fX = a.fX;
315    } else if (AlmostBequalUlps(b.fX, c.fX)) {
316        b.fX = c.fX;
317    }
318    if (AlmostBequalUlps(b.fY, a.fY)) {
319        b.fY = a.fY;
320    } else if (AlmostBequalUlps(b.fY, c.fY)) {
321        b.fY = c.fY;
322    }
323    return b;
324}
325
326/* classic one t subdivision */
327static void interp_quad_coords(const double* src, double* dst, double t) {
328    double ab = SkDInterp(src[0], src[2], t);
329    double bc = SkDInterp(src[2], src[4], t);
330    dst[0] = src[0];
331    dst[2] = ab;
332    dst[4] = SkDInterp(ab, bc, t);
333    dst[6] = bc;
334    dst[8] = src[4];
335}
336
337SkDQuadPair SkDQuad::chopAt(double t) const
338{
339    SkDQuadPair dst;
340    interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
341    interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
342    return dst;
343}
344
345static int valid_unit_divide(double numer, double denom, double* ratio)
346{
347    if (numer < 0) {
348        numer = -numer;
349        denom = -denom;
350    }
351    if (denom == 0 || numer == 0 || numer >= denom) {
352        return 0;
353    }
354    double r = numer / denom;
355    if (r == 0) {  // catch underflow if numer <<<< denom
356        return 0;
357    }
358    *ratio = r;
359    return 1;
360}
361
362/** Quad'(t) = At + B, where
363    A = 2(a - 2b + c)
364    B = 2(b - a)
365    Solve for t, only if it fits between 0 < t < 1
366*/
367int SkDQuad::FindExtrema(const double src[], double tValue[1]) {
368    /*  At + B == 0
369        t = -B / A
370    */
371    double a = src[0];
372    double b = src[2];
373    double c = src[4];
374    return valid_unit_divide(a - b, a - b - b + c, tValue);
375}
376
377/* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
378 *
379 * a = A - 2*B +   C
380 * b =     2*B - 2*C
381 * c =             C
382 */
383void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) {
384    *a = quad[0];      // a = A
385    *b = 2 * quad[2];  // b =     2*B
386    *c = quad[4];      // c =             C
387    *b -= *c;          // b =     2*B -   C
388    *a -= *b;          // a = A - 2*B +   C
389    *b -= *c;          // b =     2*B - 2*C
390}
391