如何生成等间距插值答案

【问题标题】：How to generate equispaced interpolating values如何生成等间距插值
【发布时间】：2013-10-01 13:35:18
【问题描述】：

我有一个不均匀间隔的 (x,y) 值列表。 Here 是本题使用的存档。

我能够在值之间进行插值，但我得到的不是等间距的插值点。这是我的工作：

x_data = [0.613,0.615,0.615,...]
y_data = [5.919,5.349,5.413,...]

# Interpolate values for x and y.
t = np.linspace(0, 1, len(x_data))
t2 = np.linspace(0, 1, 100)
# One-dimensional linear interpolation.
x2 = np.interp(t2, t, x_data)
y2 = np.interp(t2, t, y_data)

# Plot x,y data.
plt.scatter(x_data, y_data, marker='o', color='k', s=40, lw=0.)

# Plot interpolated points.
plt.scatter(x2, y2, marker='o', color='r', s=10, lw=0.5)

结果：

可以看出，在图中原始点分布更密集的部分，红点更靠近。

我需要一种方法来根据给定的步长值（比如 0.1）在 x、y 中生成插值等距

正如 askewchan 正确指出的那样，当我的意思是“等距 in x, y”时，我的意思是曲线中的两个连续插值点应该彼此远离（欧几里得直线距离）由相同的值。

我尝试了 unubtu 的答案，它适用于平滑曲线，但似乎不适合不太平滑的曲线：

发生这种情况是因为代码以欧几里得方式而不是直接在曲线上计算点距离，并且我需要 曲线上 的距离在点之间相同。能否以某种方式解决此问题？

【问题讨论】：

