контур разреженных (регулярных) данных показывает артефакты

Question

Я хотел бы очертить данные, которые являются довольно редкими и где максимум проходит по диагонали через изображение;контурная функция matplotlib изобретает минимумы между выборочными максимумами.

Начиная с плотно отобранного случая, когда все выглядит так, как ожидалось:

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

x_1D = np.linspace(0., 10., 100)
y_1D = np.linspace(0., 10., 100)
x, y = np.meshgrid(x_1D, y_1D)
z = np.empty_like(x)
def peak(y, y0):
    return np.exp(-(y-y0)**2)
for i in range(x_1D.size):
    z[:,i] = peak(y_1D, i/x_1D.size*y_1D.max())

fig, ax = plt.subplots(ncols=3)

ax[0].set_title('measured data')
ax[0].scatter(x, y, marker='s', c=z, cmap=plt.cm.jet, s=25)

ax[1].set_title('contour')
ax[1].contourf(x, y, z, levels=14, cmap=plt.cm.jet)

# define grid
xi = np.linspace(x_1D.min()-0.1, x_1D.max()+0.1, 1000)
yi = np.linspace(y_1D.min()-0.1, y_1D.max()+0.1, 1000)

# grid the data
triang = tri.Triangulation(x.flatten(), y.flatten())
interpolator = tri.LinearTriInterpolator(triang, z.flatten())
Xi, Yi = np.meshgrid(xi, yi)
zi = interpolator(Xi, Yi)

ax[2].set_title('interpolated')
ax[2].contourf(xi, yi, zi, levels=14, cmap=plt.cm.jet)

plt.show()

выход

Когда x дискретизируется меньше с коэффициентом 10, то есть x_1D = np.linspace(0., 10., 10), минимумы появляются между выборочными максимумами на контурном графике.

Есть ли способ, как избежать этого артефакта и сделать контур разреженных данных похожим на контур данных с плотной выборкой?

РЕДАКТИРОВАТЬ: Спасибо за ответ, который работает оченьхорошо на примере, который я привел.К сожалению, я слишком упростил проблему.Вместо того, чтобы говорить об одной диагональной линии, я должен был узнать о произвольном количестве пиков, движущихся в произвольных направлениях по всей картине;например, замените генерацию пика на

z = np.zeros_like(x)
def peak(y, y0):
    return np.exp(-(y-y0)**2)
for i in range(x_1D.size):
    z[:,i] += peak(y_1D, np.cos(i/x_1D.size*np.pi)*y_1D.max()*0.05+y_1D.max()*0.8)
for i in range(x_1D.size):
    z[:,i] += peak(y_1D, np.sin(i/x_1D.size*np.pi/2.)*y_1D.max()*0.5)

, что приведет к

Answer 1

Основная проблема вашего подхода заключается в том, что алгоритм триангуляции не знает, что пики должны соединяться друг с другом между «x-кусочками» (ваша линия точек плотных данных для постоянной x).

Немного упрощая, алгоритм триангуляции будет смотреть на соседей в направлении x и y и соединяться с ними. Затем, при попытке интерполировать с использованием этой триангуляции, точки между пиками будут примерно равны средним значениям ближайших точек в направлении x, и, следовательно, появятся минимумы. Лучшее решение - сделать свою собственную триангуляцию, подключив пики напрямую.

К счастью, мы можем на самом деле взломать триангуляцию, чтобы соединить ее с пиками, сместив координаты в направлении y так, чтобы все пики были выровнены по горизонтали. Это работает, потому что алгоритм триангуляции использует координаты, которые вы передаете. В вашем примере это легко сделать, потому что мы можем просто применить простой сдвиг y_s = y - x. Как правило, вам нужно получить уравнение для вашего пика (назовите его y_p(x)), а затем вычесть это из y, чтобы получить y_s.

Теперь, когда у вас есть смещенная триангуляция, вы можете создать новую более плотную сетку (как вы сделали) и применить тот же сдвиг. Затем вы интерполируете в смещенной сетке смещенную плотную сетку, чтобы получить значения z, правильно интерполированные. Наконец, вы снимаете плотную сетку, чтобы получить правильные значения y, и строите ее.

Ниже приведен код применения этой концепции к вашему коду и окончательный результат. Как вы видете. Это хорошо работает для этого случая.

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

def peak(y, y0):
    return np.exp(-(y-y0)**2)

x_1D = np.linspace(0., 10., 10)
y_1D = np.linspace(0., 10., 100)
x, y = np.meshgrid(x_1D, y_1D)
z = np.empty_like(x)

for i in range(x_1D.size):
    z[:,i] = peak(y_1D, i/x_1D.size*y_1D.max())

fig, ax = plt.subplots(ncols=4)

ax[0].set_title('measured data')
ax[0].scatter(x, y, marker='s', c=z, cmap=plt.cm.jet, s=25)

ax[1].set_title('contour')
ax[1].contourf(x, y, z, levels=14, cmap=plt.cm.jet)

# define output grid
xi_1D = np.linspace(x_1D.min()-0.1, x_1D.max()+0.1, 1000)
yi_1D = np.linspace(y_1D.min()-0.1, y_1D.max()+0.1, 1000)
xi, yi = np.meshgrid(xi_1D, yi_1D)

# Old Linear Interpolation
triang = tri.Triangulation(x.flatten(), y.flatten())
interpolator = tri.LinearTriInterpolator(triang, z.flatten())
zi = interpolator(xi, yi)

ax[2].set_title('interpolated')
ax[2].contourf(xi, yi, zi, levels=14, cmap=plt.cm.jet)

# === SHIFTED LINEAR INTERPOLATION ===

# make shifted interpolating mesh for the data
y_s=y-x
triang_s = tri.Triangulation(x.flatten(), y_s.flatten())
interpolator_s = tri.LinearTriInterpolator(triang_s, z.flatten())

# interpolate in the shifted state
yi_s = yi-xi
zi_s = interpolator_s(xi, yi_s)

# unshift the fine mesh
yi_us = yi_s+xi

ax[3].set_title('interpolated (shifted)')
ax[3].contourf(xi, yi_us, zi_s, levels=14, cmap=plt.cm.jet)


plt.show()