如何找到和绘制等高线形状的交点

Question

我试图找到通过点 V 的线与圆锥曲线之间的交点。二次曲线相对于 y（或 x）是不可解的，所以它是用等高线描绘的。有没有一种方法可以找到等高线图的交点？ enter image description here

代码如下：

import numpy as np
import matplotlib.pyplot as plt

p = (input("choose point P on axis OX: "))
print(p)
q = (input("choose point Q on axis OX: "))
print(q)
v = (input("choose point V on axis OX: "))
print(v)

k=3
X = np.arange(-50, 50, 0.05)
Y = k*X
plt.plot(X,Y)

plt.plot(0,0)
plt.scatter(0.0, 0.0, color='white', marker='o')
plt.text(0.0, 0.0, "O", horizontalalignment="center")

plt.plot(-v,0)
plt.scatter(-v, 0, color='red', marker='o')
plt.text(-v, 0.8, "V", horizontalalignment="center")

xmin= -10
xmax= 10
ymin= -10
ymax= 10
ax = plt.gca()
ax.get_xlim()
ax.set_xlim([xmin,xmax])
ax.set_ylim([ymin,ymax])
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

#Create random point B1
b1=4
plt.plot(0.0,b1)
plt.scatter(0.0, b1, color='blue', marker='o')
plt.text(0.8, b1, "B1", horizontalalignment="center")

x, y = np.meshgrid(X, X)

#Create VB1
l3 = b1*x+b1*v - v*y
vb = plt.contour(x,y, l3, [0], colors='k')
# l3 = b1*X/v + b1
# plt.plot(X,l3)

#Create conic
conic = x*x*b1*2*p*k-x*x*b1*2*q*k+x*x*k*k+y*y-b1*2*y+2*b1*q*x*y
cnc = plt.contour(x, y, (conic), [0], colors='k')

我试过这样做：

  c = cnc.intersection(vb)
  print(c)

或

# 
    idx = np.argwhere(np.diff(np.sign(cnc - vb))).flatten()
    plt.plot(x[idx], y[idx], 'ro')

我最后一次尝试：

import numpy as np
import matplotlib.pyplot as plt
p,q,v,k,b=5,7,2,3,4
X = np.arange(-50, 50, 0.05)

plt.plot(-v,0)
plt.scatter(-v, 0, color='red', marker='o')
plt.text(-v, 0.8, "V", horizontalalignment="center")
xmin,xmax,ymin,ymax=-10,10,-10,10
ax = plt.gca()
ax.get_xlim()
ax.set_xlim([xmin,xmax])
ax.set_ylim([ymin,ymax])
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.plot(0.0,b)
plt.scatter(0.0, 1, color='blue', marker='o')
x, y = np.meshgrid(X, X)

l = b*x+b*v-v*y
vb = plt.contour(x,y, l, [0], colors='k')

conic = x*x*b*2*p*k-x*x*b*2*q*k+x*x*k*k+y*y-b*2*y+2*b*q*x*y
cnc = plt.contour(x, y, (conic), [0], colors='k')

c = cnc.collections[0].get_paths()[1]
v = c.vertices
x1 = v[:,0]
y1 = v[:,1]
plt.plot(x1,y1)

vb1 = vb.collections[0].get_paths()[0]
v1 = vb1.vertices
x2 = v1[:,0]
y2 = v1[:,1]
plt.plot(x2,y2,color='red')

# def find_roots(x,y):
#      s = np.abs(np.diff(np.sign(y))).astype(bool)
#      return x[:-1][s] + np.diff(x)[s]/(np.abs(y[1:][s]/y[:-1][s])+1)
#
# z = find_roots(x1-x2,y1-y2)
# plt.plot(z, np.zeros(len(z)), marker="o", ls="", ms=4)
plt.show()

enter image description here

Answer 1

有点复杂。问题是 (a) 轮廓中的点不一定是排序的，并且 (b) 两个轮廓没有共同的支持。

因此需要 (a) 沿 x 对点进行排序，以及 (b) 创建一个公共值数组，首先对其进行插值。

import numpy as np
import matplotlib.pyplot as plt

p,q,v,k,b=5,7,2,3,4
X = np.arange(-50, 50, 0.05)

plt.plot(-v,0)
plt.scatter(-v, 0, color='red', marker='o')
plt.text(-v, 0.8, "V", horizontalalignment="center")
xmin,xmax,ymin,ymax=-10,10,-10,10
ax = plt.gca()
ax.get_xlim()
ax.set_xlim([xmin,xmax])
ax.set_ylim([ymin,ymax])
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)


x, y = np.meshgrid(X, X)

l = b*x+b*v-v*y
vb = plt.contour(x,y, l, [0], colors='red')

conic = x*x*b*2*p*k-x*x*b*2*q*k+x*x*k*k+y*y-b*2*y+2*b*q*x*y
cnc = plt.contour(x, y, (conic), [0], colors='blue')

c = cnc.collections[0].get_paths()[1]
v = c.vertices
x1 = np.sort(v[:,0])
y1 = v[np.argsort(v[:,0]),1]

vb1 = vb.collections[0].get_paths()[0]
v1 = vb1.vertices
x2 = np.sort(v1[:,0])
y2 = v1[np.argsort(v1[:,0]),1]


def find_roots(x,y):
    s = np.abs(np.diff(np.sign(y))).astype(bool)
    return x[:-1][s] + np.diff(x)[s]/(np.abs(y[1:][s]/y[:-1][s])+1)

x = np.linspace(max(x1.min(), x2.min()), min(x1.max(), x2.max()), 1000)

y1i = np.interp(x, x1, y1 )  # blue
y2i = np.interp(x, x2, y2 )  # red


x_intersect = find_roots(x,y2i-y1i)
y_intersect = np.interp(x_intersect, x, y2i)

plt.plot(x_intersect, y_intersect, marker="X", ms=5, color="limegreen")
plt.show()

交点就是绿点。

当然，如果需要的话，需要对圆锥轮廓 (cnc.collections[0].get_paths()[0]) 的另一条臂执行相同的操作。

如何找到和绘制等高线形状的交点

How to find and draw the intersection points of contour shapes

matplotlib

contour

python-2.7