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