Numerical calculation of imprecision
I want to calculate the intersections between a ray and a segment. To do this, I form a linear equation and look for the intersection. Now I am facing a numerical problem for one example. Shrinking 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), since this moment is the same on both the ray and the segment. But the actual return in my case is (443.19674427919506, 103.21654844284058)
The second number contains an error already in the sixth decimal place. I am guessing the error is because the rayAX and rayBX values are very close to each other. My question is, can I get a more accurate result when calculating the intersection?
source to share
Here's a more numerically robust way of getting an intersection (note that this is actually the intersection of two strings ... it seems like your original code didn't check if the intersection was within a segment):
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 X
and Y
be the endpoints of the segment. Let be P
the intersection point we are looking for. Then the ratio PX
to PY
is equal to the ratio of area ABX
to area ABY
. You can calculate area using cross products, which is what the code above does. Note that this procedure uses only one division, which helps to minimize numerical instability.
source to share
As far as I know, the best numerical stability is achieved by the Gauss or Gauss-Jordan method with full rotation.
You need to solve this 2x2 linear system for R
and S
.
(Brx - Arx).R - (Bsx - Asx).S = Asx - Arx (Bxy - Ary).R - (Bsx - Asx).S = Asy - Ary
The general turn prompts you to select the LHS ratio with the largest modulus. There are four possible options, so you have to implement four versions of the algorithm.
For example, if the top left coefficient dominates 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 precise arithmetic, this is indeed equivalent to Cramer's rule, but arguably better numerically.
Other cases are handled symmetrically.
source to share