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 "src/pathops/SkIntersections.h"
8#include "src/pathops/SkLineParameters.h"
9#include "src/pathops/SkPathOpsCubic.h"
10#include "src/pathops/SkPathOpsCurve.h"
11#include "src/pathops/SkPathOpsQuad.h"
12#include "src/pathops/SkPathOpsRect.h"
13
14// from blackpawn.com/texts/pointinpoly
15static 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
35static 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*/
47bool 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
85bool SkDQuad::hullIntersects(const SkDConic& conic, bool* isLinear) const {
86 return conic.hullIntersects(*this, isLinear);
87}
88
89bool 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)
94oddMan 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*/
102void 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
110int 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 }
127nextRoot:
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
138int 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
145static 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/*
155Numeric Solutions (5.6) suggests to solve the quadratic by computing
156 Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
157and using the roots
158 t1 = Q / A
159 t2 = C / Q
160*/
161// this does not discard real roots <= 0 or >= 1
162int 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
185bool 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 = std::min(std::min(std::min(std::min(std::min(fPts[0].fX, fPts[0].fY),
192 fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY);
193 double largest = std::max(std::max(std::max(std::max(std::max(fPts[0].fX, fPts[0].fY),
194 fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY);
195 largest = std::max(largest, -tiniest);
196 return approximately_zero_when_compared_to(distance, largest);
197}
198
199SkDVector 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 ?
217SkDPoint 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
233static 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
246bool SkDQuad::monotonicInX() const {
247 return between(fPts[0].fX, fPts[1].fX, fPts[2].fX);
248}
249
250bool SkDQuad::monotonicInY() const {
251 return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
252}
253
254/*
255Given a quadratic q, t1, and t2, find a small quadratic segment.
256
257The 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
261To find B, compute the point halfway between t1 and t2:
262
263q(at (t1 + t2)/2) == D
264
265Next, compute where D must be if we know the value of B:
266
267_12 = A/2 + B/2
26812_ = B/2 + C/2
269123 = A/4 + B/2 + C/4
270 = D
271
272Group the known values on one side:
273
274B = D*2 - A/2 - C/2
275*/
276
277// OPTIMIZE? : special case t1 = 1 && t2 = 0
278SkDQuad 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
294void 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
303SkDPoint 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 */
337static 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
347SkDQuadPair 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
355static 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*/
377int 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 */
393void 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
402int SkTQuad::intersectRay(SkIntersections* i, const SkDLine& line) const {
403 return i->intersectRay(fQuad, line);
404}
405
406bool SkTQuad::hullIntersects(const SkDConic& conic, bool* isLinear) const {
407 return conic.hullIntersects(fQuad, isLinear);
408}
409
410bool SkTQuad::hullIntersects(const SkDCubic& cubic, bool* isLinear) const {
411 return cubic.hullIntersects(fQuad, isLinear);
412}
413
414void SkTQuad::setBounds(SkDRect* rect) const {
415 rect->setBounds(fQuad);
416}
417