1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/gfx/skia/trunk/src/pathops/SkPathOpsQuad.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,355 @@
1.4 +/*
1.5 + * Copyright 2012 Google Inc.
1.6 + *
1.7 + * Use of this source code is governed by a BSD-style license that can be
1.8 + * found in the LICENSE file.
1.9 + */
1.10 +#include "SkIntersections.h"
1.11 +#include "SkLineParameters.h"
1.12 +#include "SkPathOpsCubic.h"
1.13 +#include "SkPathOpsQuad.h"
1.14 +#include "SkPathOpsTriangle.h"
1.15 +
1.16 +// from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
1.17 +// (currently only used by testing)
1.18 +double SkDQuad::nearestT(const SkDPoint& pt) const {
1.19 + SkDVector pos = fPts[0] - pt;
1.20 + // search points P of bezier curve with PM.(dP / dt) = 0
1.21 + // a calculus leads to a 3d degree equation :
1.22 + SkDVector A = fPts[1] - fPts[0];
1.23 + SkDVector B = fPts[2] - fPts[1];
1.24 + B -= A;
1.25 + double a = B.dot(B);
1.26 + double b = 3 * A.dot(B);
1.27 + double c = 2 * A.dot(A) + pos.dot(B);
1.28 + double d = pos.dot(A);
1.29 + double ts[3];
1.30 + int roots = SkDCubic::RootsValidT(a, b, c, d, ts);
1.31 + double d0 = pt.distanceSquared(fPts[0]);
1.32 + double d2 = pt.distanceSquared(fPts[2]);
1.33 + double distMin = SkTMin(d0, d2);
1.34 + int bestIndex = -1;
1.35 + for (int index = 0; index < roots; ++index) {
1.36 + SkDPoint onQuad = ptAtT(ts[index]);
1.37 + double dist = pt.distanceSquared(onQuad);
1.38 + if (distMin > dist) {
1.39 + distMin = dist;
1.40 + bestIndex = index;
1.41 + }
1.42 + }
1.43 + if (bestIndex >= 0) {
1.44 + return ts[bestIndex];
1.45 + }
1.46 + return d0 < d2 ? 0 : 1;
1.47 +}
1.48 +
1.49 +bool SkDQuad::pointInHull(const SkDPoint& pt) const {
1.50 + return ((const SkDTriangle&) fPts).contains(pt);
1.51 +}
1.52 +
1.53 +SkDPoint SkDQuad::top(double startT, double endT) const {
1.54 + SkDQuad sub = subDivide(startT, endT);
1.55 + SkDPoint topPt = sub[0];
1.56 + if (topPt.fY > sub[2].fY || (topPt.fY == sub[2].fY && topPt.fX > sub[2].fX)) {
1.57 + topPt = sub[2];
1.58 + }
1.59 + if (!between(sub[0].fY, sub[1].fY, sub[2].fY)) {
1.60 + double extremeT;
1.61 + if (FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, &extremeT)) {
1.62 + extremeT = startT + (endT - startT) * extremeT;
1.63 + SkDPoint test = ptAtT(extremeT);
1.64 + if (topPt.fY > test.fY || (topPt.fY == test.fY && topPt.fX > test.fX)) {
1.65 + topPt = test;
1.66 + }
1.67 + }
1.68 + }
1.69 + return topPt;
1.70 +}
1.71 +
1.72 +int SkDQuad::AddValidTs(double s[], int realRoots, double* t) {
1.73 + int foundRoots = 0;
1.74 + for (int index = 0; index < realRoots; ++index) {
1.75 + double tValue = s[index];
1.76 + if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
1.77 + if (approximately_less_than_zero(tValue)) {
1.78 + tValue = 0;
1.79 + } else if (approximately_greater_than_one(tValue)) {
1.80 + tValue = 1;
1.81 + }
1.82 + for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
1.83 + if (approximately_equal(t[idx2], tValue)) {
1.84 + goto nextRoot;
1.85 + }
1.86 + }
1.87 + t[foundRoots++] = tValue;
1.88 + }
1.89 +nextRoot:
1.90 + {}
1.91 + }
1.92 + return foundRoots;
1.93 +}
1.94 +
1.95 +// note: caller expects multiple results to be sorted smaller first
1.96 +// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
1.97 +// analysis of the quadratic equation, suggesting why the following looks at
1.98 +// the sign of B -- and further suggesting that the greatest loss of precision
1.99 +// is in b squared less two a c
1.100 +int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) {
1.101 + double s[2];
1.102 + int realRoots = RootsReal(A, B, C, s);
1.103 + int foundRoots = AddValidTs(s, realRoots, t);
1.104 + return foundRoots;
1.105 +}
1.106 +
1.107 +/*
1.108 +Numeric Solutions (5.6) suggests to solve the quadratic by computing
1.109 + Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
1.110 +and using the roots
1.111 + t1 = Q / A
1.112 + t2 = C / Q
1.113 +*/
1.114 +// this does not discard real roots <= 0 or >= 1
1.115 +int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) {
1.116 + const double p = B / (2 * A);
1.117 + const double q = C / A;
1.118 + if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) {
1.119 + if (approximately_zero(B)) {
1.120 + s[0] = 0;
1.121 + return C == 0;
1.122 + }
1.123 + s[0] = -C / B;
1.124 + return 1;
1.125 + }
1.126 + /* normal form: x^2 + px + q = 0 */
1.127 + const double p2 = p * p;
1.128 + if (!AlmostDequalUlps(p2, q) && p2 < q) {
1.129 + return 0;
1.130 + }
1.131 + double sqrt_D = 0;
1.132 + if (p2 > q) {
1.133 + sqrt_D = sqrt(p2 - q);
1.134 + }
1.135 + s[0] = sqrt_D - p;
1.136 + s[1] = -sqrt_D - p;
1.137 + return 1 + !AlmostDequalUlps(s[0], s[1]);
1.138 +}
1.139 +
1.140 +bool SkDQuad::isLinear(int startIndex, int endIndex) const {
1.141 + SkLineParameters lineParameters;
1.142 + lineParameters.quadEndPoints(*this, startIndex, endIndex);
1.143 + // FIXME: maybe it's possible to avoid this and compare non-normalized
1.144 + lineParameters.normalize();
1.145 + double distance = lineParameters.controlPtDistance(*this);
1.146 + return approximately_zero(distance);
1.147 +}
1.148 +
1.149 +SkDCubic SkDQuad::toCubic() const {
1.150 + SkDCubic cubic;
1.151 + cubic[0] = fPts[0];
1.152 + cubic[2] = fPts[1];
1.153 + cubic[3] = fPts[2];
1.154 + cubic[1].fX = (cubic[0].fX + cubic[2].fX * 2) / 3;
1.155 + cubic[1].fY = (cubic[0].fY + cubic[2].fY * 2) / 3;
1.156 + cubic[2].fX = (cubic[3].fX + cubic[2].fX * 2) / 3;
1.157 + cubic[2].fY = (cubic[3].fY + cubic[2].fY * 2) / 3;
1.158 + return cubic;
1.159 +}
1.160 +
1.161 +SkDVector SkDQuad::dxdyAtT(double t) const {
1.162 + double a = t - 1;
1.163 + double b = 1 - 2 * t;
1.164 + double c = t;
1.165 + SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
1.166 + a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
1.167 + return result;
1.168 +}
1.169 +
1.170 +// OPTIMIZE: assert if caller passes in t == 0 / t == 1 ?
1.171 +SkDPoint SkDQuad::ptAtT(double t) const {
1.172 + if (0 == t) {
1.173 + return fPts[0];
1.174 + }
1.175 + if (1 == t) {
1.176 + return fPts[2];
1.177 + }
1.178 + double one_t = 1 - t;
1.179 + double a = one_t * one_t;
1.180 + double b = 2 * one_t * t;
1.181 + double c = t * t;
1.182 + SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
1.183 + a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
1.184 + return result;
1.185 +}
1.186 +
1.187 +/*
1.188 +Given a quadratic q, t1, and t2, find a small quadratic segment.
1.189 +
1.190 +The new quadratic is defined by A, B, and C, where
1.191 + A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
1.192 + C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
1.193 +
1.194 +To find B, compute the point halfway between t1 and t2:
1.195 +
1.196 +q(at (t1 + t2)/2) == D
1.197 +
1.198 +Next, compute where D must be if we know the value of B:
1.199 +
1.200 +_12 = A/2 + B/2
1.201 +12_ = B/2 + C/2
1.202 +123 = A/4 + B/2 + C/4
1.203 + = D
1.204 +
1.205 +Group the known values on one side:
1.206 +
1.207 +B = D*2 - A/2 - C/2
1.208 +*/
1.209 +
1.210 +static double interp_quad_coords(const double* src, double t) {
1.211 + double ab = SkDInterp(src[0], src[2], t);
1.212 + double bc = SkDInterp(src[2], src[4], t);
1.213 + double abc = SkDInterp(ab, bc, t);
1.214 + return abc;
1.215 +}
1.216 +
1.217 +bool SkDQuad::monotonicInY() const {
1.218 + return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
1.219 +}
1.220 +
1.221 +SkDQuad SkDQuad::subDivide(double t1, double t2) const {
1.222 + SkDQuad dst;
1.223 + double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
1.224 + double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
1.225 + double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
1.226 + double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
1.227 + double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
1.228 + double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
1.229 + /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2;
1.230 + /* by = */ dst[1].fY = 2*dy - (ay + cy)/2;
1.231 + return dst;
1.232 +}
1.233 +
1.234 +void SkDQuad::align(int endIndex, SkDPoint* dstPt) const {
1.235 + if (fPts[endIndex].fX == fPts[1].fX) {
1.236 + dstPt->fX = fPts[endIndex].fX;
1.237 + }
1.238 + if (fPts[endIndex].fY == fPts[1].fY) {
1.239 + dstPt->fY = fPts[endIndex].fY;
1.240 + }
1.241 +}
1.242 +
1.243 +SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
1.244 + SkASSERT(t1 != t2);
1.245 + SkDPoint b;
1.246 +#if 0
1.247 + // this approach assumes that the control point computed directly is accurate enough
1.248 + double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
1.249 + double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
1.250 + b.fX = 2 * dx - (a.fX + c.fX) / 2;
1.251 + b.fY = 2 * dy - (a.fY + c.fY) / 2;
1.252 +#else
1.253 + SkDQuad sub = subDivide(t1, t2);
1.254 + SkDLine b0 = {{a, sub[1] + (a - sub[0])}};
1.255 + SkDLine b1 = {{c, sub[1] + (c - sub[2])}};
1.256 + SkIntersections i;
1.257 + i.intersectRay(b0, b1);
1.258 + if (i.used() == 1) {
1.259 + b = i.pt(0);
1.260 + } else {
1.261 + SkASSERT(i.used() == 2 || i.used() == 0);
1.262 + b = SkDPoint::Mid(b0[1], b1[1]);
1.263 + }
1.264 +#endif
1.265 + if (t1 == 0 || t2 == 0) {
1.266 + align(0, &b);
1.267 + }
1.268 + if (t1 == 1 || t2 == 1) {
1.269 + align(2, &b);
1.270 + }
1.271 + if (precisely_subdivide_equal(b.fX, a.fX)) {
1.272 + b.fX = a.fX;
1.273 + } else if (precisely_subdivide_equal(b.fX, c.fX)) {
1.274 + b.fX = c.fX;
1.275 + }
1.276 + if (precisely_subdivide_equal(b.fY, a.fY)) {
1.277 + b.fY = a.fY;
1.278 + } else if (precisely_subdivide_equal(b.fY, c.fY)) {
1.279 + b.fY = c.fY;
1.280 + }
1.281 + return b;
1.282 +}
1.283 +
1.284 +/* classic one t subdivision */
1.285 +static void interp_quad_coords(const double* src, double* dst, double t) {
1.286 + double ab = SkDInterp(src[0], src[2], t);
1.287 + double bc = SkDInterp(src[2], src[4], t);
1.288 + dst[0] = src[0];
1.289 + dst[2] = ab;
1.290 + dst[4] = SkDInterp(ab, bc, t);
1.291 + dst[6] = bc;
1.292 + dst[8] = src[4];
1.293 +}
1.294 +
1.295 +SkDQuadPair SkDQuad::chopAt(double t) const
1.296 +{
1.297 + SkDQuadPair dst;
1.298 + interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
1.299 + interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
1.300 + return dst;
1.301 +}
1.302 +
1.303 +static int valid_unit_divide(double numer, double denom, double* ratio)
1.304 +{
1.305 + if (numer < 0) {
1.306 + numer = -numer;
1.307 + denom = -denom;
1.308 + }
1.309 + if (denom == 0 || numer == 0 || numer >= denom) {
1.310 + return 0;
1.311 + }
1.312 + double r = numer / denom;
1.313 + if (r == 0) { // catch underflow if numer <<<< denom
1.314 + return 0;
1.315 + }
1.316 + *ratio = r;
1.317 + return 1;
1.318 +}
1.319 +
1.320 +/** Quad'(t) = At + B, where
1.321 + A = 2(a - 2b + c)
1.322 + B = 2(b - a)
1.323 + Solve for t, only if it fits between 0 < t < 1
1.324 +*/
1.325 +int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) {
1.326 + /* At + B == 0
1.327 + t = -B / A
1.328 + */
1.329 + return valid_unit_divide(a - b, a - b - b + c, tValue);
1.330 +}
1.331 +
1.332 +/* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
1.333 + *
1.334 + * a = A - 2*B + C
1.335 + * b = 2*B - 2*C
1.336 + * c = C
1.337 + */
1.338 +void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) {
1.339 + *a = quad[0]; // a = A
1.340 + *b = 2 * quad[2]; // b = 2*B
1.341 + *c = quad[4]; // c = C
1.342 + *b -= *c; // b = 2*B - C
1.343 + *a -= *b; // a = A - 2*B + C
1.344 + *b -= *c; // b = 2*B - 2*C
1.345 +}
1.346 +
1.347 +#ifdef SK_DEBUG
1.348 +void SkDQuad::dump() {
1.349 + SkDebugf("{{");
1.350 + int index = 0;
1.351 + do {
1.352 + fPts[index].dump();
1.353 + SkDebugf(", ");
1.354 + } while (++index < 2);
1.355 + fPts[index].dump();
1.356 + SkDebugf("}}\n");
1.357 +}
1.358 +#endif
