Formula for the intersection of straight lines

[davidgraham]

This is not a specific fortran problem but I know a lot of fortran users are mathematicians.

I'm using a cotangent formula to calculate the intersection of two lines from the coordinates of the four points.
The formula is failing when one of both of the lines is near vertical or horizontal.

Can anyone point me in the right direction?

Thanks

[rase]

Take a look at John Burkhard's collection of Fortran, C/C++ and other libraries, especially the geometry library
http://people.sc.fsu.edu/~jburkardt/f_src/geometry/geometry.html
The subroutine lines_exp_int_2d is exactly what you are looking for, among other geometry jewels.
Greetings, Wolf

[SergeyKostrov]

Quoting davidgraham

...
The formula is failing when one of both of the lines is near vertical or horizontal.
...

Consider a special processing forthese cases:

Case 1: Whena line is vertical x2 -x1 = 0 and if it is used in the denominator itcreatesa divide-by-zero
situation.
Case 2: Both lines are vertical - no intersection
Case 3: Both lines are horizontal - no intersection

[GVautier]

Hello

In order to avoid problems with vertical or horizontal line, it is better to use the canonical equation of the line :
u*x + v*y + w = 0

where u + v=1

With that equation, the intersection between two lines is computed by solving the equations system :
u1*xi + v1*yi + w1 = 0
u2*xi + v2*yi + w2 = 0

if the determinant (u1*v2-u2*v1) is "nul" lines are parallel if not you can compute(xi,yi) the coordinates of the intersection point :
xi = (w2*v1-w1*v2) / (u1*v2-u2*v1)
yi = (w1*u2-w2*u1) / (u1*v2-u2*v1)

That equation has many other interesting properties like the distance from a point (xp,yp) to the line which is u*xp+v*yp+w ....

[davidgraham]

Thanks,
That looks good, I will have to give it a try.
David

[anthonyrichards]

With that formula, you will still have a problem if you work out that the slope of one of the lines is infinity (either u or v infinite). i.e. the line is parallel to the Y-axis.

If a line is defined by (x1,y1) and (x2,y2), first check whether abs(x1-x2) is less than the smallest quantity you want to consider, epsilon. If that is the case, then the equation of the line is x=x1 ('vertical'). Otherwise it could be ax+by+c OR y=y1 ('horizontal'), the latter applying if abs(y1-y2) is less than epsilon.
Do the same for the other pair of points defining the other line. Set flags according to decisions. Draw a flow diagram detailing the options available at each decision point and encode accordingly.

The possible 9 pairs of solutions to consider are :

(1) (2) (3) (4) (5) (6)
a1x+b1y+c1=0 a1x+b1y+c1=0 a1x+b1y+c1=0 x+c1=0 x+c1=0 x+c1=0
a2x+b2y+c2=0 x+c2=0 y+c2=0 a2x+b2y+c2=0 x+c2=0 y+c2=0

(7) (8) (9)
y+c1=0 y+c1=0 y+c1=0
a2x+b2y+c2=0 x+c2=0 y+c2=0

only the first pair might benefit from the matrix soution, simple substitution and re-arrangement will work for the rest, except for solution pairs 5 (both lines vertical to within epsilon) and 9 (both lines horizontal to within epsilon). Use the decision flags to determine the particular path required.

If both lines are judged to be 'vertical' to within epsilon, then you can be sure that the intersection point will be further than (x1-x2)/(2*epsilon) away in the Y-direction, from one of the points on one of the lines, if x1 - x2 is the seperation of the vertical lines.

If both lines are judged to be 'horizontal' to within epsilon, then you can be sure that the intersection point will be further than (y1-y2)/(2*epsilon) away in the X-direction, from one of the points on one of the lines, if y1 - y2 is the seperation of the horizontal lines.

...and this is just the 'simple' 2-D case!

[mecej4]

> you will still have a problem if you work out that the slope of one of the lines is infinity (either u or v infinite).

That cannot happen -- you overlook the normalization u² + v² = 1.

[anthonyrichards]

um...how do you normalise u²+v² when one of u or v or both are infinity (or Nan)?
if an equation for a line is y=mx+c, or mx-y+c=0, with m=u and 1=v and m=tan(90), how do you normalise?

P.S. in 3-D, if a 3-vector from the origin to a general point P on a line is written as P=P1+S1*v1, where P1 is a 3-vector from the origin to any starting point on the line and v1 is a unit vector parallel to the line and S1 is a scalar length, then if P is a 3-vector from the origin to a general point on a second line P=P2+S2*v2, then you can solve for the points on the two lines which have the minimum separation using the following:

If
A1=(P1-P2)*v1 (dot product)
A2=(P1-P2)*v2 (dot product)

B=v1*v2 (dot product)

(all given quantities on the RHS), then

S2=(A1*B-A2)/(B²-1)
S1=S2*B-A1

will give the values for S1 and S2 for locating the points, one on each line, having the least separation

(beware B=1, which happens when the lines are parallel, so look out for this)
The seperation is then |P1-P2| found using the solutions S1 and S2 in the above formulae for P1 and P2, and if the distance = 0 (to within epsilon), the lines intersect.

[mecej4]

> um...how do you normalise u²+v² when one of u or v or both are infinity (or Nan)?

The normalization takes care of this; u and v are never going to exceed 1 in magnitude. A horizontal line has u = 0, a vertical line has v = 0.

A slightly different explanation:

Start with a x + b y + c = 0 as the equation of the line, with a, b, c in R¹ ; either a or b may be zero, but not both. Divide by (a² + b²). For horizontal or vertical lines, take limits after doing the division; call the resulting coefficients u, v, w.