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