src/third_party/skia/src/pathops/SkDCubicLineIntersection.cpp

   1 /*
   2  * Copyright 2012 Google Inc.
   3  *
   4  * Use of this source code is governed by a BSD-style license that can be
   5  * found in the LICENSE file.
   6  */
   7 #include "SkIntersections.h"
   8 #include "SkPathOpsCubic.h"
   9 #include "SkPathOpsLine.h"
  10
  11 /*
  12 Find the interection of a line and cubic by solving for valid t values.
  13
  14 Analogous to line-quadratic intersection, solve line-cubic intersection by
  15 representing the cubic as:
  16   x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
  17   y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
  18 and the line as:
  19   y = i*x + j  (if the line is more horizontal)
  20 or:
  21   x = i*y + j  (if the line is more vertical)
  22
  23 Then using Mathematica, solve for the values of t where the cubic intersects the
  24 line:
  25
  26   (in) Resultant[
  27         a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
  28         e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
  29   (out) -e     +   j     +
  30        3 e t   - 3 f t   -
  31        3 e t^2 + 6 f t^2 - 3 g t^2 +
  32          e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
  33      i ( a     -
  34        3 a t + 3 b t +
  35        3 a t^2 - 6 b t^2 + 3 c t^2 -
  36          a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
  37
  38 if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
  39
  40   (in) Resultant[
  41         a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
  42         e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y,       y]
  43   (out)  a     -   j     -
  44        3 a t   + 3 b t   +
  45        3 a t^2 - 6 b t^2 + 3 c t^2 -
  46          a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
  47      i ( e     -
  48        3 e t   + 3 f t   +
  49        3 e t^2 - 6 f t^2 + 3 g t^2 -
  50          e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
  51
  52 Solving this with Mathematica produces an expression with hundreds of terms;
  53 instead, use Numeric Solutions recipe to solve the cubic.
  54
  55 The near-horizontal case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
  56     A =   (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d)     )
  57     B = 3*(-( e - 2*f +   g    ) + i*( a - 2*b +   c    )     )
  58     C = 3*(-(-e +   f          ) + i*(-a +   b          )     )
  59     D =   (-( e                ) + i*( a                ) + j )
  60
  61 The near-vertical case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
  62     A =   ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h)     )
  63     B = 3*( ( a - 2*b +   c    ) - i*( e - 2*f +   g    )     )
  64     C = 3*( (-a +   b          ) - i*(-e +   f          )     )
  65     D =   ( ( a                ) - i*( e                ) - j )
  66
  67 For horizontal lines:
  68 (in) Resultant[
  69       a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
  70       e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
  71 (out)  e     -   j     -
  72      3 e t   + 3 f t   +
  73      3 e t^2 - 6 f t^2 + 3 g t^2 -
  74        e t^3 + 3 f t^3 - 3 g t^3 + h t^3
  75  */
  76
  77 class LineCubicIntersections {
  78 public:
  79     enum PinTPoint {
  80         kPointUninitialized,
  81         kPointInitialized
  82     };
  83
  84     LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections* i)
  85         : fCubic(c)
  86         , fLine(l)
  87         , fIntersections(i)
  88         , fAllowNear(true) {
  89         i->setMax(3);
  90     }
  91
  92     void allowNear(bool allow) {
  93         fAllowNear = allow;
  94     }
  95
  96     // see parallel routine in line quadratic intersections
  97     int intersectRay(double roots[3]) {
  98         double adj = fLine[1].fX - fLine[0].fX;
  99         double opp = fLine[1].fY - fLine[0].fY;
 100         SkDCubic r;
 101         for (int n = 0; n < 4; ++n) {
 102             r[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp;
 103         }
 104         double A, B, C, D;
 105         SkDCubic::Coefficients(&r[0].fX, &A, &B, &C, &D);
 106         return SkDCubic::RootsValidT(A, B, C, D, roots);
 107     }
 108
 109     int intersect() {
 110         addExactEndPoints();
 111         if (fAllowNear) {
 112             addNearEndPoints();
 113         }
 114         double rootVals[3];
 115         int roots = intersectRay(rootVals);
 116         for (int index = 0; index < roots; ++index) {
 117             double cubicT = rootVals[index];
 118             double lineT = findLineT(cubicT);
 119             SkDPoint pt;
 120             if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized)) {
 121     #if ONE_OFF_DEBUG
 122                 SkDPoint cPt = fCubic.ptAtT(cubicT);
 123                 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
 124                         cPt.fX, cPt.fY);
 125     #endif
 126                 for (int inner = 0; inner < fIntersections->used(); ++inner) {
 127                     if (fIntersections->pt(inner) != pt) {
 128                         continue;
 129                     }
 130                     double existingCubicT = (*fIntersections)[0][inner];
 131                     if (cubicT == existingCubicT) {
 132                         goto skipInsert;
 133                     }
 134                     // check if midway on cubic is also same point. If so, discard this
 135                     double cubicMidT = (existingCubicT + cubicT) / 2;
 136                     SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
 137                     if (cubicMidPt.approximatelyEqual(pt)) {
 138                         goto skipInsert;
 139                     }
 140                 }
 141                 fIntersections->insert(cubicT, lineT, pt);
 142         skipInsert:
 143                 ;
 144             }
 145         }
 146         return fIntersections->used();
 147     }
 148
 149     int horizontalIntersect(double axisIntercept, double roots[3]) {
 150         double A, B, C, D;
 151         SkDCubic::Coefficients(&fCubic[0].fY, &A, &B, &C, &D);
 152         D -= axisIntercept;
 153         return SkDCubic::RootsValidT(A, B, C, D, roots);
 154     }
 155
 156     int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
 157         addExactHorizontalEndPoints(left, right, axisIntercept);
 158         if (fAllowNear) {
 159             addNearHorizontalEndPoints(left, right, axisIntercept);
 160         }
 161         double rootVals[3];
 162         int roots = horizontalIntersect(axisIntercept, rootVals);
 163         for (int index = 0; index < roots; ++index) {
 164             double cubicT = rootVals[index];
 165             SkDPoint pt = fCubic.ptAtT(cubicT);
 166             double lineT = (pt.fX - left) / (right - left);
 167             if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
 168                 fIntersections->insert(cubicT, lineT, pt);
 169             }
 170         }
 171         if (flipped) {
 172             fIntersections->flip();
 173         }
 174         return fIntersections->used();
 175     }
 176
 177     int verticalIntersect(double axisIntercept, double roots[3]) {
 178         double A, B, C, D;
 179         SkDCubic::Coefficients(&fCubic[0].fX, &A, &B, &C, &D);
 180         D -= axisIntercept;
 181         return SkDCubic::RootsValidT(A, B, C, D, roots);
 182     }
 183
 184     int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
 185         addExactVerticalEndPoints(top, bottom, axisIntercept);
 186         if (fAllowNear) {
 187             addNearVerticalEndPoints(top, bottom, axisIntercept);
 188         }
 189         double rootVals[3];
 190         int roots = verticalIntersect(axisIntercept, rootVals);
 191         for (int index = 0; index < roots; ++index) {
 192             double cubicT = rootVals[index];
 193             SkDPoint pt = fCubic.ptAtT(cubicT);
 194             double lineT = (pt.fY - top) / (bottom - top);
 195             if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
 196                 fIntersections->insert(cubicT, lineT, pt);
 197             }
 198         }
 199         if (flipped) {
 200             fIntersections->flip();
 201         }
 202         return fIntersections->used();
 203     }
 204
 205     protected:
 206
 207     void addExactEndPoints() {
 208         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
 209             double lineT = fLine.exactPoint(fCubic[cIndex]);
 210             if (lineT < 0) {
 211                 continue;
 212             }
 213             double cubicT = (double) (cIndex >> 1);
 214             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
 215         }
 216     }
 217
 218     /* Note that this does not look for endpoints of the line that are near the cubic.
 219        These points are found later when check ends looks for missing points */
 220     void addNearEndPoints() {
 221         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
 222             double cubicT = (double) (cIndex >> 1);
 223             if (fIntersections->hasT(cubicT)) {
 224                 continue;
 225             }
 226             double lineT = fLine.nearPoint(fCubic[cIndex]);
 227             if (lineT < 0) {
 228                 continue;
 229             }
 230             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
 231         }
 232     }
 233
 234     void addExactHorizontalEndPoints(double left, double right, double y) {
 235         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
 236             double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
 237             if (lineT < 0) {
 238                 continue;
 239             }
 240             double cubicT = (double) (cIndex >> 1);
 241             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
 242         }
 243     }
 244
 245     void addNearHorizontalEndPoints(double left, double right, double y) {
 246         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
 247             double cubicT = (double) (cIndex >> 1);
 248             if (fIntersections->hasT(cubicT)) {
 249                 continue;
 250             }
 251             double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
 252             if (lineT < 0) {
 253                 continue;
 254             }
 255             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
 256         }
 257         // FIXME: see if line end is nearly on cubic
 258     }
 259
 260     void addExactVerticalEndPoints(double top, double bottom, double x) {
 261         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
 262             double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
 263             if (lineT < 0) {
 264                 continue;
 265             }
 266             double cubicT = (double) (cIndex >> 1);
 267             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
 268         }
 269     }
 270
 271     void addNearVerticalEndPoints(double top, double bottom, double x) {
 272         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
 273             double cubicT = (double) (cIndex >> 1);
 274             if (fIntersections->hasT(cubicT)) {
 275                 continue;
 276             }
 277             double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
 278             if (lineT < 0) {
 279                 continue;
 280             }
 281             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
 282         }
 283         // FIXME: see if line end is nearly on cubic
 284     }
 285
 286     double findLineT(double t) {
 287         SkDPoint xy = fCubic.ptAtT(t);
 288         double dx = fLine[1].fX - fLine[0].fX;
 289         double dy = fLine[1].fY - fLine[0].fY;
 290         if (fabs(dx) > fabs(dy)) {
 291             return (xy.fX - fLine[0].fX) / dx;
 292         }
 293         return (xy.fY - fLine[0].fY) / dy;
 294     }
 295
 296     bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
 297         if (!approximately_one_or_less(*lineT)) {
 298             return false;
 299         }
 300         if (!approximately_zero_or_more(*lineT)) {
 301             return false;
 302         }
 303         double cT = *cubicT = SkPinT(*cubicT);
 304         double lT = *lineT = SkPinT(*lineT);
 305         SkDPoint lPt = fLine.ptAtT(lT);
 306         SkDPoint cPt = fCubic.ptAtT(cT);
 307         if (!lPt.moreRoughlyEqual(cPt)) {
 308             return false;
 309         }
 310         // FIXME: if points are roughly equal but not approximately equal, need to do
 311         // a binary search like quad/quad intersection to find more precise t values
 312         if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
 313             *pt = lPt;
 314         } else if (ptSet == kPointUninitialized) {
 315             *pt = cPt;
 316         }
 317         SkPoint gridPt = pt->asSkPoint();
 318         if (gridPt == fLine[0].asSkPoint()) {
 319             *lineT = 0;
 320         } else if (gridPt == fLine[1].asSkPoint()) {
 321             *lineT = 1;
 322         }
 323         if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
 324             *cubicT = 0;
 325         } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
 326             *cubicT = 1;
 327         }
 328         return true;
 329     }
 330
 331 private:
 332     const SkDCubic& fCubic;
 333     const SkDLine& fLine;
 334     SkIntersections* fIntersections;
 335     bool fAllowNear;
 336 };
 337
 338 int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
 339         bool flipped) {
 340     SkDLine line = {{{ left, y }, { right, y }}};
 341     LineCubicIntersections c(cubic, line, this);
 342     return c.horizontalIntersect(y, left, right, flipped);
 343 }
 344
 345 int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
 346         bool flipped) {
 347     SkDLine line = {{{ x, top }, { x, bottom }}};
 348     LineCubicIntersections c(cubic, line, this);
 349     return c.verticalIntersect(x, top, bottom, flipped);
 350 }
 351
 352 int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
 353     LineCubicIntersections c(cubic, line, this);
 354     c.allowNear(fAllowNear);
 355     return c.intersect();
 356 }
 357
 358 int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
 359     LineCubicIntersections c(cubic, line, this);
 360     fUsed = c.intersectRay(fT[0]);
 361     for (int index = 0; index < fUsed; ++index) {
 362         fPt[index] = cubic.ptAtT(fT[0][index]);
 363     }
 364     return fUsed;
 365 }