michael@0: /*
michael@0:  * Copyright 2012 Google Inc.
michael@0:  *
michael@0:  * Use of this source code is governed by a BSD-style license that can be
michael@0:  * found in the LICENSE file.
michael@0:  */
michael@0: #include "SkIntersections.h"
michael@0: #include "SkPathOpsCubic.h"
michael@0: #include "SkPathOpsLine.h"
michael@0: 
michael@0: /*
michael@0: Find the interection of a line and cubic by solving for valid t values.
michael@0: 
michael@0: Analogous to line-quadratic intersection, solve line-cubic intersection by
michael@0: representing the cubic as:
michael@0:   x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
michael@0:   y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
michael@0: and the line as:
michael@0:   y = i*x + j  (if the line is more horizontal)
michael@0: or:
michael@0:   x = i*y + j  (if the line is more vertical)
michael@0: 
michael@0: Then using Mathematica, solve for the values of t where the cubic intersects the
michael@0: line:
michael@0: 
michael@0:   (in) Resultant[
michael@0:         a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
michael@0:         e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
michael@0:   (out) -e     +   j     +
michael@0:        3 e t   - 3 f t   -
michael@0:        3 e t^2 + 6 f t^2 - 3 g t^2 +
michael@0:          e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
michael@0:      i ( a     -
michael@0:        3 a t + 3 b t +
michael@0:        3 a t^2 - 6 b t^2 + 3 c t^2 -
michael@0:          a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
michael@0: 
michael@0: if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
michael@0: 
michael@0:   (in) Resultant[
michael@0:         a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
michael@0:         e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y,       y]
michael@0:   (out)  a     -   j     -
michael@0:        3 a t   + 3 b t   +
michael@0:        3 a t^2 - 6 b t^2 + 3 c t^2 -
michael@0:          a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
michael@0:      i ( e     -
michael@0:        3 e t   + 3 f t   +
michael@0:        3 e t^2 - 6 f t^2 + 3 g t^2 -
michael@0:          e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
michael@0: 
michael@0: Solving this with Mathematica produces an expression with hundreds of terms;
michael@0: instead, use Numeric Solutions recipe to solve the cubic.
michael@0: 
michael@0: The near-horizontal case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
michael@0:     A =   (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d)     )
michael@0:     B = 3*(-( e - 2*f +   g    ) + i*( a - 2*b +   c    )     )
michael@0:     C = 3*(-(-e +   f          ) + i*(-a +   b          )     )
michael@0:     D =   (-( e                ) + i*( a                ) + j )
michael@0: 
michael@0: The near-vertical case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
michael@0:     A =   ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h)     )
michael@0:     B = 3*( ( a - 2*b +   c    ) - i*( e - 2*f +   g    )     )
michael@0:     C = 3*( (-a +   b          ) - i*(-e +   f          )     )
michael@0:     D =   ( ( a                ) - i*( e                ) - j )
michael@0: 
michael@0: For horizontal lines:
michael@0: (in) Resultant[
michael@0:       a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
michael@0:       e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
michael@0: (out)  e     -   j     -
michael@0:      3 e t   + 3 f t   +
michael@0:      3 e t^2 - 6 f t^2 + 3 g t^2 -
michael@0:        e t^3 + 3 f t^3 - 3 g t^3 + h t^3
michael@0:  */
michael@0: 
michael@0: class LineCubicIntersections {
michael@0: public:
michael@0:     enum PinTPoint {
michael@0:         kPointUninitialized,
michael@0:         kPointInitialized
michael@0:     };
michael@0: 
michael@0:     LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections* i)
michael@0:         : fCubic(c)
michael@0:         , fLine(l)
michael@0:         , fIntersections(i)
michael@0:         , fAllowNear(true) {
michael@0:         i->setMax(3);
michael@0:     }
michael@0: 
michael@0:     void allowNear(bool allow) {
michael@0:         fAllowNear = allow;
michael@0:     }
michael@0: 
michael@0:     // see parallel routine in line quadratic intersections
michael@0:     int intersectRay(double roots[3]) {
michael@0:         double adj = fLine[1].fX - fLine[0].fX;
michael@0:         double opp = fLine[1].fY - fLine[0].fY;
michael@0:         SkDCubic r;
michael@0:         for (int n = 0; n < 4; ++n) {
michael@0:             r[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp;
michael@0:         }
michael@0:         double A, B, C, D;
michael@0:         SkDCubic::Coefficients(&r[0].fX, &A, &B, &C, &D);
michael@0:         return SkDCubic::RootsValidT(A, B, C, D, roots);
michael@0:     }
michael@0: 
michael@0:     int intersect() {
michael@0:         addExactEndPoints();
michael@0:         if (fAllowNear) {
michael@0:             addNearEndPoints();
michael@0:         }
michael@0:         double rootVals[3];
michael@0:         int roots = intersectRay(rootVals);
michael@0:         for (int index = 0; index < roots; ++index) {
michael@0:             double cubicT = rootVals[index];
michael@0:             double lineT = findLineT(cubicT);
michael@0:             SkDPoint pt;
michael@0:             if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized)) {
michael@0:     #if ONE_OFF_DEBUG
michael@0:                 SkDPoint cPt = fCubic.ptAtT(cubicT);
michael@0:                 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
michael@0:                         cPt.fX, cPt.fY);
michael@0:     #endif
michael@0:                 for (int inner = 0; inner < fIntersections->used(); ++inner) {
michael@0:                     if (fIntersections->pt(inner) != pt) {
michael@0:                         continue;
michael@0:                     }
michael@0:                     double existingCubicT = (*fIntersections)[0][inner];
michael@0:                     if (cubicT == existingCubicT) {
michael@0:                         goto skipInsert;
michael@0:                     }
michael@0:                     // check if midway on cubic is also same point. If so, discard this
michael@0:                     double cubicMidT = (existingCubicT + cubicT) / 2;
michael@0:                     SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
michael@0:                     if (cubicMidPt.approximatelyEqual(pt)) {
michael@0:                         goto skipInsert;
michael@0:                     }
michael@0:                 }
michael@0:                 fIntersections->insert(cubicT, lineT, pt);
michael@0:         skipInsert:
michael@0:                 ;
michael@0:             }
michael@0:         }
michael@0:         return fIntersections->used();
michael@0:     }
michael@0: 
michael@0:     int horizontalIntersect(double axisIntercept, double roots[3]) {
michael@0:         double A, B, C, D;
michael@0:         SkDCubic::Coefficients(&fCubic[0].fY, &A, &B, &C, &D);
michael@0:         D -= axisIntercept;
michael@0:         return SkDCubic::RootsValidT(A, B, C, D, roots);
michael@0:     }
michael@0: 
michael@0:     int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
michael@0:         addExactHorizontalEndPoints(left, right, axisIntercept);
michael@0:         if (fAllowNear) {
michael@0:             addNearHorizontalEndPoints(left, right, axisIntercept);
michael@0:         }
michael@0:         double rootVals[3];
michael@0:         int roots = horizontalIntersect(axisIntercept, rootVals);
michael@0:         for (int index = 0; index < roots; ++index) {
michael@0:             double cubicT = rootVals[index];
michael@0:             SkDPoint pt = fCubic.ptAtT(cubicT);
michael@0:             double lineT = (pt.fX - left) / (right - left);
michael@0:             if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
michael@0:                 fIntersections->insert(cubicT, lineT, pt);
michael@0:             }
michael@0:         }
michael@0:         if (flipped) {
michael@0:             fIntersections->flip();
michael@0:         }
michael@0:         return fIntersections->used();
michael@0:     }
michael@0: 
michael@0:     int verticalIntersect(double axisIntercept, double roots[3]) {
michael@0:         double A, B, C, D;
michael@0:         SkDCubic::Coefficients(&fCubic[0].fX, &A, &B, &C, &D);
michael@0:         D -= axisIntercept;
michael@0:         return SkDCubic::RootsValidT(A, B, C, D, roots);
michael@0:     }
michael@0: 
michael@0:     int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
michael@0:         addExactVerticalEndPoints(top, bottom, axisIntercept);
michael@0:         if (fAllowNear) {
michael@0:             addNearVerticalEndPoints(top, bottom, axisIntercept);
michael@0:         }
michael@0:         double rootVals[3];
michael@0:         int roots = verticalIntersect(axisIntercept, rootVals);
michael@0:         for (int index = 0; index < roots; ++index) {
michael@0:             double cubicT = rootVals[index];
michael@0:             SkDPoint pt = fCubic.ptAtT(cubicT);
michael@0:             double lineT = (pt.fY - top) / (bottom - top);
michael@0:             if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
michael@0:                 fIntersections->insert(cubicT, lineT, pt);
michael@0:             }
michael@0:         }
michael@0:         if (flipped) {
michael@0:             fIntersections->flip();
michael@0:         }
michael@0:         return fIntersections->used();
michael@0:     }
michael@0: 
michael@0:     protected:
michael@0: 
michael@0:     void addExactEndPoints() {
michael@0:         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
michael@0:             double lineT = fLine.exactPoint(fCubic[cIndex]);
michael@0:             if (lineT < 0) {
michael@0:                 continue;
michael@0:             }
michael@0:             double cubicT = (double) (cIndex >> 1);
michael@0:             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
michael@0:         }
michael@0:     }
michael@0: 
michael@0:     /* Note that this does not look for endpoints of the line that are near the cubic.
michael@0:        These points are found later when check ends looks for missing points */
michael@0:     void addNearEndPoints() {
michael@0:         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
michael@0:             double cubicT = (double) (cIndex >> 1);
michael@0:             if (fIntersections->hasT(cubicT)) {
michael@0:                 continue;
michael@0:             }
michael@0:             double lineT = fLine.nearPoint(fCubic[cIndex]);
michael@0:             if (lineT < 0) {
michael@0:                 continue;
michael@0:             }
michael@0:             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
michael@0:         }
michael@0:     }
michael@0: 
michael@0:     void addExactHorizontalEndPoints(double left, double right, double y) {
michael@0:         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
michael@0:             double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
michael@0:             if (lineT < 0) {
michael@0:                 continue;
michael@0:             }
michael@0:             double cubicT = (double) (cIndex >> 1);
michael@0:             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
michael@0:         }
michael@0:     }
michael@0: 
michael@0:     void addNearHorizontalEndPoints(double left, double right, double y) {
michael@0:         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
michael@0:             double cubicT = (double) (cIndex >> 1);
michael@0:             if (fIntersections->hasT(cubicT)) {
michael@0:                 continue;
michael@0:             }
michael@0:             double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
michael@0:             if (lineT < 0) {
michael@0:                 continue;
michael@0:             }
michael@0:             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
michael@0:         }
michael@0:         // FIXME: see if line end is nearly on cubic
michael@0:     }
michael@0: 
michael@0:     void addExactVerticalEndPoints(double top, double bottom, double x) {
michael@0:         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
michael@0:             double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
michael@0:             if (lineT < 0) {
michael@0:                 continue;
michael@0:             }
michael@0:             double cubicT = (double) (cIndex >> 1);
michael@0:             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
michael@0:         }
michael@0:     }
michael@0: 
michael@0:     void addNearVerticalEndPoints(double top, double bottom, double x) {
michael@0:         for (int cIndex = 0; cIndex < 4; cIndex += 3) {
michael@0:             double cubicT = (double) (cIndex >> 1);
michael@0:             if (fIntersections->hasT(cubicT)) {
michael@0:                 continue;
michael@0:             }
michael@0:             double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
michael@0:             if (lineT < 0) {
michael@0:                 continue;
michael@0:             }
michael@0:             fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
michael@0:         }
michael@0:         // FIXME: see if line end is nearly on cubic
michael@0:     }
michael@0: 
michael@0:     double findLineT(double t) {
michael@0:         SkDPoint xy = fCubic.ptAtT(t);
michael@0:         double dx = fLine[1].fX - fLine[0].fX;
michael@0:         double dy = fLine[1].fY - fLine[0].fY;
michael@0:         if (fabs(dx) > fabs(dy)) {
michael@0:             return (xy.fX - fLine[0].fX) / dx;
michael@0:         }
michael@0:         return (xy.fY - fLine[0].fY) / dy;
michael@0:     }
michael@0: 
michael@0:     bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
michael@0:         if (!approximately_one_or_less(*lineT)) {
michael@0:             return false;
michael@0:         }
michael@0:         if (!approximately_zero_or_more(*lineT)) {
michael@0:             return false;
michael@0:         }
michael@0:         double cT = *cubicT = SkPinT(*cubicT);
michael@0:         double lT = *lineT = SkPinT(*lineT);
michael@0:         if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
michael@0:             *pt = fLine.ptAtT(lT);
michael@0:         } else if (ptSet == kPointUninitialized) {
michael@0:             *pt = fCubic.ptAtT(cT);
michael@0:         }
michael@0:         SkPoint gridPt = pt->asSkPoint();
michael@0:         if (gridPt == fLine[0].asSkPoint()) {
michael@0:             *lineT = 0;
michael@0:         } else if (gridPt == fLine[1].asSkPoint()) {
michael@0:             *lineT = 1;
michael@0:         }
michael@0:         if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
michael@0:             *cubicT = 0;
michael@0:         } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
michael@0:             *cubicT = 1;
michael@0:         }
michael@0:         return true;
michael@0:     }
michael@0: 
michael@0: private:
michael@0:     const SkDCubic& fCubic;
michael@0:     const SkDLine& fLine;
michael@0:     SkIntersections* fIntersections;
michael@0:     bool fAllowNear;
michael@0: };
michael@0: 
michael@0: int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
michael@0:         bool flipped) {
michael@0:     SkDLine line = {{{ left, y }, { right, y }}};
michael@0:     LineCubicIntersections c(cubic, line, this);
michael@0:     return c.horizontalIntersect(y, left, right, flipped);
michael@0: }
michael@0: 
michael@0: int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
michael@0:         bool flipped) {
michael@0:     SkDLine line = {{{ x, top }, { x, bottom }}};
michael@0:     LineCubicIntersections c(cubic, line, this);
michael@0:     return c.verticalIntersect(x, top, bottom, flipped);
michael@0: }
michael@0: 
michael@0: int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
michael@0:     LineCubicIntersections c(cubic, line, this);
michael@0:     c.allowNear(fAllowNear);
michael@0:     return c.intersect();
michael@0: }
michael@0: 
michael@0: int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
michael@0:     LineCubicIntersections c(cubic, line, this);
michael@0:     fUsed = c.intersectRay(fT[0]);
michael@0:     for (int index = 0; index < fUsed; ++index) {
michael@0:         fPt[index] = cubic.ptAtT(fT[0][index]);
michael@0:     }
michael@0:     return fUsed;
michael@0: }