寻找直线和轮廓之间的交点

Question

我试图找到一条直线（红色虚线）与以红色突出显示的轮廓线的交点（见图表）。我在第二个图中使用 .get_paths 将所述等高线与其他图（第二个图）隔离开来。

我查看了轮廓相交问题How to find all the intersection points between two contour-set in an efficient way，并尝试将其用作基础，但无法重现任何有用的东西。

http://postimg.org/image/hz01fouvn/

http://postimg.org/image/m6utofwb7/

有人有什么想法吗？

重建绘图的相关函数，

#for contour 
def p_0(num,t) :
    esc_p = np.sum((((-1)**n)*(np.exp(t)**n)*((math.factorial(n)*((n+1)**0.5))**-1)) for n in range(1,num,1))
    return esc_p+1

tau = np.arange(-2,3,0.1)
r=[]

p1 = p_0(51,tau)
p2 = p_0(51,tau)

for i in p1:
    temp_r=i/p2
    r.append(temp_r)

x,y= np.meshgrid(tau,tau)
cs = plt.contour(x, y, np.log(r),50,colors='k')
whichContour =20
pa = CS.collections[whichContour].get_paths()[0]
v = pa.vertices
xx = v[:, 0]
yy = v[:, 1]
plt.plot(xx, yy, 'r-', label='Crossing contour')

#straight line 
p=0.75
logp = (np.log(p*np.exp(tau)))
plt.plot(tau,logp)

当前尝试，

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

def intercepting_line() :
matplotlib.rcParams['xtick.direction'] = 'out'
matplotlib.rcParams['ytick.direction'] = 'out'

#fake data

delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z1 = mlab.bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = mlab.bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
Z = 10.0 * (Z2 - Z1)

#plot
cs = plt.contour(X,Y,Z)
whichContour = 2 # change this to find the right contour lines

#get the vertices to calculate an intercept with a line
p = cs.collections[whichContour].get_paths()[0]
#see: http://matplotlib.org/api/path_api.html#module-matplotlib.path
v = p.vertices
xx = v[:, 0]
yy = v[:, 1]

#this shows the innermost ring now
plt.plot(xx, yy, 'r--', label='inner ring')

#fake line
x = np.arange(-2, 3.0, 0.1)
y=lambda x,m:(m*x)
y=y(x,0.9)
lineMesh = np.meshgrid(x,y)
plt.plot(x,y,'r' ,label='line')

#get the intercepts, two in this case 
x, y = find_intersections(v, lineMesh[1])
print x
print y
#plot the intercepting points
plt.plot(x[0], y[0], 'bo', label='first intercept')
#plt.plot(x[1], y[1], 'rs', label='second intercept')
plt.legend(shadow=True, fancybox=True, numpoints=1, loc='best')
plt.show()

#now we need to calculate the intercept of the vertices and whatever line
#this is pseudo code but works in case of two intercepting contour vertices

def find_intersections(A, B):
# min, max and all for arrays
amin = lambda x1, x2: np.where(x1<x2, x1, x2)
amax = lambda x1, x2: np.where(x1>x2, x1, x2)
aall = lambda abools: np.dstack(abools).all(axis=2)
slope = lambda line: (lambda d: d[:,1]/d[:,0])(np.diff(line, axis=0))

x11, x21 = np.meshgrid(A[:-1, 0], B[:-1, 0])
x12, x22 = np.meshgrid(A[1:, 0], B[1:, 0])
y11, y21 = np.meshgrid(A[:-1, 1], B[:-1, 1])
y12, y22 = np.meshgrid(A[1:, 1], B[1:, 1])
m1, m2 = np.meshgrid(slope(A), slope(B))
m1inv, m2inv = 1/m1, 1/m2

yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
xi = (yi - y21)*m2inv + x21

xconds = (amin(x11, x12) < xi, xi <= amax(x11, x12),
          amin(x21, x22) < xi, xi <= amax(x21, x22) )
yconds = (amin(y11, y12) < yi, yi <= amax(y11, y12),
          amin(y21, y22) < yi, yi <= amax(y21, y22) )

