Straight line intersection: 2 by 2 system of linear equations, what solver ?

[GianLuigi Piacentini]

Hi all,

so far, I solved this kind of problem as a 2 by 2 system of linear
(parametric) equation, using a routine from a math library (SLATEC). This
solver makes acrobacies to ensure numerical stability, minimize round-off
errors, and so on. That's fine, but seem to me an overkill for a 2 by 2.
Mantyla in his 'Introduction to Solid Modeling', solves directly (I think
it's Cramer's rule, but equations are so simple...) without too much worry
about numerical aspects. He only checks for a non-zero denominator.
So far I do not have speed or memory constraint, I can go well with SLATEC,
but I would like to hear experts opinion on the subject.
If you will point me to alternative solvers, consider that it's an hobby
project and only public domain stuff is eligible (possibly Fortran,
alternatively Ada, Basic, c or Pascal - languages that I hope to be able to
understand) .

Thanks in advance.

Gigi

[Noskosteve]

Post by GianLuigi Piacentini
Hi all,
so far, I solved this kind of problem as a 2 by 2 system of linear
(parametric) equation, using a routine from a math library (SLATEC). This
solver makes acrobacies to ensure numerical stability, minimize round-off
errors, and so on. That's fine, but seem to me an overkill for a 2 by 2.
Mantyla in his 'Introduction to Solid Modeling', solves directly (I think
it's Cramer's rule, but equations are so simple...) without too much worry
about numerical aspects. He only checks for a non-zero denominator.
So far I do not have speed or memory constraint, I can go well with SLATEC,
but I would like to hear experts opinion on the subject.
If you will point me to alternative solvers, consider that it's an hobby
project and only public domain stuff is eligible (possibly Fortran,
alternatively Ada, Basic, c or Pascal - languages that I hope to be able to
understand) .
Thanks in advance.
Gigi

I know this is quite late and I do not know if this is what you are after, but....
From the C.G.A. FAQ:

Subject 1.03: How do I find intersections of 2 2D line segments?
This problem can be extremely easy or extremely difficult depends on your applications. If all you want is the intersection point, the following should work:

You can consider the word "vector" to mean "point" or "location".

Let A,B,C,D be 2-space position vectors. This means that A is located at Ax, Ay ; B is located at Bx, By; etc. Then the directed line segments AB & CD are given by:
AB=A+r(B-A), r in [0,1]
CD=C+s(D-C), s in [0,1]

This means that A & B are two points on line AB and C & D are two points on line CD.

If AB & CD intersect, then
A+r(B-A)=C+s(D-C), or

Ax+r(Bx-Ax)=Cx+s(Dx-Cx)
Ay+r(By-Ay)=Cy+s(Dy-Cy) for some r,s in [0,1]

