如何在Python中插值两条直线之间的一条直线

Question

注意：我之前问过这个问题，但它被关闭为重复问题。然而，我和其他几个人认为它被错误地关闭了，我在我的原始帖子中解释了原因。所以我想在这里重新提出这个问题。

有没有人知道一个可以在两条直线之间进行插值的Python库？例如，给定下面的两条实线，我想要生成中间的虚线。换句话说，我想要得到中心线。输入仅为分别具有大小为N x 2和M x 2的坐标的两个numpy数组。

此外，我想知道是否有人已经在一些优化的Python库中编写了此函数。尽管优化并不是必需的。

以下是我可能拥有的两条线的示例，您可以假设它们彼此不重叠，并且x / y可以具有多个y / x坐标。

array([[ 1233.87375018,  1230.07095987],
       [ 1237.63559365,  1253.90749041],
       [ 1240.87500801,  1264.43925132],
       [ 1245.30875975,  1274.63795396],
       [ 1256.1449357 ,  1294.48254424],
       [ 1264.33600095,  1304.47893299],
       [ 1273.38192911,  1313.71468591],
       [ 1283.12411536,  1322.35942538],
       [ 1293.2559388 ,  1330.55873344],
       [ 1309.4817002 ,  1342.53074698],
       [ 1325.7074616 ,  1354.50276051],
       [ 1341.93322301,  1366.47477405],
       [ 1358.15898441,  1378.44678759],
       [ 1394.38474581,  1390.41880113]])

array([[ 1152.27115094,  1281.52899302],
       [ 1155.53345506,  1295.30515742],
       [ 1163.56506781,  1318.41642169],
       [ 1168.03497425,  1330.03181319],
       [ 1173.26135672,  1341.30559949],
       [ 1184.07110925,  1356.54121651],
       [ 1194.88086178,  1371.77683353],
       [ 1202.58908737,  1381.41765447],
       [ 1210.72465255,  1390.65097106],
       [ 1227.81309742,  1403.2904646 ],
       [ 1244.90154229,  1415.92995815],
       [ 1261.98998716,  1428.56945169],
       [ 1275.89219696,  1438.21626352],
       [ 1289.79440676,  1447.86307535],
       [ 1303.69661656,  1457.50988719],
       [ 1323.80994319,  1470.41028655],
       [ 1343.92326983,  1488.31068591],
       [ 1354.31738934,  1499.33260989],
       [ 1374.48879779,  1516.93734053],
       [ 1394.66020624,  1534.54207116]])

将其可视化后，我们有：

所以我尝试使用skimage.morphology库中的skeletonize函数，首先将坐标栅格化为填充多边形。 但是，我在末端遇到了分支，如下所示：

- jlcv

---

这可能是一个未明确说明的问题 - 即这里的“中间”如何定义？（或更糟糕的是，当曲线不像这样漂亮和凸起时。） - Oliver Charlesworth

这些线是否相交？任何 x 坐标是否有多个 y 坐标？ - user3483203

我们可以使用仅直线吗？ - IMCoins

我已经更新了问题以回答你的一些问题。关于“中间”的定义，我自己也不太确定如何定义。有多种解释方式吗？ - jlcv

Answer 1

首先，抱歉有些过头了；我很享受你的问题。如果说明太长，请随意跳到底部，我定义了一个可以完成我描述的所有操作的函数。

如果您的数组长度相同，那么您的问题将相对简单。在这种情况下，您只需要找到每个数组中对应x值和对应y值之间的平均值即可。

因此，我们可以创建长度相同的、更或多或少是原始数组的良好估计值的数组。我们可以通过对您拥有的数组进行多项式配合来实现这一点。正如评论和其他答案中所指出的那样，原始数组的中线没有明确定义，因此一个好的估计应该满足您的需求。

注意：在所有这些示例中，我已经命名了你发布的两个数组为a1和a2。

第一步：创建新数组以估算旧线路

查看您发布的数据：

这些不是特别复杂的函数，看起来一个三次多项式会非常适合它们。我们可以使用numpy创建它们：

import numpy as np