return xi[aall(xconds)], yi[aall(yconds)]

目前它找到相交点，但只在直线一致的地方，我在这里找不到解决方案的主要原因是我不明白原作者的思路，

yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
xi = (yi - y21)*m2inv + x21     

Answer 1

使用shapely可以找到交点，而不是使用该点作为fsolve()的初始猜测值来找到真正的解决方案：

#for contour 
def p_0(num,t) :
    esc_p = np.sum((((-1)**n)*(np.exp(t)**n)*((math.factorial(n)*((n+1)**0.5))**-1)) for n in range(1,num,1))
    return esc_p+1

tau = np.arange(-2,3,0.1)

x,y= np.meshgrid(tau,tau)
cs = plt.contour(x, y, np.log(p_0(51, y)/p_0(51, x)),[0.2],colors='k')

p=0.75
logp = (np.log(p*np.exp(tau)))
plt.plot(tau,logp)

from shapely.geometry import LineString
v1 = cs.collections[0].get_paths()[0].vertices

ls1 = LineString(v1)
ls2 = LineString(np.c_[tau, logp])
points = ls1.intersection(ls2)
x, y = points.x, points.y

from scipy import optimize

def f(p):
    x, y = p
    e1 = np.log(0.75*np.exp(x)) - y
    e2 = np.log(p_0(51, y)/p_0(51, x)) - 0.2
    return e1, e2

x2, y2 = optimize.fsolve(f, (x, y))

plt.plot(x, y, "ro")
plt.plot(x2, y2, "gx")

print x, y
print x2, y2

这是输出：

0.273616328952 -0.0140657435002
0.275317387697 -0.0123646847549

还有剧情：

Answer 2

将等高线视为折线，并将顶点坐标插入隐式线方程 (F(P) = a.X + b.Y + c = 0)。符号的每一次变化都是一个交点，通过求解 2x2 线性方程来计算。您不需要复杂的求解器。

如果您需要同时检测等高线，这并不复杂：考虑地形的剖面由通过该线的垂直平面。您将通过沿交叉网格图块边缘的线性插值获得高度。寻找与网格的交点与 Bresenham 画线算法密切相关。

那么你得到的是一个配置文件，即单个变量的函数。也可以通过检测符号的变化来定位与水平面的交点（iso-values）。

Answer 3

这是我用来解决这个问题的一种方法

def straight_intersection(straight1, straight2):
    p1x = straight1[0][0]
    p1y = straight1[0][1]
    p2x = straight1[1][0]
    p2y = straight1[1][1]
    p3x = straight2[0][0]
    p3y = straight2[0][1]
    p4x = straight2[1][0]
    p4y = straight2[1][1]
    x = p1y * p2x * p3x - p1y * p2x * p4x - p1x * p2y * p4x + p1x * p2y * p3x - p2x * p3x * p4y + p2x * p3y * p4x + p1x * p3x * p4y - p1x * p3y * p4x
    x = x / (p2x * p3y - p2x * p4y - p1x * p3y + p1x * p4y + p4x * p2y - p4x * p1y - p3x * p2y + p3x * p1y)
    y = ((p2y - p1y) * x + p1y * p2x - p1x * p2y) / (p2x - p1x)
    return (x, y)

Answer 4

虽然这个问题现在已经很老了，但我想为子孙后代分享我对这个问题的答案。 其实对于这个问题，我还发现了另外两种比较有效的解决方案。

第一个解决方案是像这样在 opencv 中递归使用pointPolygonTest。

# Enumerate the line segment points between 2 points
for pt in zip(*line(*p1, *p2)):
    if cv2.pointPolygonTest(conts[0], pt, False) == 0:  # If the point is on the contour
        return pt

第二种解决方案，你可以简单地画出轮廓和线，然后打个np.logical_and()就可以得到答案了

blank = np.zeros((1000, 1000))
blank_contour = drawContour(blank.copy(), cnt[0], 0, 1, 1)
blank_line = cv2.line(blank.copy(), line[0], line[1], 1, 1)
intersections = np.logical_and(blank_contour, black_line)

points = np.where(intersections == 1)