Intersection work in progress
Review URL: https://codereview.appspot.com/5576043
git-svn-id: http://skia.googlecode.com/svn/trunk@3087 2bbb7eff-a529-9590-31e7-b0007b416f81
diff --git a/experimental/Intersection/LineCubicIntersection.cpp b/experimental/Intersection/LineCubicIntersection.cpp
new file mode 100644
index 0000000..f37507b
--- /dev/null
+++ b/experimental/Intersection/LineCubicIntersection.cpp
@@ -0,0 +1,139 @@
+#include "CubicIntersection.h"
+#include "CubicUtilities.h"
+#include "Intersections.h"
+#include "LineUtilities.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 )
+ */
+
+class LineCubicIntersections : public Intersections {
+public:
+
+LineCubicIntersections(const Cubic& c, const _Line& l, Intersections& i)
+ : cubic(c)
+ , line(l)
+ , intersections(i) {
+}
+
+bool intersect() {
+ double slope;
+ double axisIntercept;
+ moreHorizontal = implicitLine(line, slope, axisIntercept);
+ double A = cubic[3].x; // d
+ double B = cubic[2].x * 3; // 3*c
+ double C = cubic[1].x * 3; // 3*b
+ double D = cubic[0].x; // a
+ A -= D - C + B; // A = -a + 3*b - 3*c + d
+ B += 3 * D - 2 * C; // B = 3*a - 6*b + 3*c
+ C -= 3 * D; // C = -3*a + 3*b
+ double E = cubic[3].y; // h
+ double F = cubic[2].y * 3; // 3*g
+ double G = cubic[1].y * 3; // 3*f
+ double H = cubic[0].y; // e
+ E -= H - G + F; // E = -e + 3*f - 3*g + h
+ F += 3 * H - 2 * G; // F = 3*e - 6*f + 3*g
+ G -= 3 * H; // G = -3*e + 3*f
+ if (moreHorizontal) {
+ A = A * slope - E;
+ B = B * slope - F;
+ C = C * slope - G;
+ D = D * slope - H + axisIntercept;
+ } else {
+ A = A - E * slope;
+ B = B - F * slope;
+ C = C - G * slope;
+ D = D - H * slope - axisIntercept;
+ }
+ double t[3];
+ int roots = cubicRoots(A, B, C, D, t);
+ for (int x = 0; x < roots; ++x) {
+ intersections.add(t[x], findLineT(t[x]));
+ }
+ return roots > 0;
+}
+
+protected:
+
+double findLineT(double t) {
+ const double* cPtr;
+ const double* lPtr;
+ if (moreHorizontal) {
+ cPtr = &cubic[0].x;
+ lPtr = &line[0].x;
+ } else {
+ cPtr = &cubic[0].y;
+ lPtr = &line[0].y;
+ }
+ double s = 1 - t;
+ double cubicVal = cPtr[0] * s * s * s + 3 * cPtr[2] * s * s * t
+ + 3 * cPtr[4] * s * t * t + cPtr[6] * t * t * t;
+ return (cubicVal - lPtr[0]) / (lPtr[2] - lPtr[0]);
+}
+
+private:
+
+const Cubic& cubic;
+const _Line& line;
+Intersections& intersections;
+bool moreHorizontal;
+
+};
+
+bool intersectStart(const Cubic& cubic, const _Line& line, Intersections& i) {
+ LineCubicIntersections c(cubic, line, i);
+ return c.intersect();
+}