Wed, 31 Dec 2014 06:09:35 +0100
Cloned upstream origin tor-browser at tor-browser-31.3.0esr-4.5-1-build1
revision ID fc1c9ff7c1b2defdbc039f12214767608f46423f for hacking purpose.
michael@0 | 1 | /* |
michael@0 | 2 | * Copyright 2012 Google Inc. |
michael@0 | 3 | * |
michael@0 | 4 | * Use of this source code is governed by a BSD-style license that can be |
michael@0 | 5 | * found in the LICENSE file. |
michael@0 | 6 | */ |
michael@0 | 7 | #include "SkIntersections.h" |
michael@0 | 8 | #include "SkLineParameters.h" |
michael@0 | 9 | #include "SkPathOpsCubic.h" |
michael@0 | 10 | #include "SkPathOpsQuad.h" |
michael@0 | 11 | #include "SkPathOpsTriangle.h" |
michael@0 | 12 | |
michael@0 | 13 | // from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html |
michael@0 | 14 | // (currently only used by testing) |
michael@0 | 15 | double SkDQuad::nearestT(const SkDPoint& pt) const { |
michael@0 | 16 | SkDVector pos = fPts[0] - pt; |
michael@0 | 17 | // search points P of bezier curve with PM.(dP / dt) = 0 |
michael@0 | 18 | // a calculus leads to a 3d degree equation : |
michael@0 | 19 | SkDVector A = fPts[1] - fPts[0]; |
michael@0 | 20 | SkDVector B = fPts[2] - fPts[1]; |
michael@0 | 21 | B -= A; |
michael@0 | 22 | double a = B.dot(B); |
michael@0 | 23 | double b = 3 * A.dot(B); |
michael@0 | 24 | double c = 2 * A.dot(A) + pos.dot(B); |
michael@0 | 25 | double d = pos.dot(A); |
michael@0 | 26 | double ts[3]; |
michael@0 | 27 | int roots = SkDCubic::RootsValidT(a, b, c, d, ts); |
michael@0 | 28 | double d0 = pt.distanceSquared(fPts[0]); |
michael@0 | 29 | double d2 = pt.distanceSquared(fPts[2]); |
michael@0 | 30 | double distMin = SkTMin(d0, d2); |
michael@0 | 31 | int bestIndex = -1; |
michael@0 | 32 | for (int index = 0; index < roots; ++index) { |
michael@0 | 33 | SkDPoint onQuad = ptAtT(ts[index]); |
michael@0 | 34 | double dist = pt.distanceSquared(onQuad); |
michael@0 | 35 | if (distMin > dist) { |
michael@0 | 36 | distMin = dist; |
michael@0 | 37 | bestIndex = index; |
michael@0 | 38 | } |
michael@0 | 39 | } |
michael@0 | 40 | if (bestIndex >= 0) { |
michael@0 | 41 | return ts[bestIndex]; |
michael@0 | 42 | } |
michael@0 | 43 | return d0 < d2 ? 0 : 1; |
michael@0 | 44 | } |
michael@0 | 45 | |
michael@0 | 46 | bool SkDQuad::pointInHull(const SkDPoint& pt) const { |
michael@0 | 47 | return ((const SkDTriangle&) fPts).contains(pt); |
michael@0 | 48 | } |
michael@0 | 49 | |
michael@0 | 50 | SkDPoint SkDQuad::top(double startT, double endT) const { |
michael@0 | 51 | SkDQuad sub = subDivide(startT, endT); |
michael@0 | 52 | SkDPoint topPt = sub[0]; |
michael@0 | 53 | if (topPt.fY > sub[2].fY || (topPt.fY == sub[2].fY && topPt.fX > sub[2].fX)) { |
michael@0 | 54 | topPt = sub[2]; |
michael@0 | 55 | } |
michael@0 | 56 | if (!between(sub[0].fY, sub[1].fY, sub[2].fY)) { |
michael@0 | 57 | double extremeT; |
michael@0 | 58 | if (FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, &extremeT)) { |
michael@0 | 59 | extremeT = startT + (endT - startT) * extremeT; |
michael@0 | 60 | SkDPoint test = ptAtT(extremeT); |
michael@0 | 61 | if (topPt.fY > test.fY || (topPt.fY == test.fY && topPt.fX > test.fX)) { |
michael@0 | 62 | topPt = test; |
michael@0 | 63 | } |
michael@0 | 64 | } |
michael@0 | 65 | } |
michael@0 | 66 | return topPt; |
michael@0 | 67 | } |
michael@0 | 68 | |
michael@0 | 69 | int SkDQuad::AddValidTs(double s[], int realRoots, double* t) { |
michael@0 | 70 | int foundRoots = 0; |
michael@0 | 71 | for (int index = 0; index < realRoots; ++index) { |
michael@0 | 72 | double tValue = s[index]; |
michael@0 | 73 | if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) { |
michael@0 | 74 | if (approximately_less_than_zero(tValue)) { |
michael@0 | 75 | tValue = 0; |
michael@0 | 76 | } else if (approximately_greater_than_one(tValue)) { |
michael@0 | 77 | tValue = 1; |
michael@0 | 78 | } |
michael@0 | 79 | for (int idx2 = 0; idx2 < foundRoots; ++idx2) { |
michael@0 | 80 | if (approximately_equal(t[idx2], tValue)) { |
michael@0 | 81 | goto nextRoot; |
michael@0 | 82 | } |
michael@0 | 83 | } |
michael@0 | 84 | t[foundRoots++] = tValue; |
michael@0 | 85 | } |
michael@0 | 86 | nextRoot: |
michael@0 | 87 | {} |
michael@0 | 88 | } |
michael@0 | 89 | return foundRoots; |
michael@0 | 90 | } |
michael@0 | 91 | |
michael@0 | 92 | // note: caller expects multiple results to be sorted smaller first |
michael@0 | 93 | // note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting |
michael@0 | 94 | // analysis of the quadratic equation, suggesting why the following looks at |
michael@0 | 95 | // the sign of B -- and further suggesting that the greatest loss of precision |
michael@0 | 96 | // is in b squared less two a c |
michael@0 | 97 | int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) { |
michael@0 | 98 | double s[2]; |
michael@0 | 99 | int realRoots = RootsReal(A, B, C, s); |
michael@0 | 100 | int foundRoots = AddValidTs(s, realRoots, t); |
michael@0 | 101 | return foundRoots; |
michael@0 | 102 | } |
michael@0 | 103 | |
michael@0 | 104 | /* |
michael@0 | 105 | Numeric Solutions (5.6) suggests to solve the quadratic by computing |
michael@0 | 106 | Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C)) |
michael@0 | 107 | and using the roots |
michael@0 | 108 | t1 = Q / A |
michael@0 | 109 | t2 = C / Q |
michael@0 | 110 | */ |
michael@0 | 111 | // this does not discard real roots <= 0 or >= 1 |
michael@0 | 112 | int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) { |
michael@0 | 113 | const double p = B / (2 * A); |
michael@0 | 114 | const double q = C / A; |
michael@0 | 115 | if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) { |
michael@0 | 116 | if (approximately_zero(B)) { |
michael@0 | 117 | s[0] = 0; |
michael@0 | 118 | return C == 0; |
michael@0 | 119 | } |
michael@0 | 120 | s[0] = -C / B; |
michael@0 | 121 | return 1; |
michael@0 | 122 | } |
michael@0 | 123 | /* normal form: x^2 + px + q = 0 */ |
michael@0 | 124 | const double p2 = p * p; |
michael@0 | 125 | if (!AlmostDequalUlps(p2, q) && p2 < q) { |
michael@0 | 126 | return 0; |
michael@0 | 127 | } |
michael@0 | 128 | double sqrt_D = 0; |
michael@0 | 129 | if (p2 > q) { |
michael@0 | 130 | sqrt_D = sqrt(p2 - q); |
michael@0 | 131 | } |
michael@0 | 132 | s[0] = sqrt_D - p; |
michael@0 | 133 | s[1] = -sqrt_D - p; |
michael@0 | 134 | return 1 + !AlmostDequalUlps(s[0], s[1]); |
michael@0 | 135 | } |
michael@0 | 136 | |
michael@0 | 137 | bool SkDQuad::isLinear(int startIndex, int endIndex) const { |
michael@0 | 138 | SkLineParameters lineParameters; |
michael@0 | 139 | lineParameters.quadEndPoints(*this, startIndex, endIndex); |
michael@0 | 140 | // FIXME: maybe it's possible to avoid this and compare non-normalized |
michael@0 | 141 | lineParameters.normalize(); |
michael@0 | 142 | double distance = lineParameters.controlPtDistance(*this); |
michael@0 | 143 | return approximately_zero(distance); |
michael@0 | 144 | } |
michael@0 | 145 | |
michael@0 | 146 | SkDCubic SkDQuad::toCubic() const { |
michael@0 | 147 | SkDCubic cubic; |
michael@0 | 148 | cubic[0] = fPts[0]; |
michael@0 | 149 | cubic[2] = fPts[1]; |
michael@0 | 150 | cubic[3] = fPts[2]; |
michael@0 | 151 | cubic[1].fX = (cubic[0].fX + cubic[2].fX * 2) / 3; |
michael@0 | 152 | cubic[1].fY = (cubic[0].fY + cubic[2].fY * 2) / 3; |
michael@0 | 153 | cubic[2].fX = (cubic[3].fX + cubic[2].fX * 2) / 3; |
michael@0 | 154 | cubic[2].fY = (cubic[3].fY + cubic[2].fY * 2) / 3; |
michael@0 | 155 | return cubic; |
michael@0 | 156 | } |
michael@0 | 157 | |
michael@0 | 158 | SkDVector SkDQuad::dxdyAtT(double t) const { |
michael@0 | 159 | double a = t - 1; |
michael@0 | 160 | double b = 1 - 2 * t; |
michael@0 | 161 | double c = t; |
michael@0 | 162 | SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX, |
michael@0 | 163 | a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY }; |
michael@0 | 164 | return result; |
michael@0 | 165 | } |
michael@0 | 166 | |
michael@0 | 167 | // OPTIMIZE: assert if caller passes in t == 0 / t == 1 ? |
michael@0 | 168 | SkDPoint SkDQuad::ptAtT(double t) const { |
michael@0 | 169 | if (0 == t) { |
michael@0 | 170 | return fPts[0]; |
michael@0 | 171 | } |
michael@0 | 172 | if (1 == t) { |
michael@0 | 173 | return fPts[2]; |
michael@0 | 174 | } |
michael@0 | 175 | double one_t = 1 - t; |
michael@0 | 176 | double a = one_t * one_t; |
michael@0 | 177 | double b = 2 * one_t * t; |
michael@0 | 178 | double c = t * t; |
michael@0 | 179 | SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX, |
michael@0 | 180 | a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY }; |
michael@0 | 181 | return result; |
michael@0 | 182 | } |
michael@0 | 183 | |
michael@0 | 184 | /* |
michael@0 | 185 | Given a quadratic q, t1, and t2, find a small quadratic segment. |
michael@0 | 186 | |
michael@0 | 187 | The new quadratic is defined by A, B, and C, where |
michael@0 | 188 | A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1 |
michael@0 | 189 | C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1 |
michael@0 | 190 | |
michael@0 | 191 | To find B, compute the point halfway between t1 and t2: |
michael@0 | 192 | |
michael@0 | 193 | q(at (t1 + t2)/2) == D |
michael@0 | 194 | |
michael@0 | 195 | Next, compute where D must be if we know the value of B: |
michael@0 | 196 | |
michael@0 | 197 | _12 = A/2 + B/2 |
michael@0 | 198 | 12_ = B/2 + C/2 |
michael@0 | 199 | 123 = A/4 + B/2 + C/4 |
michael@0 | 200 | = D |
michael@0 | 201 | |
michael@0 | 202 | Group the known values on one side: |
michael@0 | 203 | |
michael@0 | 204 | B = D*2 - A/2 - C/2 |
michael@0 | 205 | */ |
michael@0 | 206 | |
michael@0 | 207 | static double interp_quad_coords(const double* src, double t) { |
michael@0 | 208 | double ab = SkDInterp(src[0], src[2], t); |
michael@0 | 209 | double bc = SkDInterp(src[2], src[4], t); |
michael@0 | 210 | double abc = SkDInterp(ab, bc, t); |
michael@0 | 211 | return abc; |
michael@0 | 212 | } |
michael@0 | 213 | |
michael@0 | 214 | bool SkDQuad::monotonicInY() const { |
michael@0 | 215 | return between(fPts[0].fY, fPts[1].fY, fPts[2].fY); |
michael@0 | 216 | } |
michael@0 | 217 | |
michael@0 | 218 | SkDQuad SkDQuad::subDivide(double t1, double t2) const { |
michael@0 | 219 | SkDQuad dst; |
michael@0 | 220 | double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1); |
michael@0 | 221 | double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1); |
michael@0 | 222 | double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2); |
michael@0 | 223 | double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2); |
michael@0 | 224 | double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2); |
michael@0 | 225 | double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2); |
michael@0 | 226 | /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2; |
michael@0 | 227 | /* by = */ dst[1].fY = 2*dy - (ay + cy)/2; |
michael@0 | 228 | return dst; |
michael@0 | 229 | } |
michael@0 | 230 | |
michael@0 | 231 | void SkDQuad::align(int endIndex, SkDPoint* dstPt) const { |
michael@0 | 232 | if (fPts[endIndex].fX == fPts[1].fX) { |
michael@0 | 233 | dstPt->fX = fPts[endIndex].fX; |
michael@0 | 234 | } |
michael@0 | 235 | if (fPts[endIndex].fY == fPts[1].fY) { |
michael@0 | 236 | dstPt->fY = fPts[endIndex].fY; |
michael@0 | 237 | } |
michael@0 | 238 | } |
michael@0 | 239 | |
michael@0 | 240 | SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const { |
michael@0 | 241 | SkASSERT(t1 != t2); |
michael@0 | 242 | SkDPoint b; |
michael@0 | 243 | #if 0 |
michael@0 | 244 | // this approach assumes that the control point computed directly is accurate enough |
michael@0 | 245 | double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2); |
michael@0 | 246 | double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2); |
michael@0 | 247 | b.fX = 2 * dx - (a.fX + c.fX) / 2; |
michael@0 | 248 | b.fY = 2 * dy - (a.fY + c.fY) / 2; |
michael@0 | 249 | #else |
michael@0 | 250 | SkDQuad sub = subDivide(t1, t2); |
michael@0 | 251 | SkDLine b0 = {{a, sub[1] + (a - sub[0])}}; |
michael@0 | 252 | SkDLine b1 = {{c, sub[1] + (c - sub[2])}}; |
michael@0 | 253 | SkIntersections i; |
michael@0 | 254 | i.intersectRay(b0, b1); |
michael@0 | 255 | if (i.used() == 1) { |
michael@0 | 256 | b = i.pt(0); |
michael@0 | 257 | } else { |
michael@0 | 258 | SkASSERT(i.used() == 2 || i.used() == 0); |
michael@0 | 259 | b = SkDPoint::Mid(b0[1], b1[1]); |
michael@0 | 260 | } |
michael@0 | 261 | #endif |
michael@0 | 262 | if (t1 == 0 || t2 == 0) { |
michael@0 | 263 | align(0, &b); |
michael@0 | 264 | } |
michael@0 | 265 | if (t1 == 1 || t2 == 1) { |
michael@0 | 266 | align(2, &b); |
michael@0 | 267 | } |
michael@0 | 268 | if (precisely_subdivide_equal(b.fX, a.fX)) { |
michael@0 | 269 | b.fX = a.fX; |
michael@0 | 270 | } else if (precisely_subdivide_equal(b.fX, c.fX)) { |
michael@0 | 271 | b.fX = c.fX; |
michael@0 | 272 | } |
michael@0 | 273 | if (precisely_subdivide_equal(b.fY, a.fY)) { |
michael@0 | 274 | b.fY = a.fY; |
michael@0 | 275 | } else if (precisely_subdivide_equal(b.fY, c.fY)) { |
michael@0 | 276 | b.fY = c.fY; |
michael@0 | 277 | } |
michael@0 | 278 | return b; |
michael@0 | 279 | } |
michael@0 | 280 | |
michael@0 | 281 | /* classic one t subdivision */ |
michael@0 | 282 | static void interp_quad_coords(const double* src, double* dst, double t) { |
michael@0 | 283 | double ab = SkDInterp(src[0], src[2], t); |
michael@0 | 284 | double bc = SkDInterp(src[2], src[4], t); |
michael@0 | 285 | dst[0] = src[0]; |
michael@0 | 286 | dst[2] = ab; |
michael@0 | 287 | dst[4] = SkDInterp(ab, bc, t); |
michael@0 | 288 | dst[6] = bc; |
michael@0 | 289 | dst[8] = src[4]; |
michael@0 | 290 | } |
michael@0 | 291 | |
michael@0 | 292 | SkDQuadPair SkDQuad::chopAt(double t) const |
michael@0 | 293 | { |
michael@0 | 294 | SkDQuadPair dst; |
michael@0 | 295 | interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t); |
michael@0 | 296 | interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t); |
michael@0 | 297 | return dst; |
michael@0 | 298 | } |
michael@0 | 299 | |
michael@0 | 300 | static int valid_unit_divide(double numer, double denom, double* ratio) |
michael@0 | 301 | { |
michael@0 | 302 | if (numer < 0) { |
michael@0 | 303 | numer = -numer; |
michael@0 | 304 | denom = -denom; |
michael@0 | 305 | } |
michael@0 | 306 | if (denom == 0 || numer == 0 || numer >= denom) { |
michael@0 | 307 | return 0; |
michael@0 | 308 | } |
michael@0 | 309 | double r = numer / denom; |
michael@0 | 310 | if (r == 0) { // catch underflow if numer <<<< denom |
michael@0 | 311 | return 0; |
michael@0 | 312 | } |
michael@0 | 313 | *ratio = r; |
michael@0 | 314 | return 1; |
michael@0 | 315 | } |
michael@0 | 316 | |
michael@0 | 317 | /** Quad'(t) = At + B, where |
michael@0 | 318 | A = 2(a - 2b + c) |
michael@0 | 319 | B = 2(b - a) |
michael@0 | 320 | Solve for t, only if it fits between 0 < t < 1 |
michael@0 | 321 | */ |
michael@0 | 322 | int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) { |
michael@0 | 323 | /* At + B == 0 |
michael@0 | 324 | t = -B / A |
michael@0 | 325 | */ |
michael@0 | 326 | return valid_unit_divide(a - b, a - b - b + c, tValue); |
michael@0 | 327 | } |
michael@0 | 328 | |
michael@0 | 329 | /* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t) |
michael@0 | 330 | * |
michael@0 | 331 | * a = A - 2*B + C |
michael@0 | 332 | * b = 2*B - 2*C |
michael@0 | 333 | * c = C |
michael@0 | 334 | */ |
michael@0 | 335 | void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) { |
michael@0 | 336 | *a = quad[0]; // a = A |
michael@0 | 337 | *b = 2 * quad[2]; // b = 2*B |
michael@0 | 338 | *c = quad[4]; // c = C |
michael@0 | 339 | *b -= *c; // b = 2*B - C |
michael@0 | 340 | *a -= *b; // a = A - 2*B + C |
michael@0 | 341 | *b -= *c; // b = 2*B - 2*C |
michael@0 | 342 | } |
michael@0 | 343 | |
michael@0 | 344 | #ifdef SK_DEBUG |
michael@0 | 345 | void SkDQuad::dump() { |
michael@0 | 346 | SkDebugf("{{"); |
michael@0 | 347 | int index = 0; |
michael@0 | 348 | do { |
michael@0 | 349 | fPts[index].dump(); |
michael@0 | 350 | SkDebugf(", "); |
michael@0 | 351 | } while (++index < 2); |
michael@0 | 352 | fPts[index].dump(); |
michael@0 | 353 | SkDebugf("}}\n"); |
michael@0 | 354 | } |
michael@0 | 355 | #endif |