1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/gfx/skia/trunk/src/pathops/SkDQuadLineIntersection.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,413 @@
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 "SkPathOpsLine.h"
1.12 +#include "SkPathOpsQuad.h"
1.13 +
1.14 +/*
1.15 +Find the interection of a line and quadratic by solving for valid t values.
1.16 +
1.17 +From http://stackoverflow.com/questions/1853637/how-to-find-the-mathematical-function-defining-a-bezier-curve
1.18 +
1.19 +"A Bezier curve is a parametric function. A quadratic Bezier curve (i.e. three
1.20 +control points) can be expressed as: F(t) = A(1 - t)^2 + B(1 - t)t + Ct^2 where
1.21 +A, B and C are points and t goes from zero to one.
1.22 +
1.23 +This will give you two equations:
1.24 +
1.25 + x = a(1 - t)^2 + b(1 - t)t + ct^2
1.26 + y = d(1 - t)^2 + e(1 - t)t + ft^2
1.27 +
1.28 +If you add for instance the line equation (y = kx + m) to that, you'll end up
1.29 +with three equations and three unknowns (x, y and t)."
1.30 +
1.31 +Similar to above, the quadratic is represented as
1.32 + x = a(1-t)^2 + 2b(1-t)t + ct^2
1.33 + y = d(1-t)^2 + 2e(1-t)t + ft^2
1.34 +and the line as
1.35 + y = g*x + h
1.36 +
1.37 +Using Mathematica, solve for the values of t where the quadratic intersects the
1.38 +line:
1.39 +
1.40 + (in) t1 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - x,
1.41 + d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - g*x - h, x]
1.42 + (out) -d + h + 2 d t - 2 e t - d t^2 + 2 e t^2 - f t^2 +
1.43 + g (a - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2)
1.44 + (in) Solve[t1 == 0, t]
1.45 + (out) {
1.46 + {t -> (-2 d + 2 e + 2 a g - 2 b g -
1.47 + Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 -
1.48 + 4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) /
1.49 + (2 (-d + 2 e - f + a g - 2 b g + c g))
1.50 + },
1.51 + {t -> (-2 d + 2 e + 2 a g - 2 b g +
1.52 + Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 -
1.53 + 4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) /
1.54 + (2 (-d + 2 e - f + a g - 2 b g + c g))
1.55 + }
1.56 + }
1.57 +
1.58 +Using the results above (when the line tends towards horizontal)
1.59 + A = (-(d - 2*e + f) + g*(a - 2*b + c) )
1.60 + B = 2*( (d - e ) - g*(a - b ) )
1.61 + C = (-(d ) + g*(a ) + h )
1.62 +
1.63 +If g goes to infinity, we can rewrite the line in terms of x.
1.64 + x = g'*y + h'
1.65 +
1.66 +And solve accordingly in Mathematica:
1.67 +
1.68 + (in) t2 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - g'*y - h',
1.69 + d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - y, y]
1.70 + (out) a - h' - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2 -
1.71 + g' (d - 2 d t + 2 e t + d t^2 - 2 e t^2 + f t^2)
1.72 + (in) Solve[t2 == 0, t]
1.73 + (out) {
1.74 + {t -> (2 a - 2 b - 2 d g' + 2 e g' -
1.75 + Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 -
1.76 + 4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')]) /
1.77 + (2 (a - 2 b + c - d g' + 2 e g' - f g'))
1.78 + },
1.79 + {t -> (2 a - 2 b - 2 d g' + 2 e g' +
1.80 + Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 -
1.81 + 4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')])/
1.82 + (2 (a - 2 b + c - d g' + 2 e g' - f g'))
1.83 + }
1.84 + }
1.85 +
1.86 +Thus, if the slope of the line tends towards vertical, we use:
1.87 + A = ( (a - 2*b + c) - g'*(d - 2*e + f) )
1.88 + B = 2*(-(a - b ) + g'*(d - e ) )
1.89 + C = ( (a ) - g'*(d ) - h' )
1.90 + */
1.91 +
1.92 +
1.93 +class LineQuadraticIntersections {
1.94 +public:
1.95 + enum PinTPoint {
1.96 + kPointUninitialized,
1.97 + kPointInitialized
1.98 + };
1.99 +
1.100 + LineQuadraticIntersections(const SkDQuad& q, const SkDLine& l, SkIntersections* i)
1.101 + : fQuad(q)
1.102 + , fLine(l)
1.103 + , fIntersections(i)
1.104 + , fAllowNear(true) {
1.105 + i->setMax(2);
1.106 + }
1.107 +
1.108 + void allowNear(bool allow) {
1.109 + fAllowNear = allow;
1.110 + }
1.111 +
1.112 + int intersectRay(double roots[2]) {
1.113 + /*
1.114 + solve by rotating line+quad so line is horizontal, then finding the roots
1.115 + set up matrix to rotate quad to x-axis
1.116 + |cos(a) -sin(a)|
1.117 + |sin(a) cos(a)|
1.118 + note that cos(a) = A(djacent) / Hypoteneuse
1.119 + sin(a) = O(pposite) / Hypoteneuse
1.120 + since we are computing Ts, we can ignore hypoteneuse, the scale factor:
1.121 + | A -O |
1.122 + | O A |
1.123 + A = line[1].fX - line[0].fX (adjacent side of the right triangle)
1.124 + O = line[1].fY - line[0].fY (opposite side of the right triangle)
1.125 + for each of the three points (e.g. n = 0 to 2)
1.126 + quad[n].fY' = (quad[n].fY - line[0].fY) * A - (quad[n].fX - line[0].fX) * O
1.127 + */
1.128 + double adj = fLine[1].fX - fLine[0].fX;
1.129 + double opp = fLine[1].fY - fLine[0].fY;
1.130 + double r[3];
1.131 + for (int n = 0; n < 3; ++n) {
1.132 + r[n] = (fQuad[n].fY - fLine[0].fY) * adj - (fQuad[n].fX - fLine[0].fX) * opp;
1.133 + }
1.134 + double A = r[2];
1.135 + double B = r[1];
1.136 + double C = r[0];
1.137 + A += C - 2 * B; // A = a - 2*b + c
1.138 + B -= C; // B = -(b - c)
1.139 + return SkDQuad::RootsValidT(A, 2 * B, C, roots);
1.140 + }
1.141 +
1.142 + int intersect() {
1.143 + addExactEndPoints();
1.144 + if (fAllowNear) {
1.145 + addNearEndPoints();
1.146 + }
1.147 + if (fIntersections->used() == 2) {
1.148 + // FIXME : need sharable code that turns spans into coincident if middle point is on
1.149 + } else {
1.150 + double rootVals[2];
1.151 + int roots = intersectRay(rootVals);
1.152 + for (int index = 0; index < roots; ++index) {
1.153 + double quadT = rootVals[index];
1.154 + double lineT = findLineT(quadT);
1.155 + SkDPoint pt;
1.156 + if (pinTs(&quadT, &lineT, &pt, kPointUninitialized)) {
1.157 + fIntersections->insert(quadT, lineT, pt);
1.158 + }
1.159 + }
1.160 + }
1.161 + return fIntersections->used();
1.162 + }
1.163 +
1.164 + int horizontalIntersect(double axisIntercept, double roots[2]) {
1.165 + double D = fQuad[2].fY; // f
1.166 + double E = fQuad[1].fY; // e
1.167 + double F = fQuad[0].fY; // d
1.168 + D += F - 2 * E; // D = d - 2*e + f
1.169 + E -= F; // E = -(d - e)
1.170 + F -= axisIntercept;
1.171 + return SkDQuad::RootsValidT(D, 2 * E, F, roots);
1.172 + }
1.173 +
1.174 + int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
1.175 + addExactHorizontalEndPoints(left, right, axisIntercept);
1.176 + if (fAllowNear) {
1.177 + addNearHorizontalEndPoints(left, right, axisIntercept);
1.178 + }
1.179 + double rootVals[2];
1.180 + int roots = horizontalIntersect(axisIntercept, rootVals);
1.181 + for (int index = 0; index < roots; ++index) {
1.182 + double quadT = rootVals[index];
1.183 + SkDPoint pt = fQuad.ptAtT(quadT);
1.184 + double lineT = (pt.fX - left) / (right - left);
1.185 + if (pinTs(&quadT, &lineT, &pt, kPointInitialized)) {
1.186 + fIntersections->insert(quadT, lineT, pt);
1.187 + }
1.188 + }
1.189 + if (flipped) {
1.190 + fIntersections->flip();
1.191 + }
1.192 + return fIntersections->used();
1.193 + }
1.194 +
1.195 + int verticalIntersect(double axisIntercept, double roots[2]) {
1.196 + double D = fQuad[2].fX; // f
1.197 + double E = fQuad[1].fX; // e
1.198 + double F = fQuad[0].fX; // d
1.199 + D += F - 2 * E; // D = d - 2*e + f
1.200 + E -= F; // E = -(d - e)
1.201 + F -= axisIntercept;
1.202 + return SkDQuad::RootsValidT(D, 2 * E, F, roots);
1.203 + }
1.204 +
1.205 + int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
1.206 + addExactVerticalEndPoints(top, bottom, axisIntercept);
1.207 + if (fAllowNear) {
1.208 + addNearVerticalEndPoints(top, bottom, axisIntercept);
1.209 + }
1.210 + double rootVals[2];
1.211 + int roots = verticalIntersect(axisIntercept, rootVals);
1.212 + for (int index = 0; index < roots; ++index) {
1.213 + double quadT = rootVals[index];
1.214 + SkDPoint pt = fQuad.ptAtT(quadT);
1.215 + double lineT = (pt.fY - top) / (bottom - top);
1.216 + if (pinTs(&quadT, &lineT, &pt, kPointInitialized)) {
1.217 + fIntersections->insert(quadT, lineT, pt);
1.218 + }
1.219 + }
1.220 + if (flipped) {
1.221 + fIntersections->flip();
1.222 + }
1.223 + return fIntersections->used();
1.224 + }
1.225 +
1.226 +protected:
1.227 + // add endpoints first to get zero and one t values exactly
1.228 + void addExactEndPoints() {
1.229 + for (int qIndex = 0; qIndex < 3; qIndex += 2) {
1.230 + double lineT = fLine.exactPoint(fQuad[qIndex]);
1.231 + if (lineT < 0) {
1.232 + continue;
1.233 + }
1.234 + double quadT = (double) (qIndex >> 1);
1.235 + fIntersections->insert(quadT, lineT, fQuad[qIndex]);
1.236 + }
1.237 + }
1.238 +
1.239 + void addNearEndPoints() {
1.240 + for (int qIndex = 0; qIndex < 3; qIndex += 2) {
1.241 + double quadT = (double) (qIndex >> 1);
1.242 + if (fIntersections->hasT(quadT)) {
1.243 + continue;
1.244 + }
1.245 + double lineT = fLine.nearPoint(fQuad[qIndex]);
1.246 + if (lineT < 0) {
1.247 + continue;
1.248 + }
1.249 + fIntersections->insert(quadT, lineT, fQuad[qIndex]);
1.250 + }
1.251 + // FIXME: see if line end is nearly on quad
1.252 + }
1.253 +
1.254 + void addExactHorizontalEndPoints(double left, double right, double y) {
1.255 + for (int qIndex = 0; qIndex < 3; qIndex += 2) {
1.256 + double lineT = SkDLine::ExactPointH(fQuad[qIndex], left, right, y);
1.257 + if (lineT < 0) {
1.258 + continue;
1.259 + }
1.260 + double quadT = (double) (qIndex >> 1);
1.261 + fIntersections->insert(quadT, lineT, fQuad[qIndex]);
1.262 + }
1.263 + }
1.264 +
1.265 + void addNearHorizontalEndPoints(double left, double right, double y) {
1.266 + for (int qIndex = 0; qIndex < 3; qIndex += 2) {
1.267 + double quadT = (double) (qIndex >> 1);
1.268 + if (fIntersections->hasT(quadT)) {
1.269 + continue;
1.270 + }
1.271 + double lineT = SkDLine::NearPointH(fQuad[qIndex], left, right, y);
1.272 + if (lineT < 0) {
1.273 + continue;
1.274 + }
1.275 + fIntersections->insert(quadT, lineT, fQuad[qIndex]);
1.276 + }
1.277 + // FIXME: see if line end is nearly on quad
1.278 + }
1.279 +
1.280 + void addExactVerticalEndPoints(double top, double bottom, double x) {
1.281 + for (int qIndex = 0; qIndex < 3; qIndex += 2) {
1.282 + double lineT = SkDLine::ExactPointV(fQuad[qIndex], top, bottom, x);
1.283 + if (lineT < 0) {
1.284 + continue;
1.285 + }
1.286 + double quadT = (double) (qIndex >> 1);
1.287 + fIntersections->insert(quadT, lineT, fQuad[qIndex]);
1.288 + }
1.289 + }
1.290 +
1.291 + void addNearVerticalEndPoints(double top, double bottom, double x) {
1.292 + for (int qIndex = 0; qIndex < 3; qIndex += 2) {
1.293 + double quadT = (double) (qIndex >> 1);
1.294 + if (fIntersections->hasT(quadT)) {
1.295 + continue;
1.296 + }
1.297 + double lineT = SkDLine::NearPointV(fQuad[qIndex], top, bottom, x);
1.298 + if (lineT < 0) {
1.299 + continue;
1.300 + }
1.301 + fIntersections->insert(quadT, lineT, fQuad[qIndex]);
1.302 + }
1.303 + // FIXME: see if line end is nearly on quad
1.304 + }
1.305 +
1.306 + double findLineT(double t) {
1.307 + SkDPoint xy = fQuad.ptAtT(t);
1.308 + double dx = fLine[1].fX - fLine[0].fX;
1.309 + double dy = fLine[1].fY - fLine[0].fY;
1.310 + if (fabs(dx) > fabs(dy)) {
1.311 + return (xy.fX - fLine[0].fX) / dx;
1.312 + }
1.313 + return (xy.fY - fLine[0].fY) / dy;
1.314 + }
1.315 +
1.316 + bool pinTs(double* quadT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
1.317 + if (!approximately_one_or_less(*lineT)) {
1.318 + return false;
1.319 + }
1.320 + if (!approximately_zero_or_more(*lineT)) {
1.321 + return false;
1.322 + }
1.323 + double qT = *quadT = SkPinT(*quadT);
1.324 + double lT = *lineT = SkPinT(*lineT);
1.325 + if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && qT != 0 && qT != 1)) {
1.326 + *pt = fLine.ptAtT(lT);
1.327 + } else if (ptSet == kPointUninitialized) {
1.328 + *pt = fQuad.ptAtT(qT);
1.329 + }
1.330 + SkPoint gridPt = pt->asSkPoint();
1.331 + if (gridPt == fLine[0].asSkPoint()) {
1.332 + *lineT = 0;
1.333 + } else if (gridPt == fLine[1].asSkPoint()) {
1.334 + *lineT = 1;
1.335 + }
1.336 + if (gridPt == fQuad[0].asSkPoint()) {
1.337 + *quadT = 0;
1.338 + } else if (gridPt == fQuad[2].asSkPoint()) {
1.339 + *quadT = 1;
1.340 + }
1.341 + return true;
1.342 + }
1.343 +
1.344 +private:
1.345 + const SkDQuad& fQuad;
1.346 + const SkDLine& fLine;
1.347 + SkIntersections* fIntersections;
1.348 + bool fAllowNear;
1.349 +};
1.350 +
1.351 +// utility for pairs of coincident quads
1.352 +static double horizontalIntersect(const SkDQuad& quad, const SkDPoint& pt) {
1.353 + LineQuadraticIntersections q(quad, *(static_cast<SkDLine*>(0)),
1.354 + static_cast<SkIntersections*>(0));
1.355 + double rootVals[2];
1.356 + int roots = q.horizontalIntersect(pt.fY, rootVals);
1.357 + for (int index = 0; index < roots; ++index) {
1.358 + double t = rootVals[index];
1.359 + SkDPoint qPt = quad.ptAtT(t);
1.360 + if (AlmostEqualUlps(qPt.fX, pt.fX)) {
1.361 + return t;
1.362 + }
1.363 + }
1.364 + return -1;
1.365 +}
1.366 +
1.367 +static double verticalIntersect(const SkDQuad& quad, const SkDPoint& pt) {
1.368 + LineQuadraticIntersections q(quad, *(static_cast<SkDLine*>(0)),
1.369 + static_cast<SkIntersections*>(0));
1.370 + double rootVals[2];
1.371 + int roots = q.verticalIntersect(pt.fX, rootVals);
1.372 + for (int index = 0; index < roots; ++index) {
1.373 + double t = rootVals[index];
1.374 + SkDPoint qPt = quad.ptAtT(t);
1.375 + if (AlmostEqualUlps(qPt.fY, pt.fY)) {
1.376 + return t;
1.377 + }
1.378 + }
1.379 + return -1;
1.380 +}
1.381 +
1.382 +double SkIntersections::Axial(const SkDQuad& q1, const SkDPoint& p, bool vertical) {
1.383 + if (vertical) {
1.384 + return verticalIntersect(q1, p);
1.385 + }
1.386 + return horizontalIntersect(q1, p);
1.387 +}
1.388 +
1.389 +int SkIntersections::horizontal(const SkDQuad& quad, double left, double right, double y,
1.390 + bool flipped) {
1.391 + SkDLine line = {{{ left, y }, { right, y }}};
1.392 + LineQuadraticIntersections q(quad, line, this);
1.393 + return q.horizontalIntersect(y, left, right, flipped);
1.394 +}
1.395 +
1.396 +int SkIntersections::vertical(const SkDQuad& quad, double top, double bottom, double x,
1.397 + bool flipped) {
1.398 + SkDLine line = {{{ x, top }, { x, bottom }}};
1.399 + LineQuadraticIntersections q(quad, line, this);
1.400 + return q.verticalIntersect(x, top, bottom, flipped);
1.401 +}
1.402 +
1.403 +int SkIntersections::intersect(const SkDQuad& quad, const SkDLine& line) {
1.404 + LineQuadraticIntersections q(quad, line, this);
1.405 + q.allowNear(fAllowNear);
1.406 + return q.intersect();
1.407 +}
1.408 +
1.409 +int SkIntersections::intersectRay(const SkDQuad& quad, const SkDLine& line) {
1.410 + LineQuadraticIntersections q(quad, line, this);
1.411 + fUsed = q.intersectRay(fT[0]);
1.412 + for (int index = 0; index < fUsed; ++index) {
1.413 + fPt[index] = quad.ptAtT(fT[0][index]);
1.414 + }
1.415 + return fUsed;
1.416 +}
