diff -r 000000000000 -r 6474c204b198 gfx/skia/trunk/src/pathops/SkPathOpsQuad.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gfx/skia/trunk/src/pathops/SkPathOpsQuad.cpp Wed Dec 31 06:09:35 2014 +0100 @@ -0,0 +1,355 @@ +/* + * Copyright 2012 Google Inc. + * + * Use of this source code is governed by a BSD-style license that can be + * found in the LICENSE file. + */ +#include "SkIntersections.h" +#include "SkLineParameters.h" +#include "SkPathOpsCubic.h" +#include "SkPathOpsQuad.h" +#include "SkPathOpsTriangle.h" + +// from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html +// (currently only used by testing) +double SkDQuad::nearestT(const SkDPoint& pt) const { + SkDVector pos = fPts[0] - pt; + // search points P of bezier curve with PM.(dP / dt) = 0 + // a calculus leads to a 3d degree equation : + SkDVector A = fPts[1] - fPts[0]; + SkDVector B = fPts[2] - fPts[1]; + B -= A; + double a = B.dot(B); + double b = 3 * A.dot(B); + double c = 2 * A.dot(A) + pos.dot(B); + double d = pos.dot(A); + double ts[3]; + int roots = SkDCubic::RootsValidT(a, b, c, d, ts); + double d0 = pt.distanceSquared(fPts[0]); + double d2 = pt.distanceSquared(fPts[2]); + double distMin = SkTMin(d0, d2); + int bestIndex = -1; + for (int index = 0; index < roots; ++index) { + SkDPoint onQuad = ptAtT(ts[index]); + double dist = pt.distanceSquared(onQuad); + if (distMin > dist) { + distMin = dist; + bestIndex = index; + } + } + if (bestIndex >= 0) { + return ts[bestIndex]; + } + return d0 < d2 ? 0 : 1; +} + +bool SkDQuad::pointInHull(const SkDPoint& pt) const { + return ((const SkDTriangle&) fPts).contains(pt); +} + +SkDPoint SkDQuad::top(double startT, double endT) const { + SkDQuad sub = subDivide(startT, endT); + SkDPoint topPt = sub[0]; + if (topPt.fY > sub[2].fY || (topPt.fY == sub[2].fY && topPt.fX > sub[2].fX)) { + topPt = sub[2]; + } + if (!between(sub[0].fY, sub[1].fY, sub[2].fY)) { + double extremeT; + if (FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, &extremeT)) { + extremeT = startT + (endT - startT) * extremeT; + SkDPoint test = ptAtT(extremeT); + if (topPt.fY > test.fY || (topPt.fY == test.fY && topPt.fX > test.fX)) { + topPt = test; + } + } + } + return topPt; +} + +int SkDQuad::AddValidTs(double s[], int realRoots, double* t) { + int foundRoots = 0; + for (int index = 0; index < realRoots; ++index) { + double tValue = s[index]; + if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) { + if (approximately_less_than_zero(tValue)) { + tValue = 0; + } else if (approximately_greater_than_one(tValue)) { + tValue = 1; + } + for (int idx2 = 0; idx2 < foundRoots; ++idx2) { + if (approximately_equal(t[idx2], tValue)) { + goto nextRoot; + } + } + t[foundRoots++] = tValue; + } +nextRoot: + {} + } + return foundRoots; +} + +// note: caller expects multiple results to be sorted smaller first +// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting +// analysis of the quadratic equation, suggesting why the following looks at +// the sign of B -- and further suggesting that the greatest loss of precision +// is in b squared less two a c +int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) { + double s[2]; + int realRoots = RootsReal(A, B, C, s); + int foundRoots = AddValidTs(s, realRoots, t); + return foundRoots; +} + +/* +Numeric Solutions (5.6) suggests to solve the quadratic by computing + Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C)) +and using the roots + t1 = Q / A + t2 = C / Q +*/ +// this does not discard real roots <= 0 or >= 1 +int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) { + const double p = B / (2 * A); + const double q = C / A; + if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) { + if (approximately_zero(B)) { + s[0] = 0; + return C == 0; + } + s[0] = -C / B; + return 1; + } + /* normal form: x^2 + px + q = 0 */ + const double p2 = p * p; + if (!AlmostDequalUlps(p2, q) && p2 < q) { + return 0; + } + double sqrt_D = 0; + if (p2 > q) { + sqrt_D = sqrt(p2 - q); + } + s[0] = sqrt_D - p; + s[1] = -sqrt_D - p; + return 1 + !AlmostDequalUlps(s[0], s[1]); +} + +bool SkDQuad::isLinear(int startIndex, int endIndex) const { + SkLineParameters lineParameters; + lineParameters.quadEndPoints(*this, startIndex, endIndex); + // FIXME: maybe it's possible to avoid this and compare non-normalized + lineParameters.normalize(); + double distance = lineParameters.controlPtDistance(*this); + return approximately_zero(distance); +} + +SkDCubic SkDQuad::toCubic() const { + SkDCubic cubic; + cubic[0] = fPts[0]; + cubic[2] = fPts[1]; + cubic[3] = fPts[2]; + cubic[1].fX = (cubic[0].fX + cubic[2].fX * 2) / 3; + cubic[1].fY = (cubic[0].fY + cubic[2].fY * 2) / 3; + cubic[2].fX = (cubic[3].fX + cubic[2].fX * 2) / 3; + cubic[2].fY = (cubic[3].fY + cubic[2].fY * 2) / 3; + return cubic; +} + +SkDVector SkDQuad::dxdyAtT(double t) const { + double a = t - 1; + double b = 1 - 2 * t; + double c = t; + SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX, + a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY }; + return result; +} + +// OPTIMIZE: assert if caller passes in t == 0 / t == 1 ? +SkDPoint SkDQuad::ptAtT(double t) const { + if (0 == t) { + return fPts[0]; + } + if (1 == t) { + return fPts[2]; + } + double one_t = 1 - t; + double a = one_t * one_t; + double b = 2 * one_t * t; + double c = t * t; + SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX, + a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY }; + return result; +} + +/* +Given a quadratic q, t1, and t2, find a small quadratic segment. + +The new quadratic is defined by A, B, and C, where + A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1 + C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1 + +To find B, compute the point halfway between t1 and t2: + +q(at (t1 + t2)/2) == D + +Next, compute where D must be if we know the value of B: + +_12 = A/2 + B/2 +12_ = B/2 + C/2 +123 = A/4 + B/2 + C/4 + = D + +Group the known values on one side: + +B = D*2 - A/2 - C/2 +*/ + +static double interp_quad_coords(const double* src, double t) { + double ab = SkDInterp(src[0], src[2], t); + double bc = SkDInterp(src[2], src[4], t); + double abc = SkDInterp(ab, bc, t); + return abc; +} + +bool SkDQuad::monotonicInY() const { + return between(fPts[0].fY, fPts[1].fY, fPts[2].fY); +} + +SkDQuad SkDQuad::subDivide(double t1, double t2) const { + SkDQuad dst; + double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1); + double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1); + double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2); + double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2); + double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2); + double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2); + /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2; + /* by = */ dst[1].fY = 2*dy - (ay + cy)/2; + return dst; +} + +void SkDQuad::align(int endIndex, SkDPoint* dstPt) const { + if (fPts[endIndex].fX == fPts[1].fX) { + dstPt->fX = fPts[endIndex].fX; + } + if (fPts[endIndex].fY == fPts[1].fY) { + dstPt->fY = fPts[endIndex].fY; + } +} + +SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const { + SkASSERT(t1 != t2); + SkDPoint b; +#if 0 + // this approach assumes that the control point computed directly is accurate enough + double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2); + double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2); + b.fX = 2 * dx - (a.fX + c.fX) / 2; + b.fY = 2 * dy - (a.fY + c.fY) / 2; +#else + SkDQuad sub = subDivide(t1, t2); + SkDLine b0 = {{a, sub[1] + (a - sub[0])}}; + SkDLine b1 = {{c, sub[1] + (c - sub[2])}}; + SkIntersections i; + i.intersectRay(b0, b1); + if (i.used() == 1) { + b = i.pt(0); + } else { + SkASSERT(i.used() == 2 || i.used() == 0); + b = SkDPoint::Mid(b0[1], b1[1]); + } +#endif + if (t1 == 0 || t2 == 0) { + align(0, &b); + } + if (t1 == 1 || t2 == 1) { + align(2, &b); + } + if (precisely_subdivide_equal(b.fX, a.fX)) { + b.fX = a.fX; + } else if (precisely_subdivide_equal(b.fX, c.fX)) { + b.fX = c.fX; + } + if (precisely_subdivide_equal(b.fY, a.fY)) { + b.fY = a.fY; + } else if (precisely_subdivide_equal(b.fY, c.fY)) { + b.fY = c.fY; + } + return b; +} + +/* classic one t subdivision */ +static void interp_quad_coords(const double* src, double* dst, double t) { + double ab = SkDInterp(src[0], src[2], t); + double bc = SkDInterp(src[2], src[4], t); + dst[0] = src[0]; + dst[2] = ab; + dst[4] = SkDInterp(ab, bc, t); + dst[6] = bc; + dst[8] = src[4]; +} + +SkDQuadPair SkDQuad::chopAt(double t) const +{ + SkDQuadPair dst; + interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t); + interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t); + return dst; +} + +static int valid_unit_divide(double numer, double denom, double* ratio) +{ + if (numer < 0) { + numer = -numer; + denom = -denom; + } + if (denom == 0 || numer == 0 || numer >= denom) { + return 0; + } + double r = numer / denom; + if (r == 0) { // catch underflow if numer <<<< denom + return 0; + } + *ratio = r; + return 1; +} + +/** Quad'(t) = At + B, where + A = 2(a - 2b + c) + B = 2(b - a) + Solve for t, only if it fits between 0 < t < 1 +*/ +int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) { + /* At + B == 0 + t = -B / A + */ + return valid_unit_divide(a - b, a - b - b + c, tValue); +} + +/* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t) + * + * a = A - 2*B + C + * b = 2*B - 2*C + * c = C + */ +void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) { + *a = quad[0]; // a = A + *b = 2 * quad[2]; // b = 2*B + *c = quad[4]; // c = C + *b -= *c; // b = 2*B - C + *a -= *b; // a = A - 2*B + C + *b -= *c; // b = 2*B - 2*C +} + +#ifdef SK_DEBUG +void SkDQuad::dump() { + SkDebugf("{{"); + int index = 0; + do { + fPts[index].dump(); + SkDebugf(", "); + } while (++index < 2); + fPts[index].dump(); + SkDebugf("}}\n"); +} +#endif