1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/gfx/skia/trunk/src/pathops/SkDCubicLineIntersection.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,358 @@
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 "SkPathOpsCubic.h"
1.12 +#include "SkPathOpsLine.h"
1.13 +
1.14 +/*
1.15 +Find the interection of a line and cubic by solving for valid t values.
1.16 +
1.17 +Analogous to line-quadratic intersection, solve line-cubic intersection by
1.18 +representing the cubic as:
1.19 + x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
1.20 + y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
1.21 +and the line as:
1.22 + y = i*x + j (if the line is more horizontal)
1.23 +or:
1.24 + x = i*y + j (if the line is more vertical)
1.25 +
1.26 +Then using Mathematica, solve for the values of t where the cubic intersects the
1.27 +line:
1.28 +
1.29 + (in) Resultant[
1.30 + a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
1.31 + e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
1.32 + (out) -e + j +
1.33 + 3 e t - 3 f t -
1.34 + 3 e t^2 + 6 f t^2 - 3 g t^2 +
1.35 + e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
1.36 + i ( a -
1.37 + 3 a t + 3 b t +
1.38 + 3 a t^2 - 6 b t^2 + 3 c t^2 -
1.39 + a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
1.40 +
1.41 +if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
1.42 +
1.43 + (in) Resultant[
1.44 + a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
1.45 + e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
1.46 + (out) a - j -
1.47 + 3 a t + 3 b t +
1.48 + 3 a t^2 - 6 b t^2 + 3 c t^2 -
1.49 + a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
1.50 + i ( e -
1.51 + 3 e t + 3 f t +
1.52 + 3 e t^2 - 6 f t^2 + 3 g t^2 -
1.53 + e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
1.54 +
1.55 +Solving this with Mathematica produces an expression with hundreds of terms;
1.56 +instead, use Numeric Solutions recipe to solve the cubic.
1.57 +
1.58 +The near-horizontal case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
1.59 + A = (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d) )
1.60 + B = 3*(-( e - 2*f + g ) + i*( a - 2*b + c ) )
1.61 + C = 3*(-(-e + f ) + i*(-a + b ) )
1.62 + D = (-( e ) + i*( a ) + j )
1.63 +
1.64 +The near-vertical case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
1.65 + A = ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h) )
1.66 + B = 3*( ( a - 2*b + c ) - i*( e - 2*f + g ) )
1.67 + C = 3*( (-a + b ) - i*(-e + f ) )
1.68 + D = ( ( a ) - i*( e ) - j )
1.69 +
1.70 +For horizontal lines:
1.71 +(in) Resultant[
1.72 + a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
1.73 + e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
1.74 +(out) e - j -
1.75 + 3 e t + 3 f t +
1.76 + 3 e t^2 - 6 f t^2 + 3 g t^2 -
1.77 + e t^3 + 3 f t^3 - 3 g t^3 + h t^3
1.78 + */
1.79 +
1.80 +class LineCubicIntersections {
1.81 +public:
1.82 + enum PinTPoint {
1.83 + kPointUninitialized,
1.84 + kPointInitialized
1.85 + };
1.86 +
1.87 + LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections* i)
1.88 + : fCubic(c)
1.89 + , fLine(l)
1.90 + , fIntersections(i)
1.91 + , fAllowNear(true) {
1.92 + i->setMax(3);
1.93 + }
1.94 +
1.95 + void allowNear(bool allow) {
1.96 + fAllowNear = allow;
1.97 + }
1.98 +
1.99 + // see parallel routine in line quadratic intersections
1.100 + int intersectRay(double roots[3]) {
1.101 + double adj = fLine[1].fX - fLine[0].fX;
1.102 + double opp = fLine[1].fY - fLine[0].fY;
1.103 + SkDCubic r;
1.104 + for (int n = 0; n < 4; ++n) {
1.105 + r[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp;
1.106 + }
1.107 + double A, B, C, D;
1.108 + SkDCubic::Coefficients(&r[0].fX, &A, &B, &C, &D);
1.109 + return SkDCubic::RootsValidT(A, B, C, D, roots);
1.110 + }
1.111 +
1.112 + int intersect() {
1.113 + addExactEndPoints();
1.114 + if (fAllowNear) {
1.115 + addNearEndPoints();
1.116 + }
1.117 + double rootVals[3];
1.118 + int roots = intersectRay(rootVals);
1.119 + for (int index = 0; index < roots; ++index) {
1.120 + double cubicT = rootVals[index];
1.121 + double lineT = findLineT(cubicT);
1.122 + SkDPoint pt;
1.123 + if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized)) {
1.124 + #if ONE_OFF_DEBUG
1.125 + SkDPoint cPt = fCubic.ptAtT(cubicT);
1.126 + SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
1.127 + cPt.fX, cPt.fY);
1.128 + #endif
1.129 + for (int inner = 0; inner < fIntersections->used(); ++inner) {
1.130 + if (fIntersections->pt(inner) != pt) {
1.131 + continue;
1.132 + }
1.133 + double existingCubicT = (*fIntersections)[0][inner];
1.134 + if (cubicT == existingCubicT) {
1.135 + goto skipInsert;
1.136 + }
1.137 + // check if midway on cubic is also same point. If so, discard this
1.138 + double cubicMidT = (existingCubicT + cubicT) / 2;
1.139 + SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
1.140 + if (cubicMidPt.approximatelyEqual(pt)) {
1.141 + goto skipInsert;
1.142 + }
1.143 + }
1.144 + fIntersections->insert(cubicT, lineT, pt);
1.145 + skipInsert:
1.146 + ;
1.147 + }
1.148 + }
1.149 + return fIntersections->used();
1.150 + }
1.151 +
1.152 + int horizontalIntersect(double axisIntercept, double roots[3]) {
1.153 + double A, B, C, D;
1.154 + SkDCubic::Coefficients(&fCubic[0].fY, &A, &B, &C, &D);
1.155 + D -= axisIntercept;
1.156 + return SkDCubic::RootsValidT(A, B, C, D, roots);
1.157 + }
1.158 +
1.159 + int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
1.160 + addExactHorizontalEndPoints(left, right, axisIntercept);
1.161 + if (fAllowNear) {
1.162 + addNearHorizontalEndPoints(left, right, axisIntercept);
1.163 + }
1.164 + double rootVals[3];
1.165 + int roots = horizontalIntersect(axisIntercept, rootVals);
1.166 + for (int index = 0; index < roots; ++index) {
1.167 + double cubicT = rootVals[index];
1.168 + SkDPoint pt = fCubic.ptAtT(cubicT);
1.169 + double lineT = (pt.fX - left) / (right - left);
1.170 + if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
1.171 + fIntersections->insert(cubicT, lineT, pt);
1.172 + }
1.173 + }
1.174 + if (flipped) {
1.175 + fIntersections->flip();
1.176 + }
1.177 + return fIntersections->used();
1.178 + }
1.179 +
1.180 + int verticalIntersect(double axisIntercept, double roots[3]) {
1.181 + double A, B, C, D;
1.182 + SkDCubic::Coefficients(&fCubic[0].fX, &A, &B, &C, &D);
1.183 + D -= axisIntercept;
1.184 + return SkDCubic::RootsValidT(A, B, C, D, roots);
1.185 + }
1.186 +
1.187 + int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
1.188 + addExactVerticalEndPoints(top, bottom, axisIntercept);
1.189 + if (fAllowNear) {
1.190 + addNearVerticalEndPoints(top, bottom, axisIntercept);
1.191 + }
1.192 + double rootVals[3];
1.193 + int roots = verticalIntersect(axisIntercept, rootVals);
1.194 + for (int index = 0; index < roots; ++index) {
1.195 + double cubicT = rootVals[index];
1.196 + SkDPoint pt = fCubic.ptAtT(cubicT);
1.197 + double lineT = (pt.fY - top) / (bottom - top);
1.198 + if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
1.199 + fIntersections->insert(cubicT, lineT, pt);
1.200 + }
1.201 + }
1.202 + if (flipped) {
1.203 + fIntersections->flip();
1.204 + }
1.205 + return fIntersections->used();
1.206 + }
1.207 +
1.208 + protected:
1.209 +
1.210 + void addExactEndPoints() {
1.211 + for (int cIndex = 0; cIndex < 4; cIndex += 3) {
1.212 + double lineT = fLine.exactPoint(fCubic[cIndex]);
1.213 + if (lineT < 0) {
1.214 + continue;
1.215 + }
1.216 + double cubicT = (double) (cIndex >> 1);
1.217 + fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
1.218 + }
1.219 + }
1.220 +
1.221 + /* Note that this does not look for endpoints of the line that are near the cubic.
1.222 + These points are found later when check ends looks for missing points */
1.223 + void addNearEndPoints() {
1.224 + for (int cIndex = 0; cIndex < 4; cIndex += 3) {
1.225 + double cubicT = (double) (cIndex >> 1);
1.226 + if (fIntersections->hasT(cubicT)) {
1.227 + continue;
1.228 + }
1.229 + double lineT = fLine.nearPoint(fCubic[cIndex]);
1.230 + if (lineT < 0) {
1.231 + continue;
1.232 + }
1.233 + fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
1.234 + }
1.235 + }
1.236 +
1.237 + void addExactHorizontalEndPoints(double left, double right, double y) {
1.238 + for (int cIndex = 0; cIndex < 4; cIndex += 3) {
1.239 + double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
1.240 + if (lineT < 0) {
1.241 + continue;
1.242 + }
1.243 + double cubicT = (double) (cIndex >> 1);
1.244 + fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
1.245 + }
1.246 + }
1.247 +
1.248 + void addNearHorizontalEndPoints(double left, double right, double y) {
1.249 + for (int cIndex = 0; cIndex < 4; cIndex += 3) {
1.250 + double cubicT = (double) (cIndex >> 1);
1.251 + if (fIntersections->hasT(cubicT)) {
1.252 + continue;
1.253 + }
1.254 + double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
1.255 + if (lineT < 0) {
1.256 + continue;
1.257 + }
1.258 + fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
1.259 + }
1.260 + // FIXME: see if line end is nearly on cubic
1.261 + }
1.262 +
1.263 + void addExactVerticalEndPoints(double top, double bottom, double x) {
1.264 + for (int cIndex = 0; cIndex < 4; cIndex += 3) {
1.265 + double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
1.266 + if (lineT < 0) {
1.267 + continue;
1.268 + }
1.269 + double cubicT = (double) (cIndex >> 1);
1.270 + fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
1.271 + }
1.272 + }
1.273 +
1.274 + void addNearVerticalEndPoints(double top, double bottom, double x) {
1.275 + for (int cIndex = 0; cIndex < 4; cIndex += 3) {
1.276 + double cubicT = (double) (cIndex >> 1);
1.277 + if (fIntersections->hasT(cubicT)) {
1.278 + continue;
1.279 + }
1.280 + double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
1.281 + if (lineT < 0) {
1.282 + continue;
1.283 + }
1.284 + fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
1.285 + }
1.286 + // FIXME: see if line end is nearly on cubic
1.287 + }
1.288 +
1.289 + double findLineT(double t) {
1.290 + SkDPoint xy = fCubic.ptAtT(t);
1.291 + double dx = fLine[1].fX - fLine[0].fX;
1.292 + double dy = fLine[1].fY - fLine[0].fY;
1.293 + if (fabs(dx) > fabs(dy)) {
1.294 + return (xy.fX - fLine[0].fX) / dx;
1.295 + }
1.296 + return (xy.fY - fLine[0].fY) / dy;
1.297 + }
1.298 +
1.299 + bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
1.300 + if (!approximately_one_or_less(*lineT)) {
1.301 + return false;
1.302 + }
1.303 + if (!approximately_zero_or_more(*lineT)) {
1.304 + return false;
1.305 + }
1.306 + double cT = *cubicT = SkPinT(*cubicT);
1.307 + double lT = *lineT = SkPinT(*lineT);
1.308 + if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
1.309 + *pt = fLine.ptAtT(lT);
1.310 + } else if (ptSet == kPointUninitialized) {
1.311 + *pt = fCubic.ptAtT(cT);
1.312 + }
1.313 + SkPoint gridPt = pt->asSkPoint();
1.314 + if (gridPt == fLine[0].asSkPoint()) {
1.315 + *lineT = 0;
1.316 + } else if (gridPt == fLine[1].asSkPoint()) {
1.317 + *lineT = 1;
1.318 + }
1.319 + if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
1.320 + *cubicT = 0;
1.321 + } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
1.322 + *cubicT = 1;
1.323 + }
1.324 + return true;
1.325 + }
1.326 +
1.327 +private:
1.328 + const SkDCubic& fCubic;
1.329 + const SkDLine& fLine;
1.330 + SkIntersections* fIntersections;
1.331 + bool fAllowNear;
1.332 +};
1.333 +
1.334 +int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
1.335 + bool flipped) {
1.336 + SkDLine line = {{{ left, y }, { right, y }}};
1.337 + LineCubicIntersections c(cubic, line, this);
1.338 + return c.horizontalIntersect(y, left, right, flipped);
1.339 +}
1.340 +
1.341 +int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
1.342 + bool flipped) {
1.343 + SkDLine line = {{{ x, top }, { x, bottom }}};
1.344 + LineCubicIntersections c(cubic, line, this);
1.345 + return c.verticalIntersect(x, top, bottom, flipped);
1.346 +}
1.347 +
1.348 +int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
1.349 + LineCubicIntersections c(cubic, line, this);
1.350 + c.allowNear(fAllowNear);
1.351 + return c.intersect();
1.352 +}
1.353 +
1.354 +int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
1.355 + LineCubicIntersections c(cubic, line, this);
1.356 + fUsed = c.intersectRay(fT[0]);
1.357 + for (int index = 0; index < fUsed; ++index) {
1.358 + fPt[index] = cubic.ptAtT(fT[0][index]);
1.359 + }
1.360 + return fUsed;
1.361 +}
