Search code examples
pythonloopspycharmintersectionset-intersection

Detect points upstream and downstream of an intersection between two curves

I have two curves, defined by

X1=[9, 10.5, 11, 12, 12, 11, 10, 8, 7, 7]
Y1=[-5, -3.5, -2.5, -0.7, 1, 3, 4, 5, 5, 5]
X2=[5, 7, 9, 9.5, 10, 11, 12]
Y2=[-2, 4, 1, 0, -0.5, -0.7, -3]

They intersect each other enter image description here

and by a function which is written in the system code I am using, I can have the coordinates of the intersection.

loop1=Loop([9, 10.5, 11, 12, 12, 11, 10, 8, 7, 7],[-5, -3.5, -2.5, -0.7, 1, 3, 4, 5, 5, 5])
loop2=Loop([5, 7, 9, 9.5, 10, 11, 12], [-2, 4, 1, 0, -0.5, -0.7, -3])
x_int, y_int = get_intersect(loop1,loop2)
Intersection = [[],[]]
Intersection.append(x_int)
Intersection.append(y_int)

for both curves, I need to find the points which are upstream and downstream the intersection identified by (x_int, y_int).

What I tried is something like:

for x_val, y_val, x, y in zip(Intersection[0], Intersection[1], loop1[0], loop1[1]):
    if  abs(x_val - x) < 0.5 and abs(y_val - y) < 0.5:
        print(x_val, x, y_val, y)

The problem is that the result is extremely affected by the delta that I decide (0.5 in this case) and this gives me wrong results especially if I work with more decimal numbers (which is actually my case).

How can I make the loop more robust and actually find all and only the points which are upstream and downstream the intersection?

Many thanks for your help

Solution

  • TL;TR: loop over poly line segments and test if the intersection is betwwen the segment end points.

    A more robust (than "delta" in OP) approach is to find a segment of the polyline, which contains the intersection (or given point in general). This segment should IMO be part of the get_intersect function, but if you do not have access to it, you have to search the segment yourself.

    Because of roundoff errors, the given point does not exactly lie on the segment, so you still have some tol parameter, but the results should be "almost-insensitive" to its (very low) value.

    The approach uses simple geometry, namely dot product and cross product and their geometric meaning:

    • dot product of vector a and b divided by |a| is projection (length) of b onto the direction of a. Once more dividing by |a| normalizes the value to the range [0;1]
    • cross product of a and b is the area of the parallelogram having a and b as sides. Dividing it by square of length make it some dimensionless factor of distance. If a point lies exactly on the segment, the cross product is zero. But a small tolerance is needed for floating point numbers.
    X1=[9, 10.5, 11, 12, 12, 11, 10, 8, 7, 7]
Y1=[-5, -3.5, -2.5, -0.7, 1, 3, 4, 5, 5, 5]
X2=[5, 7, 9, 9.5, 10, 11, 12]
Y2=[-2, 4, 1, 0, -0.5, -0.7, -3]

x_int, y_int = 11.439024390243903, -1.7097560975609765

def splitLine(X,Y,x,y,tol=1e-12):
    """Function
    X,Y ... coordinates of line points
    x,y ... point on a polyline
    tol ... tolerance of the normalized distance from the segment
    returns ... (X_upstream,Y_upstream),(X_downstream,Y_downstream)
    """
    found = False
    for i in range(len(X)-1): # loop over segments
        # segment end points
        x1,x2 = X[i], X[i+1]
        y1,y2 = Y[i], Y[i+1]
        # segment "vector"
        dx = x2 - x1
        dy = y2 - y1
        # segment length square
        d2 = dx*dx + dy*dy
        # (int,1st end point) vector
        ix = x - x1
        iy = y - y1
        # normalized dot product
        dot = (dx*ix + dy*iy) / d2
        if dot < 0 or dot > 1: # point projection is outside segment
            continue
        # normalized cross product
        cross = (dx*iy - dy*ix) / d2
        if abs(cross) > tol: # point is perpendicularly too far away
            continue
        # here, we have found the segment containing the point!
        found = True
        break
    if not found:
        raise RuntimeError("intersection not found on segments") # or return None, according to needs
    i += 1 # the "splitting point" has one higher index than the segment
    return (X[:i],Y[:i]),(X[i:],Y[i:])

# plot
import matplotlib.pyplot as plt
plt.plot(X1,Y1,'y',linewidth=8)
plt.plot(X2,Y2,'y',linewidth=8)
plt.plot([x_int],[y_int],"r*")
(X1u,Y1u),(X1d,Y1d) = splitLine(X1,Y1,x_int,y_int)
(X2u,Y2u),(X2d,Y2d) = splitLine(X2,Y2,x_int,y_int)
plt.plot(X1u,Y1u,'g',linewidth=3)
plt.plot(X1d,Y1d,'b',linewidth=3)
plt.plot(X2u,Y2u,'g',linewidth=3)
plt.plot(X2d,Y2d,'b',linewidth=3)
plt.show()

    Result:

    enter image description here