0

I am trying to implement a function to find ray-/segment intersections in python following Gareth Rees' great instructions: https://stackoverflow.com/a/14318254/7235455 and https://stackoverflow.com/a/565282/7235455

Here's my function:

from math import radians, sin, cos
import numpy as np

def find_intersection(point0, theta, point1, point2):
    # convert arguments to arrays:
    p = np.array(point0, dtype=np.float) # ray origin
    q = np.array(point1, dtype=np.float) # segment point 1
    q2 = np.array(point2, dtype=np.float) # segment point 2
    r = np.array((cos(theta),sin(theta))) # theta as vector (= ray as vector)

    s = q2 - q # vector from point1 to point2
    rxs = np.cross(r,s)
    qpxs = np.cross(q-p,s)
    qpxr = np.cross(q-p,r)
    t = qpxs/rxs
    u = qpxr/rxs

    if rxs == 0 and qpxr == 0:
        t0 = np.dot(q-p,r)/np.dot(r,r)
        t1 = np.dot(t0+s,r)/np.dot(r,r)
        return "collinear"
    elif rxs == 0 and qpxr != 0:
        return "parallel"
    elif rxs != 0 and 0 <= t and 0 <= u and u <= 1: # removed t <= 1 since ray is inifinte
        intersection = p+t*r
        return "intersection is {0}".format(intersection)
    else:
        return None

The function works fine when there is an intersection. But it does not recognize parallelism or collinearity, because the conditions rxs == 0 and qpxr == 0 are not (ever?) met. Run e.g.:

p0 = (0.0,0.0)
theta = radians(45.0)
p1 = (1.0,1.0) 
p2 = (3.0,3.0)

c = find_intersection(p0,theta,p1,p2)

which returns None. Adding a print statement for rxs and qpxr before the if-block gives

rxs =  2.22044604925e-16 qpxr =  -1.11022302463e-16

My conclusion is, the function fails to catch the conditions of the first if-statement because of floating point issues. 2.22044604925e-16 and -1.11022302463e-16 are pretty small, but unfortunately not exactly 0. I understand that floats cannot have an exact representation in binary.

Is my conclusion correct or did I miss something? Are there any ideas for an implementation avoiding this issue? Thanks a lot!

Thodor
  • 23
  • 4

2 Answers2

0

Yes, your conclusion is right, problem is in numerical stability of "parallel" predicate.

You can compare result with small number (for example, eps=1.0E-9). Its magnitude might depend on coordinates range (note that cross-product gives doubled triangle area, so normalizing eps by MaxVecLen**2 looks reasonable).

More complex but more exact option - using robust geometrical predicates like these ones. Perhaps Python/NumPy libraries for computational geometry contain some implementations for such operations.

MBo
  • 77,366
  • 5
  • 53
  • 86
  • Took me a while to dig into it. Your idea to compare with small numbers and normalization seems the most practical. The downside is, it gives some false parallel/collinear, but I guess I can live with that. – Thodor Oct 28 '17 at 10:29
0

There is a simple and safe way to address this problem.

Write the implicit equation of the ray (S(X, Y) = a X + b Y + c = 0). When you plug the coordinates of the endpoints of the segment in the function S, you get two values, let S0 and S1. If they are of opposite sign, there is an intersection between the supporting line of the ray and the segment.

In this case, the position of the intersection along the segment is given by the value of the parameter, which equals

- S0 / (S1 - S0).

This expression enjoys the property of always being computable (provided there is a change of sign) and in range [0, 1]. It allows to safely compute the intersection point.

To select only those intersections that are on the desired half-line (ray), you simply compute the sign of S(Xo, Yo), at the origin of the ray.

This procedure will not detect parallel nor collinear rays, but it doesn't matter. In any case, it produces sound results.