How to compute which way data points continue beyond an intersection?

Question

Let's say you have two arrays of data values from a calculation, that you can model with a continuos, differentiable function each. Both "lines" of data points intersect at (at least) one point and now the question is whether the functions behind these datasets are actually crossing or anticrossing.

The image below shows the situation, where I know (from the physics behind it) that at the upper two "contact points" the yellow and green lines actually should "switch color", whereas at the lower one both functions go out of each others way:

To give an easier "toy set" of data, take this code for example:

import matplotlib.pyplot as plt
import numpy as np

x=np.arange(-10,10,.5)
y1=[np.absolute(i**3)+100*np.absolute(i) for i in x]
y2=[-np.absolute(i**3)-100*np.absolute(i) for i in x][::-1]

plt.scatter(x,y1)
plt.scatter(x,y2,color='r')

plt.show()

Which should produce the following image:

Now how could I extrapolate whether the trend behind the data is crossing (so the data from the lower left continues to the upper right) or anti-crossing (as indicated with the colors above, the data from the lower left continues to the lower right)?

So far I was able to find the "contact point" between these to datasets by looking at the derivative of the Difference between them, roughly like this:

closePoints=np.where(np.diff(np.diff(array_A - array_B) > 0))[0] + 1 

(which probably would be faster to evaluate with something like scipy's cKDTree).

Should I go on and (probably very inefficiently) check the derivative on both sides of the intersection? Or can I somehow check if the extrapolation of the data on the left side fits better to crossing or anticrossing?

Answer 1

I understood your problem as:

• You have two sequences of points in a 2D plane.
• The true curves can be approximated by straight lines between consecutive points of the sequences.
• You want to know how often and where the two curves intersect (not only come into contact but really cross each other) (polygon intersection).

A potential solution is:

• You look at each combination of a line segment of one curve with a line segment of another curve.
• Combinations where the bounding boxes of the line segments have an overlap can potentially contain intersection points.
• You solve a linear equation system to compute if and where an intersection between two lines occurs
• In case of no solution to the equation system the lines are parallel but not overlapping, dismiss this case
• In case of one solution check that it is truly within the segments, if so record this crossing point
• In case of infinitely many intersections the lines are identical. This is also no real crossing and can be dismissed.
• Do this for all combinations of line segments and eliminate twin cases, i.e. where the two curves intersect at a segment start or end

Let me give some details:

How to check if two bounding-boxes (rectangles) of the segments overlap so that the segments potentially can intersect?

The minimal x/y value of one rectangle must be smaller than the maximal x/y value of the other. This must hold for both.

If you have two segments how do you solve for intersection points?

Let's say segment A has two points (x1, y1) and (x2, y2) and segment B has two points (x2, y3) and (x4, y4).

Then you simply have two parametrized line equations which have to be set equal:

(x1, y1) + t * (x2 - x1, y2 - y1) = (x3, y3) + q * (x4 - x3, y4 - y3)

And you need to find all solutions where t or q in [0, 1). The corresponding linear equation system may be rank deficient or not solvable at all, best is to use a general solver (I chose numpy.linalg.lstsq) that does everything in one go.

Curves sharing a common point

Surprisingly difficult are cases where one point is common in the segmentation of both curves. The difficulty lies then in the correct decision of real intersection vs. contact points. The solution is to compute the angle of both adjacent segments of both curves (gives 4 angles) around the common point and look at the order of the angles. If both curves come alternating when going around the equal point then it's an intersection, otherwise it isn't.

And a code example based on your data:

import math
import matplotlib.pyplot as plt
import numpy as np

def intersect_curves(x1, y1, x2, y2):
    """
        x1, y1 data vector for curve 1
        x2, y2 data vector for curve 2
    """

    # number of points in each curve, number of segments is one less, need at least one segment in each curve
    N1 = x1.shape[0]
    N2 = x2.shape[0]

    # get segment presentation (xi, xi+1; xi+1, xi+2; ..)
    xs1 = np.vstack((x1[:-1], x1[1:]))
    ys1 = np.vstack((y1[:-1], y1[1:]))
    xs2 = np.vstack((x2[:-1], x2[1:]))
    ys2 = np.vstack((y2[:-1], y2[1:]))

    # test if bounding-boxes of segments overlap
    mix1 = np.tile(np.amin(xs1, axis=0), (N2-1,1))
    max1 = np.tile(np.amax(xs1, axis=0), (N2-1,1))
    miy1 = np.tile(np.amin(ys1, axis=0), (N2-1,1))
    may1 = np.tile(np.amax(ys1, axis=0), (N2-1,1))
    mix2 = np.transpose(np.tile(np.amin(xs2, axis=0), (N1-1,1)))
    max2 = np.transpose(np.tile(np.amax(xs2, axis=0), (N1-1,1)))
    miy2 = np.transpose(np.tile(np.amin(ys2, axis=0), (N1-1,1)))
    may2 = np.transpose(np.tile(np.amax(ys2, axis=0), (N1-1,1)))
    idx = np.where((mix2 <= max1) & (max2 >= mix1) & (miy2 <= may1) & (may2 >= miy1)) # overlapping segment combinations

    # going through all the possible segments
    x0 = []
    y0 = []
    for (i, j) in zip(idx[0], idx[1]):
        # get segment coordinates
        xa = xs1[:, j]
        ya = ys1[:, j]
        xb = xs2[:, i]
        yb = ys2[:, i]
        # ax=b, prepare matrices a and b
        a = np.array([[xa[1] - xa[0], xb[0] - xb[1]], [ya[1] - ya[0], yb[0]- yb[1]]])
        b = np.array([xb[0] - xa[0], yb[0] - ya[0]])
        r, residuals, rank, s = np.linalg.lstsq(a, b)
        # if this is not a
        if rank == 2 and not residuals and r[0] >= 0 and r[0] < 1 and r[1] >= 0 and r[1] < 1:
            if r[0] == 0 and r[1] == 0 and i > 0 and j > 0:
                # super special case of one segment point (not the first) in common, need to differentiate between crossing or contact
                angle_a1 = math.atan2(ya[1] - ya[0], xa[1] - xa[0])
                angle_b1 = math.atan2(yb[1] - yb[0], xb[1] - xb[0])

                # get previous segment
                xa2 = xs1[:, j-1]
                ya2 = ys1[:, j-1]
                xb2 = xs2[:, i-1]
                yb2 = ys2[:, i-1]
                angle_a2 = math.atan2(ya2[0] - ya2[1], xa2[0] - xa2[1])
                angle_b2 = math.atan2(yb2[0] - yb2[1], xb2[0] - xb2[1])

                # determine in which order the 4 angle are
                if angle_a2 < angle_a1:
                    h = angle_a1
                    angle_a1 = angle_a2
                    angle_a2 = h
                if (angle_b1 > angle_a1 and angle_b1 < angle_a2 and (angle_b2 < angle_a1 or angle_b2 > angle_a2)) or\
                        ((angle_b1 < angle_a1 or angle_b1 > angle_a2) and angle_b2 > angle_a1 and angle_b2 < angle_a2):
                    # both in or both out, just a contact point
                    x0.append(xa[0])
                    y0.append(ya[0])
            else:
                x0.append(xa[0] + r[0] * (xa[1] - xa[0]))
                y0.append(ya[0] + r[0] * (ya[1] - ya[0]))

    return (x0, y0)

# create data

def data_A():
    # data from question (does not intersect)
    x1 = np.arange(-10, 10, .5)
    x2 = x1
    y1 = [np.absolute(x**3)+100*np.absolute(x) for x in x1]
    y2 = [-np.absolute(x**3)-100*np.absolute(x) for x in x2][::-1]
    return (x1, y1, x2, y2)

def data_B():
    # sine, cosine, should have some intersection points
    x1 = np.arange(-10, 10, .5)
    x2 = x1
    y1 = np.sin(x1)
    y2 = np.cos(x2)
    return (x1, y1, x2, y2)

def data_C():
    # a spiral and a diagonal line, showing the more general case
    t = np.arange(0, 10, .2)
    x1 = np.sin(t * 2) * t
    y1 = np.cos(t * 2) * t
    x2 = np.arange(-10, 10, .5)
    y2 = x2
    return (x1, y1, x2, y2)

def data_D():
    # parallel and overlapping, should give no intersection point
    x1 = np.array([0, 1])
    y1 = np.array([0, 0])
    x2 = np.array([-1, 3])
    y2 = np.array([0, 0])
    return (x1, y1, x2, y2)

def data_E():
    # crossing at a segment point, should give exactly one intersection point
    x1 = np.array([-1,0,1])
    y1 = np.array([0,0,0])
    x2 = np.array([0,0,0])
    y2 = np.array([-1,0,1])
    return (x1, y1, x2, y2)

def data_F():
    # contacting at one segment point, should give no intersection point
    x1 = np.array([-1,0,-1])
    y1 = np.array([-1,0,1])
    x2 = np.array([1,0,1])
    y2 = np.array([-1,0,1])
    return (x1, y1, x2, y2)

x1, y1, x2, y2 = data_F() # select the data you like here

# show example data
plt.plot(x1, y1, 'b-o')
plt.plot(x2, y2, 'r-o')

# call to intersection computation
x0, y0 = intersect_curves(x1, y1, x2, y2)
print('{} intersection points'.format(len(x0)))

# display intersection points in green
plt.plot(x0, y0, 'go')
plt.show() # zoom in to see that the algorithm is correct

I tested it extensively and should get most (all) border cases right (see data_A-F in code). Some examples:

Some Comments:

• The assumption about the line approximation is crucial. Most true curves might only be to some extent be approximable to lines locally. Because of this places where the two curves come close but to not intersect with a distance in the order of the distance of consecutive sampling points of your curve - you may obtain false positives or false negatives. The solution is then to either use more points or to use additonal knowledge about the true curves. Splines might give a lower error rate but also require more computations, better sampling of the curves would be preferable then.
• Self-intersection is trivially included when taking two times the same curve and let them intersect
• This solution has the additional advantage that it isn't restricted to curves of the form y=f(x) but it's applicable to arbitrary curves in 2D.

Answer 2

You could use a spline interpolation for the difference function g(x) = y1(x) - y(2). Finding the minimum of the square g(x)**2 would be a contact or crossing point. Looking at the first and second derivative you could decide if it is a contact point( g(x) has minimum, g'(x)==0, g''(x) != 0) or a crossing point (g(x) is a stationary point, g'(x)==0, g''(x)==0).

The following code searches for a minimum of g(x)**2 in constrained interval and then plot the derivatives. The use of a constrained interval is to find multiple points successively by excluding intervals in which previous points were.

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as sopt
import scipy.interpolate as sip

# test functions:
nocrossingTest = True
if nocrossingTest:
    f1 = lambda x: +np.absolute(x**3)+100*np.absolute(x)
    f2 = lambda x: -np.absolute(x**3)-100*np.absolute(x)
else:
    f1 = lambda x: +np.absolute(x**3)+100*x
    f2 = lambda x: -np.absolute(x**3)-100*x

xp = np.arange(-10,10,.5)
y1p, y2p = f1(xp), f2(xp) # test array


# Do Interpolation of y1-y2 to find crossing point:
g12 = sip.InterpolatedUnivariateSpline(xp, y1p - y2p) # Spline Interpolator of Difference
dg12 = g12.derivative() # spline derivative
ddg12 = dg12.derivative() # spline derivative

# Bounded least square fit to find minimal distance
gg = lambda x: g12(x)*g12(x)
rr = sopt.minimize_scalar(gg, bounds=[-1,1]) # search minium in Interval [-1,1]
x_c = rr['x'] # x value with minimum distance
print("Crossing point is at x = {} (Distance: {})".format(x_c, g12(x_c)))



fg = plt.figure(1)
fg.clf()
fg,ax = plt.subplots(1, 1,num=1)
ax.set_title("Function Values $y$")
ax.plot(xp, np.vstack([y1p,y2p]).T, 'x',)
xx = np.linspace(xp[0], xp[-1], 1000)
ax.plot(xx, np.vstack([f1(xx), f2(xx)]).T, '-', alpha=0.5)
ax.grid(True)
ax.legend(loc="best")
fg.canvas.draw()

fg = plt.figure(2)
fg.clf()
fg,axx = plt.subplots(3, 1,num=2)
axx[0].set_title("$g(x) = y_1(x) - y_2(x)$")
axx[1].set_title("$dg(x)/dx$")
axx[2].set_title("$d^2g(x)/dx^2$")
for ax,g in zip(axx, [g12, dg12, ddg12]):
    ax.plot(xx, g(xx))
    ax.plot(x_c, g(x_c), 'ro', alpha=.5)
    ax.grid(True)

fg.tight_layout()
plt.show()

The difference function show that the difference is not smooth: