public static int PointOfIntersection(Point3D La, Point3D Lb, Point3D P0, Point3D P1, Point3D P2, out Point3D intersection) { //thanks to Paul Bourke and Bryan Hanson //la,Lb constitute the line. P0,P1,P2 constitute the plane intersection = new Point3D(0, 0, 0); Vector3D N = CrossProduct(P0, P1, P2); //N is the normal of the plane double n = DotProduct(N, new Vector3D(P2.X - La.X, P2.Y - La.Y, P2.Z - La.Z)); Vector3D LbMinusLa = new Vector3D(Lb.X - La.X, Lb.Y - La.Y, Lb.Z - La.Z); double d = DotProduct(N, LbMinusLa); if (d == 0) return -1;//Line is parallel to or in the plane double u = n / d; if ((u >= 0.0) && (u <= 1.0)) {//Plane is between the two points }//can be used for checking but does not influence the outcome intersection = new Point3D(La.X + u * LbMinusLa.X, La.Y + u * LbMinusLa.Y, La.Z + u * LbMinusLa.Z); return 1; } public static Vector3D CrossProduct(Point3D p0, Point3D p1, Point3D p2) //3 points constituting the plane {//parallel planes have the same normal Vector3D A = new Vector3D(p1.X - p0.X, p1.Y - p0.Y, p1.Z - p0.Z); Vector3D B = new Vector3D(p2.X - p0.X, p2.Y - p0.Y, p2.Z - p0.Z); return CrossProduct(A, B); } public static Vector3D CrossProduct(Vector3D A, Vector3D B) //A and B are vectors ie magnitude and direction only { double X = A.Y * B.Z - A.Z * B.Y; double Y = A.Z * B.X - A.X * B.Z; double Z = A.X * B.Y - A.Y * B.X; return new Vector3D(X, Y, Z); } public static double DotProduct(Vector3D A, Vector3D B) { double answer = A.X * B.X + A.Y * B.Y + A.Z * B.Z; return answer; }