CODEWITHSUNDEEP
javageometryintersectionnumeric
Geminus
Geminus

Reputation: 156

Numerical inaccuracy calculating intersection

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?

Upvotes: 3

Views: 358

Answers (2)

user1196549
user1196549

Reputation:

As far as I know, best numerical stability is achieved by the Gauss or Gauss-Jordan method with total pivoting.

You need to solve this linear 2x2 system for R and S.

(Brx - Arx).R - (Bsx - Asx).S = Asx - Arx
(Bxy - Ary).R - (Bsx - Asx).S = Asy - Ary

Total pivoting tells you to choose the LHS coefficient with the largest module. There are four possible choices, hence you will have to implement four versions of the algorithm.

For instance, assuming the top left coefficient is dominant in the system

A.X + B.Y = C
D.X + E.Y = F

Then 

  X + (B/A).Y = (C/A)
D.X + E    .Y = F

,

(E - D.(B/A)) Y = F - D.(C/A)

and

Y = (F - D.(C/A)) / (E - D.(B/A))
X = (C/A) - (B/A).Y

Using exact arithmetic this is indeed equivalent to Cramer's rule, but possibly better from a numerical standpoint.

The other cases are treated symmetrically.

Upvotes: 0

arghbleargh
arghbleargh

Reputation: 3170

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.

Upvotes: 3

Related Questions