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