寻找直线与轮廓的交点

Question

我正在尝试找到一条直线（红色虚线）与红色轮廓线的交点（见图）。在第二张图中，我使用了.get_paths来把这条轮廓线从其他线中分离出来。

我查阅了一个关于轮廓交点的问题，如何高效找到两个轮廓集合之间的所有交点，并尝试以此为基础，但没有成功得到任何有用的结果。

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

虽然这个问题现在已经有点时间了，但我想分享我的解决方案，以便将来的人能看到。

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

# 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)

Answer 2

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

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 3

把你的轮廓线看作是多条线段，然后把每个顶点的坐标代入一个隐式线方程（F(P) = a.X + b.Y + c = 0）。每当符号发生变化，就表示有一个交点，这个交点是通过解2x2的线性方程来计算的。你不需要复杂的求解器。

如果你想同时检测轮廓线，其实也没那么复杂：想象一下用一条垂直的平面切割地形。这样你就可以通过在交叉的网格边缘进行线性插值来获得高度。找到与网格的交点和Bresenham线段绘制算法有很大的关系。

最终你得到的是一个轮廓图，也就是一个单变量的函数。找到与水平面（等值线）的交点也是通过检测符号的变化来完成的。

Answer 4

使用 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

还有这个图：