# Find the range of x values in a1
min_a1_x, max_a1_x = min(a1[:,0]), max(a1[:,0])
# Create an evenly spaced array that ranges from the minimum to the maximum
# I used 100 elements, but you can use more or fewer. 
# This will be used as your new x coordinates
new_a1_x = np.linspace(min_a1_x, max_a1_x, 100)
# Fit a 3rd degree polynomial to your data
a1_coefs = np.polyfit(a1[:,0],a1[:,1], 3)
# Get your new y coordinates from the coefficients of the above polynomial
new_a1_y = np.polyval(a1_coefs, new_a1_x)

# Repeat for array 2:
min_a2_x, max_a2_x = min(a2[:,0]), max(a2[:,0])
new_a2_x = np.linspace(min_a2_x, max_a2_x, 100)
a2_coefs = np.polyfit(a2[:,0],a2[:,1], 3)
new_a2_y = np.polyval(a2_coefs, new_a2_x)

结果：

这并不太糟糕！如果您有更复杂的函数，您将需要拟合更高次数的多项式或找到其他适合拟合您的数据的函数。

现在，您有两组长度相同的数组（我选择了长度为100，您可以根据您希望中点线平滑程度来决定使用更多或更少的数据）。这些集合代表了原始数组的估计的x和y坐标。在上面的示例中，我将它们命名为new_a1_x、new_a1_y、new_a2_x和new_a2_y。

第二步：计算每个新数组中每个x和每个y的平均值

接下来，我们要为每个估计数组找到平均x和平均y值。只需使用np.mean：

midx = [np.mean([new_a1_x[i], new_a2_x[i]]) for i in range(100)]
midy = [np.mean([new_a1_y[i], new_a2_y[i]]) for i in range(100)]

midx和midy现在表示我们两个估计数组之间的中点。现在，只需要将您的原始（非估计）数组与中点数组一起绘制：

plt.plot(a1[:,0], a1[:,1],c='black')
plt.plot(a2[:,0], a2[:,1],c='black')
plt.plot(midx, midy, '--', c='black')
plt.show()

看这里：

即可得到结果。

这种方法仍然适用于更复杂的、嘈杂的数据（但必须要谨慎地拟合函数）：

作为函数：

我将上面的代码放在一个函数中，这样你就可以轻松使用它。它会返回一个数组，其中包含你估计的中点，格式与你原始的数组相同。

参数：a1和a2是你的两个输入数组，poly_deg是你想要拟合的多项式的次数，n_points是你想要在中点数组中的点数，plot是一个布尔值，表示你是否想要绘图。

import matplotlib.pyplot as plt
import numpy as np

def interpolate(a1, a2, poly_deg=3, n_points=100, plot=True):

    min_a1_x, max_a1_x = min(a1[:,0]), max(a1[:,0])
    new_a1_x = np.linspace(min_a1_x, max_a1_x, n_points)
    a1_coefs = np.polyfit(a1[:,0],a1[:,1], poly_deg)
    new_a1_y = np.polyval(a1_coefs, new_a1_x)

    min_a2_x, max_a2_x = min(a2[:,0]), max(a2[:,0])
    new_a2_x = np.linspace(min_a2_x, max_a2_x, n_points)
    a2_coefs = np.polyfit(a2[:,0],a2[:,1], poly_deg)
    new_a2_y = np.polyval(a2_coefs, new_a2_x)

    midx = [np.mean([new_a1_x[i], new_a2_x[i]]) for i in range(n_points)]
    midy = [np.mean([new_a1_y[i], new_a2_y[i]]) for i in range(n_points)]

    if plot:
        plt.plot(a1[:,0], a1[:,1],c='black')
        plt.plot(a2[:,0], a2[:,1],c='black')
        plt.plot(midx, midy, '--', c='black')
        plt.show()

    return np.array([[x, y] for x, y in zip(midx, midy)])

[编辑]:

我回想起这个问题，发现一个更简单的方法，就是使用np.interp使两个数组都变得密集。这种方法与上面的拟合直线方法基本相似，但不使用polyfit/polyval来近似拟合线段，而是仅仅使其变得更加密集:

