CODEWITHSUNDEEP
pythonnumpygeometry
DrRonne
DrRonne

Reputation: 103

Intersections for 3D lines

I have written a function that should calculate all intersectionpoints between a line and all lines that are given to it, and this in 3D. I have these lines parametrized because that seemed to be the easiest way to work with things. The problem is that when I input the variables "t1" and "t2" back into the functions of the lines, there seems to be an inaccuracy that is too big to be acceptable for the thing that I need.

t1 is the parameter for the line of which you would like to know all intersections, so it's written in this form:

x = xo + t1 * dx
y = yo + t1 * dy
z = zo + t1 * dz

Where [xo, yo, zo] represent a point on the line that I call the "origin" and [dx, dy, dz] represents the direction of that line. The other lines are given in the same form and the function I wrote basically solves the following equation:

xo1 + t1 * dx1 = xo2 + t2 * dx2
yo1 + t1 * dy1 = yo2 + t2 * dy2
zo1 + t1 * dz1 = zo2 + t2 * dz2

Where everything is given except for t1 and t2, that's what I'm looking for here. However, I don't think actually finding t1 and t2 is the problem, I do have a solution that gives me some kind of result. As mentioned earlier, the problem is really that when I feed t1 and t2 back into these formulas to get the actual intersectionpoints, that they differ slightly from eachother. I'm talking about differences that are mostly 0.005-0.05 away from eachother in euclidean distance. But in extreme cases it could be up to 0.5 inaccuracy. I am aware that most lines in 3D do not intersect and therefore do not have a solution to these equations, but for the tests that I'm doing right now, I am 100% sure that all of the lines are within the same plane, but some might be parallel to each other. However, these inaccuracies occur for all lines, and I'm really just looking for a solution that gets it accuratly when they do intersect.

Here's the code I have for this:

def lineIntersection(self, lines):
    origins = np.zeros((len(lines), 3), dtype=np.float32)
    directions = np.zeros((len(lines), 3), dtype=np.float32)
    for i in range(0, len(lines)):
        origins[i] = lines[i].origin
        directions[i] = lines[i].direction

    ox = origins[:, 0]
    oy = origins[:, 1]
    dx = self.origin[0]
    dy = self.origin[1]
    x1 = directions[:, 0]
    y1 = directions[:, 1]
    x2 = self.direction[0]
    y2 = self.direction[1]
    t2 = np.divide((np.subtract(np.add(oy, np.multiply(np.divide(np.subtract(dx, ox), x1), y1)), dy)), np.subtract(y2, np.multiply(np.divide(x2, x1), y1)))
    t1 = np.divide((np.add(dx, np.subtract(np.multiply(t2, x2), ox))), x1)

    testx1 = np.add(ox, np.multiply(t1, x1))
    testx2 = np.add(dx, np.multiply(t2, x2))
    testy1 = np.add(oy, np.multiply(t1, y1))
    testy2 = np.add(dy, np.multiply(t2, y2))
    testz1 = np.add(origins[:, 2], np.multiply(t1, directions[:, 2]))
    testz2 = np.add(self.origin[2], np.multiply(t2, self.direction[2]))

    arr1 = np.array([testx1, testy1, testz1]).T
    arr2 = np.array([testx2, testy2, testz2]).T
    diff = np.linalg.norm(np.subtract(arr1, arr2), axis=1)
    narr = arr1[diff < 0.05] #Filtering out points that aren't actually intersecting
    nt2 = t2[diff < 0.05]

    return narr, nt2

This function is located in the "Line" class and has an origin and direction as explained earlier. The input it takes, is an array of objects from the "Line" class.

So to be clear, I'm asking why this doesn't seem to work as precise as I want it to be and how I can fix it. Or, if there are alternatives to calculating intersectionpoints that are really accurate, I would love to hear about it.

Upvotes: 0

Views: 2905

Answers (1)

MBo
MBo

Reputation: 80187

Inaccuracy is common case for intersection of lines forming small angle.
I did not checked your algo correctness, but seems you just solve system of three equation with linalg solver.
In case of almost parallel lines intermediate values (determinant) might be small causing significant errors.
Have you tried more robust numeric algorithms like SVD?

But perhaps you really don't need them:

Note that when you are sure that all lines lie in the same plane, you can exploit 2D algorithm - just check what component of dx,dy,dz have the smallest magnitude (check for some distinct lines) and ignore corresponding component - it is similar to projecting of lines onto OXY or OXZ or OYZ plane. 2D code should be much simpler.

For true 3D case there is well-tested vector approach intended to find distance (the shortest line segment) between two skew lines - it is just zero length for intersecting ones. Example here.

Note that det (determinant) magnitude is evaluated to check for parallel (and almost parallel) lines too.

Upvotes: 0

Related Questions