Numerical inaccuracy calculating intersection

Question

I want to calculate intersections between a ray and a segment. For this I form the linear equation and look for an intersection. Now I encounter a numerical problem for one example. An abbreviation of my code:

public class Test {
    public static void main(String[] args) {
        double rayAX = 443.19661703858895d;
        double rayAY = 666.3485960845833d;

        double rayBX = 443.196744279195d;
        double rayBY = 103.21654864924565d;

        double segAX = 450.0d;
        double segAY = 114.42801992127828d;

        double segBX = 443.196744279195d;
        double segBY = 103.21654864924565d;


        double a1 = (rayBY - rayAY) / (rayBX - rayAX);
        double t1 = rayAY - rayAX * a1;

        double a2 = (segBY - segAY) / (segBX - segAX);
        double t2 = segAY - segAX * a2;

        double x = (t2 - t1) / (a1 - a2);
        double y = a1 * x + t1;

        System.out.println(x);
        System.out.println(y);
    }
}

Obviously the return should be (443.196744279195, 103.21654864924565) as this point is the same on both the ray and the segment. But the actual return is in my case (443.19674427919506, 103.21654844284058)

In the second number there is an error already in the sixth decimal place. I guess the error is because the values rayAX and rayBX are very close to each other. My question is: Can I get a more precise result when calculating the intersection?

Answer 1

Here's a more numerically stable way of getting the intersection (note that it's actually the intersection of two lines... it seems like your original code didn't check if the intersection was within the segment either):

double rX = rayBX - rayAX;                                                                                                                                           
double rY = rayBY - rayAY;                                                                                                                                           

double sAX = segAX - rayAX;                                                                                                                                          
double sAY = segAY - rayAY;                                                                                                                                          
double areaA = sAX * rY - sAY * rX;                                                                                                                                  

double sBX = segBX - rayAX;                                                                                                                                          
double sBY = segBY - rayAY;                                                                                                                                          
double areaB = sBX * rY - sBY * rX;                                                                                                                                  

double t = areaA / (areaA - areaB);
// if t is not between 0 and 1, intersection is not in segment                                                                                                 
double x = (1 - t) * segAX + t * segBX;                                                                                                                              
double y = (1 - t) * segAY + t * segBY;

Rough explanation: Let A and B be the endpoints of the ray, and let X and Y be the endpoints of the segment. Let P be the intersection point we're looking for. Then, the ratio of PX to PY is equal to the ratio of the area of ABX to the area of ABY. You can calculate the area using cross products, which is what the code above is doing. Note how this procedure only uses one division, which helps to minimize the numerical instability.