“equispaced in x, y”是否表示沿曲线的切线等距？
@askewchan 对不起，我表达得模棱两可。我将重构问题以使其更清楚。
换句话说，你想要dx = step / sqrt(1 + (y')**2)
您需要将插值构建为函数y(x)，然后在某个数组x 上调用它，该数组不是线性间隔的，而是与 x (dx) 中的每一步间隔为您的所需的步长（例如 ds）除以（1 + 斜率平方）的 sqrt。这是因为曲线段的弧长的长度为 ds = sqrt(1 + slope^2)*dx（来自勾股定理 ds^2 = dx^2 + dy^2）
"...应该彼此保持距离..."，如何测量距离？是直线还是沿曲线的距离？

标签： python numpy scipy interpolation

【解决方案1】：

将您的 xy 数据转换为参数化曲线，即计算点之间的所有距离并通过累积求和生成曲线上的坐标。然后相对于新坐标独立地插值 x 和 y 坐标。

import numpy as np
from matplotlib import pyplot as plt

data = '''0.615   5.349
    0.615   5.413
    0.617   6.674
    0.617   6.616
    0.63    7.418
    0.642   7.809
    0.648   8.04
    0.673   8.789
    0.695   9.45
    0.712   9.825
    0.734   10.265
    0.748   10.516
    0.764   10.782
    0.775   10.979
    0.783   11.1
    0.808   11.479
    0.849   11.951
    0.899   12.295
    0.951   12.537
    0.972   12.675
    1.038   12.937
    1.098   13.173
    1.162   13.464
    1.228   13.789
    1.294   14.126
    1.363   14.518
    1.441   14.969
    1.545   15.538
    1.64    16.071
    1.765   16.7
    1.904   17.484
    2.027   18.36
    2.123   19.235
    2.149   19.655
    2.172   20.096
    2.198   20.528
    2.221   20.945
    2.265   21.352
    2.312   21.76
    2.365   22.228
    2.401   22.836
    2.477   23.804'''

data = np.array([line.split() for line in data.split('\n')],dtype=float)

x,y = data.T
xd = np.diff(x)
yd = np.diff(y)
dist = np.sqrt(xd**2+yd**2)
u = np.cumsum(dist)
u = np.hstack([[0],u])

t = np.linspace(0,u.max(),10)
xn = np.interp(t, u, x)
yn = np.interp(t, u, y)

f = plt.figure()
ax = f.add_subplot(111)
ax.set_aspect('equal')
ax.plot(x,y,'o', alpha=0.3)
ax.plot(xn,yn,'ro', markersize=8)
ax.set_xlim(0,5)

【讨论】：

就我而言，它在这一行给了我ValueError: object too deep for desired array - xn = numpy.interp(t, u, x)
我刚刚运行了代码，它运行良好。 numpy 是 1.10.1 版本和 scipy 0.15.1
谢谢你，这正是我也在寻找的。关于如何将其应用于更高维数据（可能使用scipy.interpolate.interpn）的任何想法？我宁愿不遍历每个维度并独立执行np.interp。
在此评论后 10 分钟，我找到了解决方案:)。在这里作为答案发布

【解决方案2】：

让我们首先考虑一个简单的案例。假设您的数据看起来像蓝线，以下。

如果您想选择距离为r 的等距点，那么r 会有一些临界值，其中 (1,2) 处的尖点是第一个等距点。

如果您希望点之间的距离大于这个临界距离，那么第一个等距点会从 (1,2) 跳到非常不同的某个地方—— 由绿色弧线与蓝色线的交点描绘。这种变化不是渐进的。

这个玩具案例表明，参数r 的微小变化可能会对解产生根本的、不连续的影响。

也建议你必须知道第i个等距点的位置在您确定第 (i+1) 个等距点的位置之前。

所以看来需要一个迭代解决方案：

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

x, y = np.genfromtxt('data', unpack=True, skip_header=1)
# find lots of points on the piecewise linear curve defined by x and y
M = 1000
t = np.linspace(0, len(x), M)
x = np.interp(t, np.arange(len(x)), x)
y = np.interp(t, np.arange(len(y)), y)
tol = 1.5
i, idx = 0, [0]
while i < len(x):
    total_dist = 0
    for j in range(i+1, len(x)):
        total_dist += math.sqrt((x[j]-x[j-1])**2 + (y[j]-y[j-1])**2)
        if total_dist > tol:
            idx.append(j)
            break
    i = j+1

xn = x[idx]
yn = y[idx]
fig, ax = plt.subplots()
ax.plot(x, y, '-')
ax.scatter(xn, yn, s=50)
ax.set_aspect('equal')
plt.show()

注意：我将纵横比设置为'equal'，以便更明显地看出这些点是等距的。

【讨论】：

这是一个很好的答案，抱歉，我花了这么长时间才回复您（遇到了另一个问题）此代码适用于相当平滑的曲线，但似乎因表现不佳而失败曲线。请查看已编辑的问题。
@Gabriel：我已经编辑了我的答案。基本相同 - 您只需在点之间跳跃时跟踪总距离。
非常感谢@unubtu，这个答案现在很完美！

【解决方案3】：

以下脚本将以x_max - x_min / len(x) = 0.04438的相等步长插入点

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt

data = np.loadtxt('data.txt')
x = data[:,0]
y = data[:,1]

f = interp1d(x, y)
x_new = np.linspace(np.min(x), np.max(x), x.shape[0])
y_new = f(x_new)

plt.plot(x,y,'o', x_new, y_new, '*r')
plt.show()

【讨论】：

此解决方案生成在 x 上等距但在 y 上不等距的点，因此在曲线中不等距。还是谢谢你的回答。

【解决方案4】：

扩展@Christian K. 的答案，以下是使用scipy.interpolate.interpn 对更高维数据执行此操作的方法。假设我们要重新采样到 10 个等距点：

import numpy as np
import scipy
# Assuming that 'data' is rows x dims (where dims is the dimensionality)
diffs = data[1:, :] - data[:-1, :]
dist = np.linalg.norm(diffs, axis=1)
u = np.cumsum(dist)
u = np.hstack([[0], u])
t = np.linspace(0, u[-1], 10)
resampled = scipy.interpolate.interpn((u,), pts, t)

【讨论】：

【解决方案5】：

可以沿曲线生成等距点。但是必须有更多的定义你想要一个真正的答案。抱歉，我为此任务编写的代码是在 MATLAB 中，但我可以描述一般的想法。有三种可能。

首先，就简单的欧几里得距离而言，这些点与邻居的距离是否真正等距？这样做需要找到曲线上任意点与固定半径圆的交点。然后沿着曲线走。

接下来，如果您打算将距离表示为沿曲线本身的距离，如果曲线是分段线性的，那么问题又很容易解决。只需沿着曲线走，因为线段上的距离很容易测量。

最后，如果您打算让曲线成为三次样条曲线，这同样不是非常困难，但需要更多的工作。这里的诀窍是：

计算沿曲线从点到点的分段线性弧长。就叫它吧。
生成一对三次样条，x(t), y(t)。
区分 x 和 y 作为 t 的函数。由于这些是立方段，这很容易。导数函数将是分段二次的。
使用 ode 求解器沿曲线移动，集成微分弧长函数。在 MATLAB 中，ODE45 运行良好。

因此，一个集成

sqrt((x')^2 + (y')^2)

同样，在 MATLAB 中，可以设置 ODE45 来识别函数穿过某些指定点的位置。

如果您的 MATLAB 技能能够胜任这项任务，您可以查看 interparc 中的代码以获取更多说明。它是相当好的注释代码。

【讨论】：