Detect points upstream and downstream of an intersection between two curves

Question

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

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

"the points" are vertices of the polyline? "upstream" and "downstream" means "before" and "after" the intersection? How does `get_intersect` look like? Don't you have information about the intersecting segment in that function? — Jan Stránský, Sep 24 '20 at 09:57
@JanStránský the points can be any point in the polyline. Yes, upstream and downstream I mean before and after. Unfortunately I don't have info about the function and I don't have access to see it. — Dario Vaccaro, Sep 24 '20 at 10:03
Please provide the values of `x_int, y_int`. You can simple iterate over the polyline segments and check distance of the intersection from the segment. If it is below some limit (`1e-12` let's say), then you have the "border" — Jan Stránský, Sep 24 '20 at 10:21
@JanStránský well in this case **x_int=11.439024390243903** and **y_int=-1.7097560975609765**. But in my case curves are defined by much more points and each point has many decimal numbers — Dario Vaccaro, Sep 24 '20 at 10:53

Jan Stránský · Accepted Answer · 2020-09-24T13:37:35.823

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:

This can actually work! I will implement it in my actual code! Many thanks! — Dario Vaccaro, Sep 24 '20 at 15:54

Detect points upstream and downstream of an intersection between two curves

1 Answers1