Accuracy of line intersection with vertical line

Question

I am calculating the intersection point of two lines with the following function:

// Functions of lines as per requested:
// f(y1) = starty1 + x * d1
// f(y2) = starty2 + x * d2
// x1 and y1 are the coordinates of the first point
// x2 and y2 are the coordinates of the second point
// d1 and d2 are the deltas of the corresponding lines
private static double[] intersect(double x1, double y1, double d1, double x2, double y2, double d2) {
    double starty1 = y1 - x1 * d1;
    double starty2 = y2 - x2 * d2;
    double rx = (starty2 - starty1) / (d1 - d2);
    double ry = starty1 + d1 * rx;

    tmpRes[0] = rx;
    tmpRes[1] = ry;

    return tmpRes;
}

// This is the same function, but takes 4 points to make two lines, 
// instead of two points and two deltas.
private static double[] tmpRes = new double[2];
private static double[] intersect(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4) {
    double d1 = (y1 - y2) / (x1 - x2);
    double d2 = (y3 - y4) / (x3 - x4);
    double starty1 = y1 - x1 * d1;
    double starty2 = y3 - x3 * d2;
    double rx = (starty2 - starty1) / (d1 - d2);
    double ry = starty1 + d1 * rx;

    tmpRes[0] = rx;
    tmpRes[1] = ry;

    return tmpRes;
}

However as d1 or d2 get bigger (for vertical lines) the results are getting much less accurate. How can I prevent this from happening?

For my case I have two lines perpendicular to each other. If the lines are rotated 45 degrees, I get accurate results. If the lines are at 0 or 90 degrees, I get inaccurate results (one axis of the intersection is correct, the other is all over the place.

Edit

Using a cross product:

private static double[] crTmp = new double[3];
public static double[] cross(double a, double b, double c, double a2, double b2, double c2){
    double newA = b*c2 - c*b2;
    double newB = c*a2 - a*c2;
    double newC = a*b2 - b*a2;
    crTmp[0] = newA;
    crTmp[1] = newB;
    crTmp[2] = newC;
    return crTmp;
}


public static double[] linesIntersect(double x1, double y1, double d1, double x2, double y2, double d2)
{
    double dd1 = 1.0 / d1;
    double dd2 = 1.0 / d2;

    double a1, b1, a2, b2, c1, c2;
    if (Math.abs(d1) < Math.abs(dd1)) {
        a1 = d1;
        b1 = -1.0;
        c1 = y1 - x1 * d1;
    } else {
        a1 = 1.0;
        b1 = dd1;
        c1 = -x1 - y1 * dd1;
    }
    if (Math.abs(d2) < Math.abs(dd2)) {
        a2 = d2;
        b2 = -1.0;
        c2 = y2 - x2 * d2;
    } else {
        a2 = 1.0;
        b2 = dd2;
        c2 = -x2 - y2 * dd2;
    }

    double[] v1 = {a1, b1, c1};
    double[] v2 = {a2, b2, c2};
    double[] res = cross(v1[0], v1[1], v1[2], v2[0], v2[1], v2[2]);
    tmpRes[0] = res[0] / res[2];
    tmpRes[1] = res[1] / res[2];
    return tmpRes;
}

Answer 1

I would like to make the accepted answer dead simple. Note that this is C#, so the syntax is similar to Java and should be easily understood.

// Setup formulas:

double A(double y2, double y1) => y1-y2;
double B(double x1,double x2) => x2-x1;
double C(double x1, double y1, double x2, double y2) => x1 * y2 - x2 * y1;

// Gives the Ax+By+C=0 form of the line
(double A, double B, double C) LineEquation(double x1, double y1, double x2, double y2) => 
(A(y2, y1), B(x1, x2), C(x1, y1, x2, y2));

//The Mobius version of a Cartesian point, called a "homogeneous point"
(double X, double Y, double Z) HomogeneousPointOfIntersection(double A1, double B1, double C1, double A2, double B2, double C2) =>
(B1 * C2 - B2 * C1, A2 * C1 - A1 * C2, A1 * B2 - A2 * B1);

Note the formula HomogeneousPointOfIntersection returns a triplet which is called (X,Y,Z). To obtain the actual Cartesian (X,Y) point, you have to do a final division of the X,Y,Z components.

BUT! If the Z component is zero, it means the lines are parallel and you should not do the final division as it will be a NaN error. So you want to return/escape early if Z is zero.

Here is the final sample code.

// Test that two sloping lines intersect


(double x, double y) p1 = (4, 2);
(double x, double y) p2 = (7, 5);

var line1 = LineEquation(p1.x,p1.y,p2.x,p2.y);

(double x, double y) p3 = (7, 2);
(double x, double y) p4 = (5, 7);

var line2 = LineEquation(p3.x, p3.y, p4.x, p4.y);

// h_p stands for homogeneous point
var h_p = HomogeneousPointOfIntersection(line1.A, line1.B, line1.C,line2.A,line2.B,line2.C);

if(h_p.Z == 0) return; // Return early as the below will be NaN.In this test example, this will be bypassed as we can get a successful final result.

var actualpointX = h_p.X / h_p.Z; // yields the Cartesian X value
var actualpointY = h_p.Y / h_p.Z; // yields the Cartesian Y value

Therefore, this code will handle all conditions with the one edge condition of parallel lines. Feel free to test with vertical and horizontal lines and you will see that these conditions will still return a point intersection.

Answer 2

It is easiest if you use homogeneous notation:

• Change your line representation from y = d*x + c to

  d*x - y + c = 0 = [d -1 c] . [x y 1]
  

  (where . means inner product)

• Using this notation, you can write your lines as two vectors: [d1 -1 y1] and [d2 -1 y2]

• Take the cross product of these two vectors, giving you a new vector:

  [d1 -1 y1] x [d2 -1 y2] = [a b c]
  

  (I'll leave you to look up how to calculate the cross product, but it's just simple multiplication)

The intersection of the two points is at (a/c, b/c). Provided the two lines are not parallel, c will be non-zero.

See: http://robotics.stanford.edu/~birch/projective/node4.html

One advantage of the a*x + b*y + c = 0 form of the line equation is that you can represent vertical lines naturally: you can't represent the line x = 1 by in the form y = m*x + c, because m would be infinity, whereas you can with 1*x + 0*y - 1 = 0.