Python — matplotlib: find intersection of lineplots
I have a probably simple question, that keeps me going already for quiet a while. Is there a simple way to return the intersection of two plotted (non-analytical) datasets in python matplotlib ? For elaboration, I have something like this:
x=[1.4,2.1,3,5.9,8,9,23] y=[2.3,3.1,1,3.9,8,9,11] x1=[1,2,3,4,6,8,9] y1=[4,12,7,1,6.3,8.5,12] plot(x1,y1,'k-',x,y,'b-') The data in this example is totaly arbitrary. I would now like to know if there is a simple build in function that I keep missing, that returns me the precise intersections between the two plots. Hope I made myself clear, and also that I didnt miss something totaly obvious.
3 Answers 3
We could use scipy.interpolate.PiecewisePolynomial to create functions which are defined by your piecewise-linear data.
p1=interpolate.PiecewisePolynomial(x1,y1[:,np.newaxis]) p2=interpolate.PiecewisePolynomial(x2,y2[:,np.newaxis]) We could then take the difference of these two functions,
and use optimize.fsolve to find the roots of pdiff :
import scipy.interpolate as interpolate import scipy.optimize as optimize import numpy as np x1=np.array([1.4,2.1,3,5.9,8,9,23]) y1=np.array([2.3,3.1,1,3.9,8,9,11]) x2=np.array([1,2,3,4,6,8,9]) y2=np.array([4,12,7,1,6.3,8.5,12]) p1=interpolate.PiecewisePolynomial(x1,y1[:,np.newaxis]) p2=interpolate.PiecewisePolynomial(x2,y2[:,np.newaxis]) def pdiff(x): return p1(x)-p2(x) xs=np.r_[x1,x2] xs.sort() x_min=xs.min() x_max=xs.max() x_mid=xs[:-1]+np.diff(xs)/2 roots=set() for val in x_mid: root,infodict,ier,mesg = optimize.fsolve(pdiff,val,full_output=True) # ier==1 indicates a root has been found if ier==1 and x_min [[ 3.85714286 1.85714286 1.85714286] [ 4.60606061 2.60606061 2.60606061]] The first column is the x-value, the second column is the y-value of the first PiecewisePolynomial evaluated at x , and the third column is the y-value for the second PiecewisePolynomial.
thank you very much for your time! although not as simple as I hoped for, this will certainly do in resolving my problems 🙂
@unutbu Please can you check out this question stackoverflow.com/questions/45200428/… and suggest me something
Parametric solution
If the sequences and define arbitrary (x,y) curves, rather than y(x) curves, we need a parametric approach to finding intersections. Since it's not entirely obvious how to do so, and because @unutbu's solution uses a defunct interpolator in SciPy, I thought it might be useful to revisit this question.
import numpy as np from numpy.linalg import norm from scipy.optimize import fsolve from scipy.interpolate import interp1d import matplotlib.pyplot as plt x1_array = np.array([1,2,3,4,6,8,9]) y1_array = np.array([4,12,7,1,6.3,8.5,12]) x2_array = np.array([1.4,2.1,3,5.9,8,9,23]) y2_array = np.array([2.3,3.1,1,3.9,8,9,11]) s1_array = np.linspace(0,1,num=len(x1_array)) s2_array = np.linspace(0,1,num=len(x2_array)) # Arguments given to interp1d: # - extrapolate: to make sure we don't get a fatal value error when fsolve searches # beyond the bounds of [0,1] # - copy: use refs to the arrays # - assume_sorted: because s_array ('x') increases monotonically across [0,1] kwargs_ = dict(fill_value='extrapolate', copy=False, assume_sorted=True) x1_interp = interp1d(s1_array,x1_array, **kwargs_) y1_interp = interp1d(s1_array,y1_array, **kwargs_) x2_interp = interp1d(s2_array,x2_array, **kwargs_) y2_interp = interp1d(s2_array,y2_array, **kwargs_) xydiff_lambda = lambda s12: (np.abs(x1_interp(s12[0])-x2_interp(s12[1])), np.abs(y1_interp(s12[0])-y2_interp(s12[1]))) s12_intercept, _, ier, mesg \ = fsolve(xydiff_lambda, [0.5, 0.3], full_output=True) xy1_intercept = x1_interp(s12_intercept[0]),y1_interp(s12_intercept[0]) xy2_intercept = x2_interp(s12_intercept[1]),y2_interp(s12_intercept[1]) plt.plot(x1_interp(s1_array),y1_interp(s1_array),'b.', ls='-', label='x1 data') plt.plot(x2_interp(s2_array),y2_interp(s2_array),'r.', ls='-', label='x2 data') if s12_intercept[0]>0 and s12_intercept[0], s2=<>\n'.format(s12_intercept[0],s12_intercept[1]), 'intercept @ xy1=<>\n'.format(np.array(xy1_intercept)), 'intercept @ xy2=<>\n'.format(np.array(xy2_intercept)), 'fsolve apparent success? <>: "<>"\n'.format(ier==1,mesg,), 'is intercept really good? <>\n'.format(s12_intercept[0]>=0 and s12_intercept[0]=0 and s12_intercept[1] which returns, for this particular choice of initial guess [0.5,0.3]:
intercept @ s1=0.4761904761904762, s2=0.3825944170771757 intercept @ xy1=[3.85714286 1.85714286] intercept @ xy2=[3.85714286 1.85714286] fsolve apparent success? True: "The solution converged." is intercept really good? True This method only finds one intersection: we would need to iterate over several initial guesses (as @unutbu's code does), check their veracity, and eliminate duplicates using np.close . Note that fsolve may falsely indicate successful detection of an intersection in the return value ier , which is why the extra checking is done here.
Here's the plot for this solution:
Python - Find all intersection points of 2 graphs
I'm trying to find all the intersection points of two graphs and display them on the final plot. I've looked around and tried multiple things, but I haven't been able to obtain what l'm looking for. Currently, I attempting to generate a list wherein the intersection points would be listed, though I keep getting the following error:
import numpy as np from scipy.optimize import fsolve import matplotlib.pyplot as plt x = np.arange(-7.0, 7.0, 0.05) def y(x): return np.sin(x)*(0.003*x**4 - 0.1*x**3 + x**2 + 4*x + 3) def g(x): return -10 * np.arctan(x) def intersection(x): if (y(x) - g(x)) == 0: print y.all(x) plt.plot(x, y(x), '-') plt.plot(x, g(x), '-') plt.show() 2 Answers 2
import numpy as np from scipy.optimize import fsolve import matplotlib.pyplot as plt x = np.arange(-7.0, 7.0, 0.05) y = np.sin(x)*(0.003*x**4 - 0.1*x**3 + x**2 + 4*x + 3) g = -10 * np.arctan(x) def intersection(): idx = np.argwhere(np.isclose(y, g, atol=10)).reshape(-1) print idx plt.plot(x, y, '-') plt.plot(x, g, '-') plt.show() intersection() edit: you don't use a function, but a list of values
from scipy.optimize import fsolve def findIntersection(fun1,fun2,x0): return fsolve(lambda x : fun1(x) - fun2(x),x0) result = findIntersection(y,g,0.0) Now, you just need to iterate through your range to get all the roots. This gives some duplicates, which you might be able to remove by using mpmath, setting the precision low enough, and using a set.
from scipy.optimize import fsolve import numpy as np rng = np.arange(-7.0, 7.0, 0.05) def y(x): return np.sin(x)*(0.003*x**4 - 0.1*x**3 + x**2 + 4*x + 3) def g(x): return -10 * np.arctan(x) def findIntersection(fun1,fun2,x0): return fsolve(lambda x : fun1(x) - fun2(x),x0) result = [] for x in rng: result.append(float(findIntersection(y,g,x))) Intersection of 2D polygons
I have two numpy arrays that are OpenCV convex hulls and I want to check for intersection without creating for loops or creating images and performing numpy.bitwise_and on them, both of which are quite slow in Python. The arrays look like this:
Considering [[x1 y1]] as one single element, I want to perform intersection between two numpy ndarrays. How can I do that? I have found a few questions of similar nature, but I could not figure out the solution to this from there.
4 Answers 4
You can use a view of the array as a single dimension to the intersect1d function like this:
def multidim_intersect(arr1, arr2): arr1_view = arr1.view([('',arr1.dtype)]*arr1.shape[1]) arr2_view = arr2.view([('',arr2.dtype)]*arr2.shape[1]) intersected = numpy.intersect1d(arr1_view, arr2_view) return intersected.view(arr1.dtype).reshape(-1, arr1.shape[1]) This creates a view of each array, changing each row to a tuple of values. It then performs the intersection, and changes the result back to the original format. Here's an example of using it:
test_arr1 = numpy.array([[0, 2], [1, 3], [4, 5], [0, 2]]) test_arr2 = numpy.array([[1, 2], [0, 2], [3, 1], [1, 3]]) print multidim_intersect(test_arr1, test_arr2) Thank you very much for the response! This would be perfect if all the points in the perimeter were in the numpy arrays. However, in convex hulls, I think only a few points are passed as guides. But intersection, as it means in this case, is common values in both areas, which may not be common inside the numpy arrays themselves. I just read my post above, and realized that I did not mention this at all. I am sorry for that.
When I apply your view on my numpy, it looks like this: [[[(x1,)(y1,)]][[(x2,)(y2,)]]. [[(xn,)(yn,)]]] Whereas what we really need is this: [(x1,y1), (x2,y2), (x3,y3), . (xn,yn)] Any ideas?
Do you just have an extra axis for some reason? Can you just reshape it first with test_arr1.reshape(len(test_arr1), 2) ? That will avoid a copy.
Actually, the array is generated by opencv function, and it does not even require me to set a dtype or dimension, so I have no control over how it is shaped. I will try your suggestion now.
that's testing for points being in common. this is about the intersection of the entire shapes described by the polygons/contours -- this question shouldn't be the top answer because it's wrong. how did 18 people think this was a correct answer, except if they didn't read the question well enough.
>>> import Polygon >>> a = Polygon.Polygon([(0,0),(1,0),(0,1)]) >>> b = Polygon.Polygon([(0.3,0.3), (0.3, 0.6), (0.6, 0.3)]) >>> a & b Polygon:
To convert the result of cv2.findContours to Polygon point format, you can:
points1 = contours[0].reshape(-1,2) This will convert the shape from (N, 1, 2) to (N, 2)
Following is a full example:
import Polygon import cv2 import numpy as np from scipy.misc import bytescale y, x = np.ogrid[-2:2:100j, -2:2:100j] f1 = bytescale(np.exp(-x**2 - y**2), low=0, high=255) f2 = bytescale(np.exp(-(x+1)**2 - y**2), low=0, high=255) c1, hierarchy = cv2.findContours((f1>120).astype(np.uint8), cv2.cv.CV_RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE) c2, hierarchy = cv2.findContours((f2>120).astype(np.uint8), cv2.cv.CV_RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE) points1 = c1[0].reshape(-1,2) # convert shape (n, 1, 2) to (n, 2) points2 = c2[0].reshape(-1,2) import pylab as pl poly1 = pl.Polygon(points1, color="blue", alpha=0.5) poly2 = pl.Polygon(points2, color="red", alpha=0.5) pl.figure(figsize=(8,3)) ax = pl.subplot(121) ax.add_artist(poly1) ax.add_artist(poly2) pl.xlim(0, 100) pl.ylim(0, 100) a = Polygon.Polygon(points1) b = Polygon.Polygon(points2) intersect = a&b # calculate the intersect polygon poly3 = pl.Polygon(intersect[0], color="green") # intersect[0] are the points of the polygon ax = pl.subplot(122) ax.add_artist(poly3) pl.xlim(0, 100) pl.ylim(0, 100) pl.show() How to find the intersection points between two plotted curves in matplotlib?
I want to find the intersection of two curves. For example below. There can be multiple intersection points. Right now I am finding the intersection points by finding the distance between x,y coordinates. But this method not giving accurate points sometimes when point of intersection is in between(17-18 x-axis) as you can see in the plot. I need to get all the points from the curve to solve this. Is there any method to get all of them?
2 Answers 2
Here is my approach. I first created two test curves, using only 12 sample points, just to illustrate the concept. The exact equations of the curves are lost after creating the arrays with sample points.
Then, the intersections between the two curves are searched. By going through the array point by point, and check when one curve goes from below the other to above (or the reverse), the intersection point can be calculated by solving a linear equation.
Afterwards the intersection points are plotted to visually check the result.
import numpy as np from matplotlib import pyplot as plt N = 12 t = np.linspace(0, 50, N) curve1 = np.sin(t*.08+1.4)*np.random.uniform(0.5, 0.9) + 1 curve2 = -np.cos(t*.07+.1)*np.random.uniform(0.7, 1.0) + 1 # note that from now on, we don't have the exact formula of the curves, as we didn't save the random numbers # we only have the points correspondent to the given t values fig, ax = plt.subplots() ax.plot(t, curve1,'b-') ax.plot(t, curve1,'bo') ax.plot(t, curve2,'r-') ax.plot(t, curve2,'ro') intersections = [] prev_dif = 0 t0, prev_c1, prev_c2 = None, None, None for t1, c1, c2 in zip(t, curve1, curve2): new_dif = c2 - c1 if np.abs(new_dif) < 1e-12: # found an exact zero, this is very unprobable intersections.append((t1, c1)) elif new_dif * prev_dif < 0: # the function changed signs between this point and the previous # do a linear interpolation to find the t between t0 and t1 where the curves would be equal # this is the intersection between the line [(t0, prev_c1), (t1, c1)] and the line [(t0, prev_c2), (t1, c2)] # because of the sign change, we know that there is an intersection between t0 and t1 denom = prev_dif - new_dif intersections.append(((-new_dif*t0 + prev_dif*t1) / denom, (c1*prev_c2 - c2*prev_c1) / denom)) t0, prev_c1, prev_c2, prev_dif = t1, c1, c2, new_dif print(intersections) ax.plot(*zip(*intersections), 'go', alpha=0.7, ms=10) plt.show()