Как создать равные интерполяционные значения

У меня есть список значений (x, y), которые не равномерно распределены. Здесь - архив, используемый в этом вопросе.

Я могу интерполировать между значениями, но то, что я получаю, не является равным интерполирующим точкам. Вот что я делаю:

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)

Результат:

Как видно, красные точки ближе друг к другу в сечениях графика, где распределение исходных точек плотнее.

Мне нужен способ генерировать интерполированные точки, равные x, y согласно заданному значению шага (скажем 0,1)

Как правильно указывает askewchan, когда я имею в виду "равномерно в x, y", я имею в виду, что две последовательные интерполированные точки на кривой должны быть дистанцированы друг от друга (эвклидовое расстояние по прямой) то же значение.

Я попытался ответить unubtu, и он хорошо работает для гладких кривых, но, похоже, сломается не так гладко:

Это происходит потому, что код вычисляет расстояние точки в евклидовом направлении вместо прямой над кривой, и мне нужно, чтобы расстояние по кривой было одинаковым между точками. Можно ли каким-то образом решить эту проблему?

Ответы

Ответ 1

Рассмотрим сначала простой случай. Предположим, что ваши данные выглядели как синяя линия, ниже.

Если вы хотите выбрать эквидистантные точки, расположенные на расстоянии 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', чтобы сделать более очевидным, что точки эквидистантны.

Ответ 2

Преобразуйте xy-данные в параметризованную кривую, т.е. вычислите все расстояния между точками и сформируйте координаты на кривой путем суммирования. Затем интерполируем x- и y-координаты независимо друг от друга относительно новых координат.

import numpy as np
from pylab import plot

data = '''    0.613   5.919
    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(),20)
xn = np.interp(t, u, x)
yn = np.interp(t, u, y)

plot(x,y,'o')
plot(xn,yn,'gs')
xlim(0,5.5)
ylim(10,17.5)

Ответ 3

Следующий script будет интерполировать точки с равным шагом 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()

Ответ 4

Возможно создание эквидистантных точек вдоль кривой. Но должно быть больше определения того, что вы хотите для реального ответа. Извините, но код, который я написал для этой задачи, находится в MATLAB, но я могу описать общие идеи. Существует три возможности.

Во-первых, точки, которые должны быть действительно равноудалены от соседей в терминах простого евклидова расстояния? Для этого потребуется найти пересечение в любой точке кривой с кругом фиксированного радиуса. Затем просто шагните по кривой.

Затем, если вы намереваетесь расстояние до среднего расстояния по самой кривой, если кривая является кусочно-линейной, проблема снова легко сделать. Просто шагните по кривой, так как расстояние на отрезке линии легко измерить.

Наконец, если вы намереваетесь, чтобы кривая была кубическим сплайном, снова это не невероятно сложно, но это немного больше работы. Вот трюк:

Вычислить кусочно-линейную длину дуги от точки к точке вдоль кривой. Назовите его t.
Создайте пару кубических сплайнов, x (t), y (t).
Дифференцируем x и y как функции от t. Поскольку это кубические сегменты, это легко. Производные функции будут кусочно квадратичными.
Используйте ode solver для перемещения по кривой, интегрируя дифференциальную функцию arclength. В MATLAB ODE45 работал красиво.

Таким образом, каждый интегрирует

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

Опять же, в MATLAB ODE45 можно установить для идентификации тех мест, где функция пересекает определенные указанные точки.

Если ваши навыки MATLAB соответствуют задаче, вы можете посмотреть код interparc для получения дополнительных пояснений. Это достаточно хорошо прокомментированный код.