在二维混合数组中查找最近点

Question

为了给你们提供一些背景，我在做这个反演技术，我试图用积分值重现一个剖面。要做到这一点，我需要在数组中沿某一行找到值。为了举例说明我的问题，我有以下代码：

fig, ax = plt.subplots(1, figsize = (10,10))


#Create the grid (different grid spacing):
X = np.arange(0,10.01,0.25)
Y = np.arange(0,10.01,1.00)

#Create the 2D array to be plotted
Z = []
for i in range(np.size(X)):
    Zaux = []
    for j in range(np.size(Y)):
        Zaux.append(i*j + j)
        
        ax.scatter(X[i],Y[j], color = 'red', s = 0.25)
    Z.append(Zaux)

#Mesh the 1D grids:
Ymesh, Xmesh = np.meshgrid(Y, X)

#Plot the color plot:
ax.pcolor(Y,X,  Z, cmap='viridis', vmin=np.nanmin(Z), vmax=np.nanmax(Z))

#Plot the points in the grid of the color plot:
for i in range(np.size(X)):
    for j in range(np.size(Y)):
        ax.scatter(Y[j],X[i], color = 'red', s = 3)

#Create a set of lines:

for i in np.linspace(0,2,5):
    X_line = np.linspace(0,10,256)
    Y_line = i*X_line*3.1415-4

    #Plot each line:
    ax.plot(X_line,Y_line, color = 'blue')

ax.set_xlim(0,10)
ax.set_ylim(0,10)

plt.show()

将输出此图表：

我需要在Z中找到每一条线所穿过的最近点。其思想是将蓝色线条交叉的Z值进行积分，并将其绘制为线条斜率的函数。有人有好的解决办法吗？我试过一套for循环，但我觉得有点笨重

无论如何，谢谢你抽出时间

Answer 1

我不确定最近的点是什么。那看起来也很“笨重”。如果它正好在两点之间传递怎么办？另外，我已经为另一个项目编写了根据四个相邻像素的接近程度来衡量它们的代码，所以我将继续。我还冒昧地不重新缩放图片

i,j = np.meshgrid(np.arange(41),np.arange(11))
Z = i*j + j

class Image_knn():
    def fit(self, image):
        self.image = image.astype('float')

    def predict(self, x, y):
        image = self.image
        weights_x = [1-(x % 1), x % 1]
        weights_y = [1-(y % 1), y % 1]
        start_x = np.floor(x).astype('int')
        start_y = np.floor(y).astype('int')
        return sum([image[np.clip(np.floor(start_x + x), 0, image.shape[0]-1).astype('int'),
                          np.clip(np.floor(start_y + y), 0, image.shape[1]-1).astype('int')] * weights_x[x]*weights_y[y] 
                    for x,y in itertools.product(range(2),range(2))])

如果我们给它坐标，它就会返回图片

image_model = Image_knn()
image_model.fit(Z)
assert np.allclose(image_model.predict(*np.where(np.ones(Z.shape, dtype='bool'))).reshape((11,41)), Z)

我生成m=100条线，并缩放这些线上的点，使它们均匀分布。这是他们每10个人的一幅图。

n = 1000
m = 100
slopes = np.linspace(1e-10,10,m)
t, slope = np.meshgrid(np.linspace(0,1,n), slopes)
x_max, y_max = Z.shape[0]-1, Z.shape[1]-1
lines_x = t
lines_y = t*slope
scales = np.broadcast_to(np.stack([x_max/lines_x[:,-1], y_max/lines_y[:,-1]]).min(axis=0), (n,m)).T
lines_x *= scales
lines_y *= scales

最后我可以得到由斜率和积分组成的点，并画出它。你可能应该仔细看看“积分”，这只是我的猜测

%%timeit
points = np.array([(slope, np.mean(image_model.predict(lines_x[i],lines_y[i]))
                    *np.linalg.norm(np.array((lines_x[i,-1],lines_y[i,-1]))))
                   for i,slope in enumerate(slopes)])

plt.scatter(points[:,0],points[:,1])

注意最后一个块中的%%timeit。这在我的机器上占用了~38.3 ms，因此未进行优化。正如Donald Knuth所说，“过早优化是万恶之源”。如果要对此进行优化，您将删除for循环，通过重新塑造和重新塑造模型中直线点的所有坐标，然后使用坡度组织它们。但我看不出有什么理由把它扔了几分钟

最后我们得到了一个漂亮的尖头作为奖励。请注意，最大值为4是有意义的，因为对于40 x 10的图片，对角线的斜率为4。尖点的直觉有点难以解释，但我想你可能已经有了。对于长度，它归结为函数(x,y) -> sqrt(x^2+y^2)，在矩形上向上和向左时具有不同的方向差