blob: 685f49e16044e235bd533c699f23fbecb0d86423 [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 */
7#include "SkLineParameters.h"
8#include "SkPathOpsCubic.h"
9#include "SkPathOpsQuad.h"
10#include "SkPathOpsTriangle.h"
11
12// from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
13// (currently only used by testing)
14double SkDQuad::nearestT(const SkDPoint& pt) const {
15 SkDVector pos = fPts[0] - pt;
16 // search points P of bezier curve with PM.(dP / dt) = 0
17 // a calculus leads to a 3d degree equation :
18 SkDVector A = fPts[1] - fPts[0];
19 SkDVector B = fPts[2] - fPts[1];
20 B -= A;
21 double a = B.dot(B);
22 double b = 3 * A.dot(B);
23 double c = 2 * A.dot(A) + pos.dot(B);
24 double d = pos.dot(A);
25 double ts[3];
26 int roots = SkDCubic::RootsValidT(a, b, c, d, ts);
27 double d0 = pt.distanceSquared(fPts[0]);
28 double d2 = pt.distanceSquared(fPts[2]);
caryclark@google.com3b97af52013-04-23 11:56:44 +000029 double distMin = SkTMin(d0, d2);
caryclark@google.com07393ca2013-04-08 11:47:37 +000030 int bestIndex = -1;
31 for (int index = 0; index < roots; ++index) {
32 SkDPoint onQuad = xyAtT(ts[index]);
33 double dist = pt.distanceSquared(onQuad);
34 if (distMin > dist) {
35 distMin = dist;
36 bestIndex = index;
37 }
38 }
39 if (bestIndex >= 0) {
40 return ts[bestIndex];
41 }
42 return d0 < d2 ? 0 : 1;
43}
44
45bool SkDQuad::pointInHull(const SkDPoint& pt) const {
46 return ((const SkDTriangle&) fPts).contains(pt);
47}
48
49SkDPoint SkDQuad::top(double startT, double endT) const {
50 SkDQuad sub = subDivide(startT, endT);
51 SkDPoint topPt = sub[0];
52 if (topPt.fY > sub[2].fY || (topPt.fY == sub[2].fY && topPt.fX > sub[2].fX)) {
53 topPt = sub[2];
54 }
55 if (!between(sub[0].fY, sub[1].fY, sub[2].fY)) {
56 double extremeT;
57 if (FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, &extremeT)) {
58 extremeT = startT + (endT - startT) * extremeT;
59 SkDPoint test = xyAtT(extremeT);
60 if (topPt.fY > test.fY || (topPt.fY == test.fY && topPt.fX > test.fX)) {
61 topPt = test;
62 }
63 }
64 }
65 return topPt;
66}
67
68int SkDQuad::AddValidTs(double s[], int realRoots, double* t) {
69 int foundRoots = 0;
70 for (int index = 0; index < realRoots; ++index) {
71 double tValue = s[index];
72 if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
73 if (approximately_less_than_zero(tValue)) {
74 tValue = 0;
75 } else if (approximately_greater_than_one(tValue)) {
76 tValue = 1;
77 }
78 for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
79 if (approximately_equal(t[idx2], tValue)) {
80 goto nextRoot;
81 }
82 }
83 t[foundRoots++] = tValue;
84 }
85nextRoot:
86 {}
87 }
88 return foundRoots;
89}
90
91// note: caller expects multiple results to be sorted smaller first
92// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
93// analysis of the quadratic equation, suggesting why the following looks at
94// the sign of B -- and further suggesting that the greatest loss of precision
95// is in b squared less two a c
96int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) {
97 double s[2];
98 int realRoots = RootsReal(A, B, C, s);
99 int foundRoots = AddValidTs(s, realRoots, t);
100 return foundRoots;
101}
102
103/*
104Numeric Solutions (5.6) suggests to solve the quadratic by computing
105 Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
106and using the roots
107 t1 = Q / A
108 t2 = C / Q
109*/
110// this does not discard real roots <= 0 or >= 1
111int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) {
112 const double p = B / (2 * A);
113 const double q = C / A;
114 if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) {
115 if (approximately_zero(B)) {
116 s[0] = 0;
117 return C == 0;
118 }
119 s[0] = -C / B;
120 return 1;
121 }
122 /* normal form: x^2 + px + q = 0 */
123 const double p2 = p * p;
124 if (!AlmostEqualUlps(p2, q) && p2 < q) {
125 return 0;
126 }
127 double sqrt_D = 0;
128 if (p2 > q) {
129 sqrt_D = sqrt(p2 - q);
130 }
131 s[0] = sqrt_D - p;
132 s[1] = -sqrt_D - p;
133 return 1 + !AlmostEqualUlps(s[0], s[1]);
134}
135
136bool SkDQuad::isLinear(int startIndex, int endIndex) const {
137 SkLineParameters lineParameters;
138 lineParameters.quadEndPoints(*this, startIndex, endIndex);
139 // FIXME: maybe it's possible to avoid this and compare non-normalized
140 lineParameters.normalize();
141 double distance = lineParameters.controlPtDistance(*this);
142 return approximately_zero(distance);
143}
144
145SkDCubic SkDQuad::toCubic() const {
146 SkDCubic cubic;
147 cubic[0] = fPts[0];
148 cubic[2] = fPts[1];
149 cubic[3] = fPts[2];
150 cubic[1].fX = (cubic[0].fX + cubic[2].fX * 2) / 3;
151 cubic[1].fY = (cubic[0].fY + cubic[2].fY * 2) / 3;
152 cubic[2].fX = (cubic[3].fX + cubic[2].fX * 2) / 3;
153 cubic[2].fY = (cubic[3].fY + cubic[2].fY * 2) / 3;
154 return cubic;
155}
156
157SkDVector SkDQuad::dxdyAtT(double t) const {
158 double a = t - 1;
159 double b = 1 - 2 * t;
160 double c = t;
161 SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
162 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
163 return result;
164}
165
166SkDPoint SkDQuad::xyAtT(double t) const {
167 double one_t = 1 - t;
168 double a = one_t * one_t;
169 double b = 2 * one_t * t;
170 double c = t * t;
171 SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
172 a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
173 return result;
174}
175
176/*
177Given a quadratic q, t1, and t2, find a small quadratic segment.
178
179The new quadratic is defined by A, B, and C, where
180 A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
181 C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
182
183To find B, compute the point halfway between t1 and t2:
184
185q(at (t1 + t2)/2) == D
186
187Next, compute where D must be if we know the value of B:
188
189_12 = A/2 + B/2
19012_ = B/2 + C/2
191123 = A/4 + B/2 + C/4
192 = D
193
194Group the known values on one side:
195
196B = D*2 - A/2 - C/2
197*/
198
199static double interp_quad_coords(const double* src, double t) {
200 double ab = SkDInterp(src[0], src[2], t);
201 double bc = SkDInterp(src[2], src[4], t);
202 double abc = SkDInterp(ab, bc, t);
203 return abc;
204}
205
206bool SkDQuad::monotonicInY() const {
207 return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
208}
209
210SkDQuad SkDQuad::subDivide(double t1, double t2) const {
211 SkDQuad dst;
212 double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
213 double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
214 double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
215 double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
216 double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
217 double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
218 /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2;
219 /* by = */ dst[1].fY = 2*dy - (ay + cy)/2;
220 return dst;
221}
222
223SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
224 SkDPoint b;
225 double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
226 double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
227 b.fX = 2 * dx - (a.fX + c.fX) / 2;
228 b.fY = 2 * dy - (a.fY + c.fY) / 2;
229 return b;
230}
231
232/* classic one t subdivision */
233static void interp_quad_coords(const double* src, double* dst, double t) {
234 double ab = SkDInterp(src[0], src[2], t);
235 double bc = SkDInterp(src[2], src[4], t);
236 dst[0] = src[0];
237 dst[2] = ab;
238 dst[4] = SkDInterp(ab, bc, t);
239 dst[6] = bc;
240 dst[8] = src[4];
241}
242
243SkDQuadPair SkDQuad::chopAt(double t) const
244{
245 SkDQuadPair dst;
246 interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
247 interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
248 return dst;
249}
250
251static int valid_unit_divide(double numer, double denom, double* ratio)
252{
253 if (numer < 0) {
254 numer = -numer;
255 denom = -denom;
256 }
257 if (denom == 0 || numer == 0 || numer >= denom) {
258 return 0;
259 }
260 double r = numer / denom;
261 if (r == 0) { // catch underflow if numer <<<< denom
262 return 0;
263 }
264 *ratio = r;
265 return 1;
266}
267
268/** Quad'(t) = At + B, where
269 A = 2(a - 2b + c)
270 B = 2(b - a)
271 Solve for t, only if it fits between 0 < t < 1
272*/
273int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) {
274 /* At + B == 0
275 t = -B / A
276 */
277 return valid_unit_divide(a - b, a - b - b + c, tValue);
278}
279
280/* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
281 *
282 * a = A - 2*B + C
283 * b = 2*B - 2*C
284 * c = C
285 */
286void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) {
287 *a = quad[0]; // a = A
288 *b = 2 * quad[2]; // b = 2*B
289 *c = quad[4]; // c = C
290 *b -= *c; // b = 2*B - C
291 *a -= *b; // a = A - 2*B + C
292 *b -= *c; // b = 2*B - 2*C
293}