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 (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) {
123 if (approximately_zero(B)) {
124 s[0] = 0;
125 return C == 0;
126 }
127 s[0] = -C / B;
128 return 1;
129 }
130 /* normal form: x^2 + px + q = 0 */
131 const double p2 = p * p;
132 if (!AlmostDequalUlps(p2, q) && p2 < q) {
133 return 0;
134 }
135 double sqrt_D = 0;
136 if (p2 > q) {
137 sqrt_D = sqrt(p2 - q);
138 }
139 s[0] = sqrt_D - p;
140 s[1] = -sqrt_D - p;
141 return 1 + !AlmostDequalUlps(s[0], s[1]);
142 }
143
isLinear(int startIndex,int endIndex) const144 bool SkDQuad::isLinear(int startIndex, int endIndex) const {
145 SkLineParameters lineParameters;
146 lineParameters.quadEndPoints(*this, startIndex, endIndex);
147 // FIXME: maybe it's possible to avoid this and compare non-normalized
148 lineParameters.normalize();
149 double distance = lineParameters.controlPtDistance(*this);
150 double tiniest = SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(fPts[0].fX, fPts[0].fY),
151 fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY);
152 double largest = SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(fPts[0].fX, fPts[0].fY),
153 fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY);
154 largest = SkTMax(largest, -tiniest);
155 return approximately_zero_when_compared_to(distance, largest);
156 }
157
dxdyAtT(double t) const158 SkDVector SkDQuad::dxdyAtT(double t) const {
159 double a = t - 1;
160 double b = 1 - 2 * t;
161 double c = t;
162 SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
163 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
164 return result;
165 }
166
167 // OPTIMIZE: assert if caller passes in t == 0 / t == 1 ?
ptAtT(double t) const168 SkDPoint SkDQuad::ptAtT(double t) const {
169 if (0 == t) {
170 return fPts[0];
171 }
172 if (1 == t) {
173 return fPts[2];
174 }
175 double one_t = 1 - t;
176 double a = one_t * one_t;
177 double b = 2 * one_t * t;
178 double c = t * t;
179 SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
180 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
181 return result;
182 }
183
interp_quad_coords(const double * src,double t)184 static double interp_quad_coords(const double* src, double t) {
185 double ab = SkDInterp(src[0], src[2], t);
186 double bc = SkDInterp(src[2], src[4], t);
187 double abc = SkDInterp(ab, bc, t);
188 return abc;
189 }
190
monotonicInX() const191 bool SkDQuad::monotonicInX() const {
192 return between(fPts[0].fX, fPts[1].fX, fPts[2].fX);
193 }
194
monotonicInY() const195 bool SkDQuad::monotonicInY() const {
196 return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
197 }
198
199 /*
200 Given a quadratic q, t1, and t2, find a small quadratic segment.
201
202 The new quadratic is defined by A, B, and C, where
203 A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
204 C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
205
206 To find B, compute the point halfway between t1 and t2:
207
208 q(at (t1 + t2)/2) == D
209
210 Next, compute where D must be if we know the value of B:
211
212 _12 = A/2 + B/2
213 12_ = B/2 + C/2
214 123 = A/4 + B/2 + C/4
215 = D
216
217 Group the known values on one side:
218
219 B = D*2 - A/2 - C/2
220 */
221
222 // OPTIMIZE : special case either or both of t1 = 0, t2 = 1
subDivide(double t1,double t2) const223 SkDQuad SkDQuad::subDivide(double t1, double t2) const {
224 SkDQuad dst;
225 double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
226 double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
227 double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
228 double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
229 double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
230 double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
231 /* bx = */ dst[1].fX = 2 * dx - (ax + cx) / 2;
232 /* by = */ dst[1].fY = 2 * dy - (ay + cy) / 2;
233 return dst;
234 }
235
align(int endIndex,SkDPoint * dstPt) const236 void SkDQuad::align(int endIndex, SkDPoint* dstPt) const {
237 if (fPts[endIndex].fX == fPts[1].fX) {
238 dstPt->fX = fPts[endIndex].fX;
239 }
240 if (fPts[endIndex].fY == fPts[1].fY) {
241 dstPt->fY = fPts[endIndex].fY;
242 }
243 }
244
subDivide(const SkDPoint & a,const SkDPoint & c,double t1,double t2) const245 SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
246 SkASSERT(t1 != t2);
247 SkDPoint b;
248 SkDQuad sub = subDivide(t1, t2);
249 SkDLine b0 = {{a, sub[1] + (a - sub[0])}};
250 SkDLine b1 = {{c, sub[1] + (c - sub[2])}};
251 SkIntersections i;
252 i.intersectRay(b0, b1);
253 if (i.used() == 1 && i[0][0] >= 0 && i[1][0] >= 0) {
254 b = i.pt(0);
255 } else {
256 SkASSERT(i.used() <= 2);
257 b = SkDPoint::Mid(b0[1], b1[1]);
258 }
259 if (t1 == 0 || t2 == 0) {
260 align(0, &b);
261 }
262 if (t1 == 1 || t2 == 1) {
263 align(2, &b);
264 }
265 if (AlmostBequalUlps(b.fX, a.fX)) {
266 b.fX = a.fX;
267 } else if (AlmostBequalUlps(b.fX, c.fX)) {
268 b.fX = c.fX;
269 }
270 if (AlmostBequalUlps(b.fY, a.fY)) {
271 b.fY = a.fY;
272 } else if (AlmostBequalUlps(b.fY, c.fY)) {
273 b.fY = c.fY;
274 }
275 return b;
276 }
277
278 /* classic one t subdivision */
interp_quad_coords(const double * src,double * dst,double t)279 static void interp_quad_coords(const double* src, double* dst, double t) {
280 double ab = SkDInterp(src[0], src[2], t);
281 double bc = SkDInterp(src[2], src[4], t);
282 dst[0] = src[0];
283 dst[2] = ab;
284 dst[4] = SkDInterp(ab, bc, t);
285 dst[6] = bc;
286 dst[8] = src[4];
287 }
288
chopAt(double t) const289 SkDQuadPair SkDQuad::chopAt(double t) const
290 {
291 SkDQuadPair dst;
292 interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
293 interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
294 return dst;
295 }
296
valid_unit_divide(double numer,double denom,double * ratio)297 static int valid_unit_divide(double numer, double denom, double* ratio)
298 {
299 if (numer < 0) {
300 numer = -numer;
301 denom = -denom;
302 }
303 if (denom == 0 || numer == 0 || numer >= denom) {
304 return 0;
305 }
306 double r = numer / denom;
307 if (r == 0) { // catch underflow if numer <<<< denom
308 return 0;
309 }
310 *ratio = r;
311 return 1;
312 }
313
314 /** Quad'(t) = At + B, where
315 A = 2(a - 2b + c)
316 B = 2(b - a)
317 Solve for t, only if it fits between 0 < t < 1
318 */
FindExtrema(const double src[],double tValue[1])319 int SkDQuad::FindExtrema(const double src[], double tValue[1]) {
320 /* At + B == 0
321 t = -B / A
322 */
323 double a = src[0];
324 double b = src[2];
325 double c = src[4];
326 return valid_unit_divide(a - b, a - b - b + c, tValue);
327 }
328
329 /* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
330 *
331 * a = A - 2*B + C
332 * b = 2*B - 2*C
333 * c = C
334 */
SetABC(const double * quad,double * a,double * b,double * c)335 void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) {
336 *a = quad[0]; // a = A
337 *b = 2 * quad[2]; // b = 2*B
338 *c = quad[4]; // c = C
339 *b -= *c; // b = 2*B - C
340 *a -= *b; // a = A - 2*B + C
341 *b -= *c; // b = 2*B - 2*C
342 }
343