Solving the above for r and s yields
(Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
r = ----------------------------- (eqn 1)
(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)

(Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
s = ----------------------------- (eqn 2)
(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)

Let P be the position/location (vector) of the intersection point, then
P=A+r(B-A) or

Px=Ax+r(Bx-Ax)
Py=Ay+r(By-Ay)

Note that "s" is not needed for position of intersection.

By examining the values of r & s, you can also determine some
other limiting conditions:

If 0<=r<=1 & 0<=s<=1, intersection exists
r<0 or r>1 or s<0 or s>1 line segments do not intersect

If the denominator in eqn 1 is zero, AB & CD are parallel
If the numerator in eqn 1 is also zero, AB & CD are collinear.

If they are collinear, then the segments may be projected to the x-
or y-axis, and overlap of the projected intervals checked.

If the intersection point of the 2 lines are needed (lines in this
context mean infinite lines) regardless whether the two line
segments intersect, then

If r>1, P is located on extension of AB
If r<0, P is located on extension of BA
If s>1, P is located on extension of CD
If s<0, P is located on extension of DC

Also note that the denominators of eqn 1 & 2 are identical and thus, only need to be calculated once.

References:

[O'Rourke (C)] pp. 249-51
[Gems III] pp. 199-202 "Faster Line Segment Intersection,"

[GianLuigi Piacentini]

Post by Noskosteve

Post by GianLuigi Piacentini
Hi all,
so far, I solved this kind of problem as a 2 by 2 system of linear
(parametric) equation, using a routine from a math library (SLATEC).
This
solver makes acrobacies to ensure numerical stability, minimize
round-off errors, and so on. That's fine, but seem to me an overkill for
a 2 by 2. Mantyla in his 'Introduction to Solid Modeling', solves
directly (I think it's Cramer's rule, but equations are so simple...)
without too much worry
about numerical aspects. He only checks for a non-zero denominator.
So far I do not have speed or memory constraint, I can go well with
SLATEC, but I would like to hear experts opinion on the subject.
If you will point me to alternative solvers, consider that it's an hobby
project and only public domain stuff is eligible (possibly Fortran,
alternatively Ada, Basic, c or Pascal - languages that I hope to be able
to understand) .
Thanks in advance.
Gigi

I know this is quite late and I do not know if this is what you are after,
Subject 1.03: How do I find intersections of 2 2D line segments?
This problem can be extremely easy or extremely difficult depends on
your applications. If all you want is the intersection point, the
You can consider the word "vector" to mean "point" or "location".
Let A,B,C,D be 2-space position vectors. This means that A is located
at Ax, Ay ; B is located at Bx, By; etc. Then the directed line
AB=A+r(B-A), r in [0,1]
CD=C+s(D-C), s in [0,1]
This means that A & B are two points on line AB and C & D are two points on line CD.
If AB & CD intersect, then
A+r(B-A)=C+s(D-C), or
Ax+r(Bx-Ax)=Cx+s(Dx-Cx)
Ay+r(By-Ay)=Cy+s(Dy-Cy) for some r,s in [0,1]
Solving the above for r and s yields
(Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
r = ----------------------------- (eqn 1)
(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
(Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
s = ----------------------------- (eqn 2)
(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
Let P be the position/location (vector) of the intersection point, then
P=A+r(B-A) or
Px=Ax+r(Bx-Ax)
Py=Ay+r(By-Ay)
Note that "s" is not needed for position of intersection.
By examining the values of r & s, you can also determine some
If 0<=r<=1 & 0<=s<=1, intersection exists
r<0 or r>1 or s<0 or s>1 line segments do not intersect
If the denominator in eqn 1 is zero, AB & CD are parallel
If the numerator in eqn 1 is also zero, AB & CD are collinear.
If they are collinear, then the segments may be projected to the x-
or y-axis, and overlap of the projected intervals checked.
If the intersection point of the 2 lines are needed (lines in this
context mean infinite lines) regardless whether the two line
segments intersect, then
If r>1, P is located on extension of AB
If r<0, P is located on extension of BA
If s>1, P is located on extension of CD
If s<0, P is located on extension of DC
Also note that the denominators of eqn 1 & 2 are identical and thus,
only need to be calculated once.
[O'Rourke (C)] pp. 249-51
[Gems III] pp. 199-202 "Faster Line Segment Intersection,"

Thanks.

Seems to me that we have the usual problems: possible cancellation in the
floating point subtraction (I think we cannot overcome this) and the
equalities to zero.
We are living in a floating point world, so probably the equalities are to
be transformed into ABS(value) < epsilon, where epsilon is a small number.
And now the problem is in properly choosing epsilon.

Any guidance ?

Gigi

[Noskosteve]

Post by GianLuigi Piacentini

Post by Noskosteve

Post by GianLuigi Piacentini
Hi all,
so far, I solved this kind of problem as a 2 by 2 system of linear
(parametric) equation, using a routine from a math library (SLATEC).
This
solver makes acrobacies to ensure numerical stability, minimize
round-off errors, and so on. That's fine, but seem to me an overkill for
a 2 by 2. Mantyla in his 'Introduction to Solid Modeling', solves
directly (I think it's Cramer's rule, but equations are so simple...)
without too much worry
about numerical aspects. He only checks for a non-zero denominator.
So far I do not have speed or memory constraint, I can go well with
SLATEC, but I would like to hear experts opinion on the subject.
If you will point me to alternative solvers, consider that it's an hobby
project and only public domain stuff is eligible (possibly Fortran,
alternatively Ada, Basic, c or Pascal - languages that I hope to be able
to understand) .
Thanks in advance.
Gigi

I know this is quite late and I do not know if this is what you are after,
Subject 1.03: How do I find intersections of 2 2D line segments?
This problem can be extremely easy or extremely difficult depends on
your applications. If all you want is the intersection point, the
You can consider the word "vector" to mean "point" or "location".
Let A,B,C,D be 2-space position vectors. This means that A is located
at Ax, Ay ; B is located at Bx, By; etc. Then the directed line
AB=A+r(B-A), r in [0,1]
CD=C+s(D-C), s in [0,1]
This means that A & B are two points on line AB and C & D are two points on line CD.
If AB & CD intersect, then
A+r(B-A)=C+s(D-C), or
Ax+r(Bx-Ax)=Cx+s(Dx-Cx)
Ay+r(By-Ay)=Cy+s(Dy-Cy) for some r,s in [0,1]
Solving the above for r and s yields
(Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
r = ----------------------------- (eqn 1)
(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
(Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
s = ----------------------------- (eqn 2)
(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
Let P be the position/location (vector) of the intersection point, then
P=A+r(B-A) or
Px=Ax+r(Bx-Ax)
Py=Ay+r(By-Ay)
Note that "s" is not needed for position of intersection.
By examining the values of r & s, you can also determine some
If 0<=r<=1 & 0<=s<=1, intersection exists
r<0 or r>1 or s<0 or s>1 line segments do not intersect
If the denominator in eqn 1 is zero, AB & CD are parallel
If the numerator in eqn 1 is also zero, AB & CD are collinear.
If they are collinear, then the segments may be projected to the x-
or y-axis, and overlap of the projected intervals checked.
If the intersection point of the 2 lines are needed (lines in this
context mean infinite lines) regardless whether the two line
segments intersect, then
If r>1, P is located on extension of AB
If r<0, P is located on extension of BA
If s>1, P is located on extension of CD
If s<0, P is located on extension of DC
Also note that the denominators of eqn 1 & 2 are identical and thus,
only need to be calculated once.
[O'Rourke (C)] pp. 249-51
[Gems III] pp. 199-202 "Faster Line Segment Intersection,"

Thanks.
Seems to me that we have the usual problems: possible cancellation in the
floating point subtraction (I think we cannot overcome this) and the
equalities to zero.
We are living in a floating point world, so probably the equalities are to
be transformed into ABS(value) < epsilon, where epsilon is a small number.
And now the problem is in properly choosing epsilon.
Any guidance ?
Gigi

I can not add anything. I used this method on a floating point platform several years ago and had no problems.
Ciao, Steve