| /* |
| * 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 = xyAtT(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 = xyAtT(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 (!AlmostEqualUlps(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 + !AlmostEqualUlps(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::xyAtT(double t) const { |
| 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 |
| } |