gfx/skia/trunk/src/pathops/SkPathOpsQuad.cpp

Thu, 22 Jan 2015 13:21:57 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Thu, 22 Jan 2015 13:21:57 +0100
branch
TOR_BUG_9701
changeset 15
b8a032363ba2
permissions
-rw-r--r--

Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6

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

mercurial