min_a1_x, max_a1_x = min(a1[:,0]), max(a1[:,0])
min_a2_x, max_a2_x = min(a2[:,0]), max(a2[:,0])

new_a1_x = np.linspace(min_a1_x, max_a1_x, 100)
new_a2_x = np.linspace(min_a2_x, max_a2_x, 100)

new_a1_y = np.interp(new_a1_x, a1[:,0], a1[:,1])
new_a2_y = np.interp(new_a2_x, a2[:,0], a2[:,1])

midx = [np.mean([new_a1_x[i], new_a2_x[i]]) for i in range(100)]
midy = [np.mean([new_a1_y[i], new_a2_y[i]]) for i in range(100)]

plt.plot(a1[:,0], a1[:,1],c='black')
plt.plot(a2[:,0], a2[:,1],c='black')
plt.plot(midx, midy, '--', c='black')
plt.show()

- sacuL

---

@sacul 你知道C#的等价物吗？我发现Math.NET包有大部分，如果不是全部的函数，但它们并不完全相同，所以我不确定哪一个可以转换。 - jbassking10

抱歉，我在C#方面不是很擅长，但从Math.NET文档中看来，MathNet.Numerics中的PolynomialFunc类型的Fit可能能够处理许多函数拟合问题，特别是Polynomial函数，它基本上可以替换我的代码中的np.polyfit。其余部分比较简单，但我在那方面不是一个好的帮助来源... - sacuL

好的。那真是太有帮助了！谢谢。 - jbassking10

Answer 2

“两条线之间的界限”并不十分明确。您可以通过在两条曲线之间进行三角形化来获得一个简单而合理的解决方案（您可以通过从顶点到顶点前进，选择产生较少倾斜三角形的对角线来进行三角形化）。

然后，插值曲线将连接两侧的中心部分。

- user1196549

---

嗯，这可能有效，我只需在每行中加密点以获得更精确的解决方案。你所说的“不那么偏斜的三角形”是什么意思？ - jlcv

@JustinLiang：如果这些曲线的方程式已知，答案将会不同。 - user1196549

@JustinLiang “我只需要加密点”：中位曲线不需要比给定的曲线更准确/平滑。因此，没有必要加密。 - user1196549

Answer 3

我从事河流工作，所以这是一个常见的问题。我的其中一个解决方案与您在问题中展示的完全相同——即使骨架化这个斑点。您会发现边界有问题，所以我所做的看起来效果不错，只是简单地镜像边界。为了使这种方法有效，斑点不能与图像的角落相交。

您可以在 RivGraph 中找到我的实现；这个特定的算法在 rivers/river_utils.py 中被称为 "mask_to_centerline"。

以下是一个示例输出，显示中心线的端点如何延伸到对象的期望边缘：

- Jon

Answer 4

sacuL的解决方案对我几乎起作用了，但我需要汇总不止两个曲线。这是我对sacuL解决方案的概括：

def interp(*axis_list):
    min_max_xs = [(min(axis[:,0]), max(axis[:,0])) for axis in axis_list]

    new_axis_xs = [np.linspace(min_x, max_x, 100) for min_x, max_x in min_max_xs]
    new_axis_ys = [np.interp(new_x_axis, axis[:,0], axis[:,1]) for axis, new_x_axis in zip(axis_list, new_axis_xs)]

    midx = [np.mean([new_axis_xs[axis_idx][i] for axis_idx in range(len(axis_list))]) for i in range(100)]
    midy = [np.mean([new_axis_ys[axis_idx][i] for axis_idx in range(len(axis_list))]) for i in range(100)]

    for axis in axis_list:
        plt.plot(axis[:,0], axis[:,1],c='black')
    plt.plot(midx, midy, '--', c='black')
    plt.show()

如果我们现在运行一个例子：

a1 = np.array([[x, x**2+5*(x%4)] for x in range(10)])
a2 = np.array([[x-0.5, x**2+6*(x%3)] for x in range(10)])
a3 = np.array([[x+0.2, x**2+7*(x%2)] for x in range(10)])
interp(a1, a2, a3)

我们得到了这个情节：。

- yosef