SkPathOpsQuad.cpp revision 8f6ef4010f6835c5ce9ede180e50a6a58512a81e 
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 "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]);
30    double distMin = SkTMin(d0, d2);
31    int bestIndex = -1;
32    for (int index = 0; index < roots; ++index) {
33        SkDPoint onQuad = xyAtT(ts[index]);
34        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;
60            SkDPoint test = xyAtT(extremeT);
61            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;
125    if (!AlmostEqualUlps(p2, q) && p2 < q) {
126        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;
134    return 1 + !AlmostEqualUlps(s[0], s[1]);
135}
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
167SkDPoint SkDQuad::xyAtT(double t) const {
168    double one_t = 1 - t;
169    double a = one_t * one_t;
170    double b = 2 * one_t * t;
171    double c = t * t;
172    SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
173            a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
174    return result;
175}
176
177/*
178Given a quadratic q, t1, and t2, find a small quadratic segment.
179
180The new quadratic is defined by A, B, and C, where
181 A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
182 C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
183
184To find B, compute the point halfway between t1 and t2:
185
186q(at (t1 + t2)/2) == D
187
188Next, compute where D must be if we know the value of B:
189
190_12 = A/2 + B/2
19112_ = B/2 + C/2
192123 = A/4 + B/2 + C/4
193    = D
194
195Group the known values on one side:
196
197B   = D*2 - A/2 - C/2
198*/
199
200static double interp_quad_coords(const double* src, double t) {
201    double ab = SkDInterp(src[0], src[2], t);
202    double bc = SkDInterp(src[2], src[4], t);
203    double abc = SkDInterp(ab, bc, t);
204    return abc;
205}
206
207bool SkDQuad::monotonicInY() const {
208    return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
209}
210
211SkDQuad SkDQuad::subDivide(double t1, double t2) const {
212    SkDQuad dst;
213    double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
214    double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
215    double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
216    double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
217    double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
218    double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
219    /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2;
220    /* by = */ dst[1].fY = 2*dy - (ay + cy)/2;
221    return dst;
222}
223
224void SkDQuad::align(int endIndex, SkDPoint* dstPt) const {
225    if (fPts[endIndex].fX == fPts[1].fX) {
226        dstPt->fX = fPts[endIndex].fX;
227    }
228    if (fPts[endIndex].fY == fPts[1].fY) {
229        dstPt->fY = fPts[endIndex].fY;
230    }
231}
232
233SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
234    SkASSERT(t1 != t2);
235    SkDPoint b;
236#if 0
237    // this approach assumes that the control point computed directly is accurate enough
238    double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
239    double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
240    b.fX = 2 * dx - (a.fX + c.fX) / 2;
241    b.fY = 2 * dy - (a.fY + c.fY) / 2;
242#else
243    SkDQuad sub = subDivide(t1, t2);
244    SkDLine b0 = {{a, sub[1] + (a - sub[0])}};
245    SkDLine b1 = {{c, sub[1] + (c - sub[2])}};
246    SkIntersections i;
247    i.intersectRay(b0, b1);
248    if (i.used() == 1) {
249        b = i.pt(0);
250    } else {
251        SkASSERT(i.used() == 2 || i.used() == 0);
252        b = SkDPoint::Mid(b0[1], b1[1]);
253    }
254#endif
255    if (t1 == 0 || t2 == 0) {
256        align(0, &b);
257    }
258    if (t1 == 1 || t2 == 1) {
259        align(2, &b);
260    }
261    if (precisely_subdivide_equal(b.fX, a.fX)) {
262        b.fX = a.fX;
263    } else if (precisely_subdivide_equal(b.fX, c.fX)) {
264        b.fX = c.fX;
265    }
266    if (precisely_subdivide_equal(b.fY, a.fY)) {
267        b.fY = a.fY;
268    } else if (precisely_subdivide_equal(b.fY, c.fY)) {
269        b.fY = c.fY;
270    }
271    return b;
272}
273
274/* classic one t subdivision */
275static void interp_quad_coords(const double* src, double* dst, double t) {
276    double ab = SkDInterp(src[0], src[2], t);
277    double bc = SkDInterp(src[2], src[4], t);
278    dst[0] = src[0];
279    dst[2] = ab;
280    dst[4] = SkDInterp(ab, bc, t);
281    dst[6] = bc;
282    dst[8] = src[4];
283}
284
285SkDQuadPair SkDQuad::chopAt(double t) const
286{
287    SkDQuadPair dst;
288    interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
289    interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
290    return dst;
291}
292
293static int valid_unit_divide(double numer, double denom, double* ratio)
294{
295    if (numer < 0) {
296        numer = -numer;
297        denom = -denom;
298    }
299    if (denom == 0 || numer == 0 || numer >= denom) {
300        return 0;
301    }
302    double r = numer / denom;
303    if (r == 0) {  // catch underflow if numer <<<< denom
304        return 0;
305    }
306    *ratio = r;
307    return 1;
308}
309
310/** Quad'(t) = At + B, where
311    A = 2(a - 2b + c)
312    B = 2(b - a)
313    Solve for t, only if it fits between 0 < t < 1
314*/
315int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) {
316    /*  At + B == 0
317        t = -B / A
318    */
319    return valid_unit_divide(a - b, a - b - b + c, tValue);
320}
321
322/* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
323 *
324 * a = A - 2*B +   C
325 * b =     2*B - 2*C
326 * c =             C
327 */
328void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) {
329    *a = quad[0];      // a = A
330    *b = 2 * quad[2];  // b =     2*B
331    *c = quad[4];      // c =             C
332    *b -= *c;          // b =     2*B -   C
333    *a -= *b;          // a = A - 2*B +   C
334    *b -= *c;          // b =     2*B - 2*C
335}
336