Line data Source code
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
14 0 : static bool pointInTriangle(const SkDPoint fPts[3], const SkDPoint& test) {
15 0 : SkDVector v0 = fPts[2] - fPts[0];
16 0 : SkDVector v1 = fPts[1] - fPts[0];
17 0 : SkDVector v2 = test - fPts[0];
18 0 : double dot00 = v0.dot(v0);
19 0 : double dot01 = v0.dot(v1);
20 0 : double dot02 = v0.dot(v2);
21 0 : double dot11 = v1.dot(v1);
22 0 : double dot12 = v1.dot(v2);
23 : // Compute barycentric coordinates
24 0 : double invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
25 0 : double u = (dot11 * dot02 - dot01 * dot12) * invDenom;
26 0 : double v = (dot00 * dot12 - dot01 * dot02) * invDenom;
27 : // Check if point is in triangle
28 0 : return u >= 0 && v >= 0 && u + v < 1;
29 : }
30 :
31 0 : static bool matchesEnd(const SkDPoint fPts[3], const SkDPoint& test) {
32 0 : 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 : */
43 0 : bool SkDQuad::hullIntersects(const SkDQuad& q2, bool* isLinear) const {
44 0 : bool linear = true;
45 0 : for (int oddMan = 0; oddMan < kPointCount; ++oddMan) {
46 : const SkDPoint* endPt[2];
47 0 : this->otherPts(oddMan, endPt);
48 0 : double origX = endPt[0]->fX;
49 0 : double origY = endPt[0]->fY;
50 0 : double adj = endPt[1]->fX - origX;
51 0 : double opp = endPt[1]->fY - origY;
52 0 : double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp;
53 0 : if (approximately_zero(sign)) {
54 0 : continue;
55 : }
56 0 : linear = false;
57 0 : bool foundOutlier = false;
58 0 : for (int n = 0; n < kPointCount; ++n) {
59 0 : double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp;
60 0 : if (test * sign > 0 && !precisely_zero(test)) {
61 0 : foundOutlier = true;
62 0 : break;
63 : }
64 : }
65 0 : if (!foundOutlier) {
66 0 : return false;
67 : }
68 : }
69 0 : 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 0 : if (pointInTriangle(fPts, q2.fPts[0]) || pointInTriangle(fPts, q2.fPts[2])) {
74 0 : linear = false;
75 : }
76 : }
77 0 : *isLinear = linear;
78 0 : return true;
79 : }
80 :
81 0 : bool SkDQuad::hullIntersects(const SkDConic& conic, bool* isLinear) const {
82 0 : return conic.hullIntersects(*this, isLinear);
83 : }
84 :
85 0 : bool SkDQuad::hullIntersects(const SkDCubic& cubic, bool* isLinear) const {
86 0 : 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)
90 : oddMan 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 : */
98 0 : void SkDQuad::otherPts(int oddMan, const SkDPoint* endPt[2]) const {
99 0 : for (int opp = 1; opp < kPointCount; ++opp) {
100 0 : int end = (oddMan ^ opp) - oddMan; // choose a value not equal to oddMan
101 0 : end &= ~(end >> 2); // if the value went negative, set it to zero
102 0 : endPt[opp - 1] = &fPts[end];
103 : }
104 0 : }
105 :
106 33 : int SkDQuad::AddValidTs(double s[], int realRoots, double* t) {
107 33 : int foundRoots = 0;
108 130 : for (int index = 0; index < realRoots; ++index) {
109 97 : double tValue = s[index];
110 97 : if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
111 33 : if (approximately_less_than_zero(tValue)) {
112 3 : tValue = 0;
113 30 : } else if (approximately_greater_than_one(tValue)) {
114 0 : tValue = 1;
115 : }
116 33 : for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
117 0 : if (approximately_equal(t[idx2], tValue)) {
118 0 : goto nextRoot;
119 : }
120 : }
121 33 : t[foundRoots++] = tValue;
122 : }
123 : nextRoot:
124 : {}
125 : }
126 33 : 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
134 0 : int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) {
135 : double s[2];
136 0 : int realRoots = RootsReal(A, B, C, s);
137 0 : int foundRoots = AddValidTs(s, realRoots, t);
138 0 : return foundRoots;
139 : }
140 :
141 : /*
142 : Numeric Solutions (5.6) suggests to solve the quadratic by computing
143 : Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
144 : and using the roots
145 : t1 = Q / A
146 : t2 = C / Q
147 : */
148 : // this does not discard real roots <= 0 or >= 1
149 1 : int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) {
150 1 : const double p = B / (2 * A);
151 1 : const double q = C / A;
152 1 : if (!A || (approximately_zero(A) && (approximately_zero_inverse(p)
153 0 : || approximately_zero_inverse(q)))) {
154 0 : if (approximately_zero(B)) {
155 0 : s[0] = 0;
156 0 : return C == 0;
157 : }
158 0 : s[0] = -C / B;
159 0 : return 1;
160 : }
161 : /* normal form: x^2 + px + q = 0 */
162 1 : const double p2 = p * p;
163 1 : if (!AlmostDequalUlps(p2, q) && p2 < q) {
164 0 : return 0;
165 : }
166 1 : double sqrt_D = 0;
167 1 : if (p2 > q) {
168 1 : sqrt_D = sqrt(p2 - q);
169 : }
170 1 : s[0] = sqrt_D - p;
171 1 : s[1] = -sqrt_D - p;
172 1 : return 1 + !AlmostDequalUlps(s[0], s[1]);
173 : }
174 :
175 0 : bool SkDQuad::isLinear(int startIndex, int endIndex) const {
176 : SkLineParameters lineParameters;
177 0 : lineParameters.quadEndPoints(*this, startIndex, endIndex);
178 : // FIXME: maybe it's possible to avoid this and compare non-normalized
179 0 : lineParameters.normalize();
180 0 : double distance = lineParameters.controlPtDistance(*this);
181 : double tiniest = SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(fPts[0].fX, fPts[0].fY),
182 0 : fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY);
183 0 : double largest = SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(fPts[0].fX, fPts[0].fY),
184 0 : fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY);
185 0 : largest = SkTMax(largest, -tiniest);
186 0 : return approximately_zero_when_compared_to(distance, largest);
187 : }
188 :
189 0 : SkDVector SkDQuad::dxdyAtT(double t) const {
190 0 : double a = t - 1;
191 0 : double b = 1 - 2 * t;
192 0 : double c = t;
193 0 : SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
194 0 : a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
195 0 : if (result.fX == 0 && result.fY == 0) {
196 0 : if (zero_or_one(t)) {
197 0 : result = fPts[2] - fPts[0];
198 : } else {
199 : // incomplete
200 0 : SkDebugf("!q");
201 : }
202 : }
203 0 : return result;
204 : }
205 :
206 : // OPTIMIZE: assert if caller passes in t == 0 / t == 1 ?
207 0 : SkDPoint SkDQuad::ptAtT(double t) const {
208 0 : if (0 == t) {
209 0 : return fPts[0];
210 : }
211 0 : if (1 == t) {
212 0 : return fPts[2];
213 : }
214 0 : double one_t = 1 - t;
215 0 : double a = one_t * one_t;
216 0 : double b = 2 * one_t * t;
217 0 : double c = t * t;
218 0 : SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
219 0 : a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
220 0 : return result;
221 : }
222 :
223 0 : static double interp_quad_coords(const double* src, double t) {
224 0 : if (0 == t) {
225 0 : return src[0];
226 : }
227 0 : if (1 == t) {
228 0 : return src[4];
229 : }
230 0 : double ab = SkDInterp(src[0], src[2], t);
231 0 : double bc = SkDInterp(src[2], src[4], t);
232 0 : double abc = SkDInterp(ab, bc, t);
233 0 : return abc;
234 : }
235 :
236 0 : bool SkDQuad::monotonicInX() const {
237 0 : return between(fPts[0].fX, fPts[1].fX, fPts[2].fX);
238 : }
239 :
240 0 : bool SkDQuad::monotonicInY() const {
241 0 : return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
242 : }
243 :
244 : /*
245 : Given a quadratic q, t1, and t2, find a small quadratic segment.
246 :
247 : The 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 :
251 : To find B, compute the point halfway between t1 and t2:
252 :
253 : q(at (t1 + t2)/2) == D
254 :
255 : Next, compute where D must be if we know the value of B:
256 :
257 : _12 = A/2 + B/2
258 : 12_ = B/2 + C/2
259 : 123 = A/4 + B/2 + C/4
260 : = D
261 :
262 : Group the known values on one side:
263 :
264 : B = D*2 - A/2 - C/2
265 : */
266 :
267 : // OPTIMIZE? : special case t1 = 1 && t2 = 0
268 0 : SkDQuad SkDQuad::subDivide(double t1, double t2) const {
269 0 : if (0 == t1 && 1 == t2) {
270 0 : return *this;
271 : }
272 : SkDQuad dst;
273 0 : double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
274 0 : double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
275 0 : double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
276 0 : double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
277 0 : double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
278 0 : double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
279 0 : /* bx = */ dst[1].fX = 2 * dx - (ax + cx) / 2;
280 0 : /* by = */ dst[1].fY = 2 * dy - (ay + cy) / 2;
281 0 : return dst;
282 : }
283 :
284 0 : void SkDQuad::align(int endIndex, SkDPoint* dstPt) const {
285 0 : if (fPts[endIndex].fX == fPts[1].fX) {
286 0 : dstPt->fX = fPts[endIndex].fX;
287 : }
288 0 : if (fPts[endIndex].fY == fPts[1].fY) {
289 0 : dstPt->fY = fPts[endIndex].fY;
290 : }
291 0 : }
292 :
293 0 : SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
294 0 : SkASSERT(t1 != t2);
295 : SkDPoint b;
296 0 : SkDQuad sub = subDivide(t1, t2);
297 0 : SkDLine b0 = {{a, sub[1] + (a - sub[0])}};
298 0 : SkDLine b1 = {{c, sub[1] + (c - sub[2])}};
299 0 : SkIntersections i;
300 0 : i.intersectRay(b0, b1);
301 0 : if (i.used() == 1 && i[0][0] >= 0 && i[1][0] >= 0) {
302 0 : b = i.pt(0);
303 : } else {
304 0 : SkASSERT(i.used() <= 2);
305 0 : return SkDPoint::Mid(b0[1], b1[1]);
306 : }
307 0 : if (t1 == 0 || t2 == 0) {
308 0 : align(0, &b);
309 : }
310 0 : if (t1 == 1 || t2 == 1) {
311 0 : align(2, &b);
312 : }
313 0 : if (AlmostBequalUlps(b.fX, a.fX)) {
314 0 : b.fX = a.fX;
315 0 : } else if (AlmostBequalUlps(b.fX, c.fX)) {
316 0 : b.fX = c.fX;
317 : }
318 0 : if (AlmostBequalUlps(b.fY, a.fY)) {
319 0 : b.fY = a.fY;
320 0 : } else if (AlmostBequalUlps(b.fY, c.fY)) {
321 0 : b.fY = c.fY;
322 : }
323 0 : return b;
324 : }
325 :
326 : /* classic one t subdivision */
327 0 : static void interp_quad_coords(const double* src, double* dst, double t) {
328 0 : double ab = SkDInterp(src[0], src[2], t);
329 0 : double bc = SkDInterp(src[2], src[4], t);
330 0 : dst[0] = src[0];
331 0 : dst[2] = ab;
332 0 : dst[4] = SkDInterp(ab, bc, t);
333 0 : dst[6] = bc;
334 0 : dst[8] = src[4];
335 0 : }
336 :
337 0 : SkDQuadPair SkDQuad::chopAt(double t) const
338 : {
339 : SkDQuadPair dst;
340 0 : interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
341 0 : interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
342 0 : return dst;
343 : }
344 :
345 0 : static int valid_unit_divide(double numer, double denom, double* ratio)
346 : {
347 0 : if (numer < 0) {
348 0 : numer = -numer;
349 0 : denom = -denom;
350 : }
351 0 : if (denom == 0 || numer == 0 || numer >= denom) {
352 0 : return 0;
353 : }
354 0 : double r = numer / denom;
355 0 : if (r == 0) { // catch underflow if numer <<<< denom
356 0 : return 0;
357 : }
358 0 : *ratio = r;
359 0 : 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 : */
367 0 : int SkDQuad::FindExtrema(const double src[], double tValue[1]) {
368 : /* At + B == 0
369 : t = -B / A
370 : */
371 0 : double a = src[0];
372 0 : double b = src[2];
373 0 : double c = src[4];
374 0 : 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 : */
383 0 : void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) {
384 0 : *a = quad[0]; // a = A
385 0 : *b = 2 * quad[2]; // b = 2*B
386 0 : *c = quad[4]; // c = C
387 0 : *b -= *c; // b = 2*B - C
388 0 : *a -= *b; // a = A - 2*B + C
389 0 : *b -= *c; // b = 2*B - 2*C
390 0 : }
|
