blob: c1d068af345701e64c5251007d172d29e83fbeb3 [file] [log] [blame]
caryclark@google.com07393ca2013-04-08 11:47:37 +00001/*
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 */
caryclark@google.comcffbcc32013-06-04 17:59:42 +00007#include "SkIntersections.h"
caryclark@google.com07393ca2013-04-08 11:47:37 +00008#include "SkLineParameters.h"
9#include "SkPathOpsCubic.h"
10#include "SkPathOpsQuad.h"
11#include "SkPathOpsTriangle.h"
12
13// from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
14// (currently only used by testing)
15double SkDQuad::nearestT(const SkDPoint& pt) const {
16 SkDVector pos = fPts[0] - pt;
17 // search points P of bezier curve with PM.(dP / dt) = 0
18 // a calculus leads to a 3d degree equation :
19 SkDVector A = fPts[1] - fPts[0];
20 SkDVector B = fPts[2] - fPts[1];
21 B -= A;
22 double a = B.dot(B);
23 double b = 3 * A.dot(B);
24 double c = 2 * A.dot(A) + pos.dot(B);
25 double d = pos.dot(A);
26 double ts[3];
27 int roots = SkDCubic::RootsValidT(a, b, c, d, ts);
28 double d0 = pt.distanceSquared(fPts[0]);
29 double d2 = pt.distanceSquared(fPts[2]);
caryclark@google.com3b97af52013-04-23 11:56:44 +000030 double distMin = SkTMin(d0, d2);
caryclark@google.com07393ca2013-04-08 11:47:37 +000031 int bestIndex = -1;
32 for (int index = 0; index < roots; ++index) {
caryclark@google.com4fdbb222013-07-23 15:27:41 +000033 SkDPoint onQuad = ptAtT(ts[index]);
caryclark@google.com07393ca2013-04-08 11:47:37 +000034 double dist = pt.distanceSquared(onQuad);
35 if (distMin > dist) {
36 distMin = dist;
37 bestIndex = index;
38 }
39 }
40 if (bestIndex >= 0) {
41 return ts[bestIndex];
42 }
43 return d0 < d2 ? 0 : 1;
44}
45
46bool SkDQuad::pointInHull(const SkDPoint& pt) const {
47 return ((const SkDTriangle&) fPts).contains(pt);
48}
49
50SkDPoint SkDQuad::top(double startT, double endT) const {
51 SkDQuad sub = subDivide(startT, endT);
52 SkDPoint topPt = sub[0];
53 if (topPt.fY > sub[2].fY || (topPt.fY == sub[2].fY && topPt.fX > sub[2].fX)) {
54 topPt = sub[2];
55 }
56 if (!between(sub[0].fY, sub[1].fY, sub[2].fY)) {
57 double extremeT;
58 if (FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, &extremeT)) {
59 extremeT = startT + (endT - startT) * extremeT;
caryclark@google.com4fdbb222013-07-23 15:27:41 +000060 SkDPoint test = ptAtT(extremeT);
caryclark@google.com07393ca2013-04-08 11:47:37 +000061 if (topPt.fY > test.fY || (topPt.fY == test.fY && topPt.fX > test.fX)) {
62 topPt = test;
63 }
64 }
65 }
66 return topPt;
67}
68
69int SkDQuad::AddValidTs(double s[], int realRoots, double* t) {
70 int foundRoots = 0;
71 for (int index = 0; index < realRoots; ++index) {
72 double tValue = s[index];
73 if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
74 if (approximately_less_than_zero(tValue)) {
75 tValue = 0;
76 } else if (approximately_greater_than_one(tValue)) {
77 tValue = 1;
78 }
79 for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
80 if (approximately_equal(t[idx2], tValue)) {
81 goto nextRoot;
82 }
83 }
84 t[foundRoots++] = tValue;
85 }
86nextRoot:
87 {}
88 }
89 return foundRoots;
90}
91
92// note: caller expects multiple results to be sorted smaller first
93// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
94// analysis of the quadratic equation, suggesting why the following looks at
95// the sign of B -- and further suggesting that the greatest loss of precision
96// is in b squared less two a c
97int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) {
98 double s[2];
99 int realRoots = RootsReal(A, B, C, s);
100 int foundRoots = AddValidTs(s, realRoots, t);
101 return foundRoots;
102}
103
104/*
105Numeric Solutions (5.6) suggests to solve the quadratic by computing
106 Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
107and using the roots
108 t1 = Q / A
109 t2 = C / Q
110*/
111// this does not discard real roots <= 0 or >= 1
112int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) {
113 const double p = B / (2 * A);
114 const double q = C / A;
115 if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) {
116 if (approximately_zero(B)) {
117 s[0] = 0;
118 return C == 0;
119 }
120 s[0] = -C / B;
121 return 1;
122 }
123 /* normal form: x^2 + px + q = 0 */
124 const double p2 = p * p;
caryclark@google.com7eaa53d2013-10-02 14:49:34 +0000125 if (!AlmostDequalUlps(p2, q) && p2 < q) {
caryclark@google.com07393ca2013-04-08 11:47:37 +0000126 return 0;
127 }
128 double sqrt_D = 0;
129 if (p2 > q) {
130 sqrt_D = sqrt(p2 - q);
131 }
132 s[0] = sqrt_D - p;
133 s[1] = -sqrt_D - p;
caryclark@google.com7eaa53d2013-10-02 14:49:34 +0000134 return 1 + !AlmostDequalUlps(s[0], s[1]);
caryclark@google.com07393ca2013-04-08 11:47:37 +0000135}
136
137bool SkDQuad::isLinear(int startIndex, int endIndex) const {
138 SkLineParameters lineParameters;
139 lineParameters.quadEndPoints(*this, startIndex, endIndex);
140 // FIXME: maybe it's possible to avoid this and compare non-normalized
141 lineParameters.normalize();
142 double distance = lineParameters.controlPtDistance(*this);
143 return approximately_zero(distance);
144}
145
146SkDCubic SkDQuad::toCubic() const {
147 SkDCubic cubic;
148 cubic[0] = fPts[0];
149 cubic[2] = fPts[1];
150 cubic[3] = fPts[2];
151 cubic[1].fX = (cubic[0].fX + cubic[2].fX * 2) / 3;
152 cubic[1].fY = (cubic[0].fY + cubic[2].fY * 2) / 3;
153 cubic[2].fX = (cubic[3].fX + cubic[2].fX * 2) / 3;
154 cubic[2].fY = (cubic[3].fY + cubic[2].fY * 2) / 3;
155 return cubic;
156}
157
158SkDVector 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
caryclark@google.comfa2aeee2013-07-15 13:29:13 +0000167// OPTIMIZE: assert if caller passes in t == 0 / t == 1 ?
caryclark@google.com4fdbb222013-07-23 15:27:41 +0000168SkDPoint SkDQuad::ptAtT(double t) const {
169 if (0 == t) {
170 return fPts[0];
171 }
172 if (1 == t) {
173 return fPts[2];
174 }
caryclark@google.com07393ca2013-04-08 11:47:37 +0000175 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
184/*
185Given a quadratic q, t1, and t2, find a small quadratic segment.
186
187The new quadratic is defined by A, B, and C, where
188 A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
189 C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
190
191To find B, compute the point halfway between t1 and t2:
192
193q(at (t1 + t2)/2) == D
194
195Next, compute where D must be if we know the value of B:
196
197_12 = A/2 + B/2
19812_ = B/2 + C/2
199123 = A/4 + B/2 + C/4
200 = D
201
202Group the known values on one side:
203
204B = D*2 - A/2 - C/2
205*/
206
207static double interp_quad_coords(const double* src, double t) {
208 double ab = SkDInterp(src[0], src[2], t);
209 double bc = SkDInterp(src[2], src[4], t);
210 double abc = SkDInterp(ab, bc, t);
211 return abc;
212}
213
214bool SkDQuad::monotonicInY() const {
215 return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
216}
217
218SkDQuad SkDQuad::subDivide(double t1, double t2) const {
219 SkDQuad dst;
220 double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
221 double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
222 double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
223 double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
224 double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
225 double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
226 /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2;
227 /* by = */ dst[1].fY = 2*dy - (ay + cy)/2;
228 return dst;
229}
230
skia.committer@gmail.com8f6ef402013-06-05 07:01:06 +0000231void SkDQuad::align(int endIndex, SkDPoint* dstPt) const {
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000232 if (fPts[endIndex].fX == fPts[1].fX) {
233 dstPt->fX = fPts[endIndex].fX;
234 }
235 if (fPts[endIndex].fY == fPts[1].fY) {
236 dstPt->fY = fPts[endIndex].fY;
237 }
238}
239
caryclark@google.com07393ca2013-04-08 11:47:37 +0000240SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000241 SkASSERT(t1 != t2);
caryclark@google.com07393ca2013-04-08 11:47:37 +0000242 SkDPoint b;
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000243#if 0
244 // this approach assumes that the control point computed directly is accurate enough
caryclark@google.com07393ca2013-04-08 11:47:37 +0000245 double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
246 double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
247 b.fX = 2 * dx - (a.fX + c.fX) / 2;
248 b.fY = 2 * dy - (a.fY + c.fY) / 2;
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000249#else
250 SkDQuad sub = subDivide(t1, t2);
251 SkDLine b0 = {{a, sub[1] + (a - sub[0])}};
252 SkDLine b1 = {{c, sub[1] + (c - sub[2])}};
253 SkIntersections i;
254 i.intersectRay(b0, b1);
commit-bot@chromium.org4431e772014-04-14 17:08:59 +0000255 if (i.used() == 1 && i[0][0] >= 0 && i[1][0] >= 0) {
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000256 b = i.pt(0);
257 } else {
commit-bot@chromium.org4431e772014-04-14 17:08:59 +0000258 SkASSERT(i.used() <= 2);
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000259 b = SkDPoint::Mid(b0[1], b1[1]);
260 }
261#endif
262 if (t1 == 0 || t2 == 0) {
263 align(0, &b);
264 }
265 if (t1 == 1 || t2 == 1) {
266 align(2, &b);
267 }
commit-bot@chromium.org4431e772014-04-14 17:08:59 +0000268 if (AlmostBequalUlps(b.fX, a.fX)) {
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000269 b.fX = a.fX;
commit-bot@chromium.org4431e772014-04-14 17:08:59 +0000270 } else if (AlmostBequalUlps(b.fX, c.fX)) {
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000271 b.fX = c.fX;
272 }
commit-bot@chromium.org4431e772014-04-14 17:08:59 +0000273 if (AlmostBequalUlps(b.fY, a.fY)) {
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000274 b.fY = a.fY;
commit-bot@chromium.org4431e772014-04-14 17:08:59 +0000275 } else if (AlmostBequalUlps(b.fY, c.fY)) {
caryclark@google.comcffbcc32013-06-04 17:59:42 +0000276 b.fY = c.fY;
277 }
caryclark@google.com07393ca2013-04-08 11:47:37 +0000278 return b;
279}
280
281/* classic one t subdivision */
282static void interp_quad_coords(const double* src, double* dst, double t) {
283 double ab = SkDInterp(src[0], src[2], t);
284 double bc = SkDInterp(src[2], src[4], t);
285 dst[0] = src[0];
286 dst[2] = ab;
287 dst[4] = SkDInterp(ab, bc, t);
288 dst[6] = bc;
289 dst[8] = src[4];
290}
291
292SkDQuadPair SkDQuad::chopAt(double t) const
293{
294 SkDQuadPair dst;
295 interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
296 interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
297 return dst;
298}
299
300static int valid_unit_divide(double numer, double denom, double* ratio)
301{
302 if (numer < 0) {
303 numer = -numer;
304 denom = -denom;
305 }
306 if (denom == 0 || numer == 0 || numer >= denom) {
307 return 0;
308 }
309 double r = numer / denom;
310 if (r == 0) { // catch underflow if numer <<<< denom
311 return 0;
312 }
313 *ratio = r;
314 return 1;
315}
316
317/** Quad'(t) = At + B, where
318 A = 2(a - 2b + c)
319 B = 2(b - a)
320 Solve for t, only if it fits between 0 < t < 1
321*/
322int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) {
323 /* At + B == 0
324 t = -B / A
325 */
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 */
335void 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}