Проблемы с 2D-интерполяцией в Scipy

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

from __future__ import division

import numpy as np
import scipy.interpolate as intp
import pylab as plt

# Defining some variables:

quadratic = -3/128
linear = 1/16
pn = np.poly1d([quadratic, linear,0])

pixels_x = 50
pixels_y = 30
frame = np.zeros((pixels_x,pixels_y))

x_width= np.concatenate((np.linspace(8,7.8,57) , np.linspace(7.8,8,pixels_y-57)))

def data(x,y):
    z = y*(np.exp(-(x-5)**2/3) + np.exp(-(x)**2/5) + np.exp(-(x+5)**2))
    return(z)

# Generating grid coordinates

yt = np.arange(380,380+pixels_y*4,4)
xt = np.linspace(-7.8,7.8,pixels_x)
X, Y = np.meshgrid(xt,yt)
Y=Y.T
X=X.T

Y_m = np.zeros((pixels_x,pixels_y))
X_m = np.zeros((pixels_x,pixels_y))

# generating distorted grid coordinates:    

for i in range(pixels_y):
    Y_m[:,i] = Y[:,i] - pn(xt)
    X_m[:,i] = np.linspace(-x_width[i],x_width[i],pixels_x)


# Sample data:
for i in range(pixels_y):
    for j in range(pixels_x):
        frame[j,i] = data(X_m[j,i],Y_m[j,i])


Y_m = Y_m.flatten()
X_m = X_m.flatten()
frame = frame.flatten()
##
Y = Y.flatten()
X = X.flatten()
ipf = intp.interp2d(X_m,Y_m,frame)
interpolated_frame = ipf(xt,yt)

На данный момент у меня есть вопросы:

  1. Код работает, но я получаю следующее предупреждение:

    Предупреждение: больше нельзя добавлять узлы, потому что количество коэффициентов B-сплайна уже превышает количество точек данных m. Вероятные причины: либо s, либо m слишком малы. (fp>s) kx,ky=1,1 nx,ny=54,31 m=1500 fp=0,000006 s=0,000000

Также появляются некоторые артефакты интерполяции, и я предполагаю, что они связаны с предупреждением - Ребята, вы знаете, что я делаю неправильно?

  1. Для моих реальных приложений кадры должны быть около 500 * 100, но при этом я получаю MemoryError. Могу ли я что-то сделать, чтобы помочь этому, кроме разделения кадра на несколько частей?

Спасибо!


person Dzz    schedule 10.10.2010    source источник


Ответы (3)


Эта проблема, скорее всего, связана с использованием bisplrep и bisplev в interp2d. Документы упоминают, что они используют коэффициент сглаживания s=0.0 и что bisplrep и bisplev следует использовать напрямую, если больше контроля над < strong>s необходим. Связанные документы упоминают, что s следует искать между (m-sqrt(2*m),m+sqrt(2*m)), где m — количество точек, используемых для построения сплайнов. У меня была похожая проблема, и я нашел ее решенной при непосредственном использовании bisplrep и bisplev, где s является необязательным.

person Andre    schedule 17.02.2016
comment
но это линейная интерполяция, почему так называются узлы сплайна? - person SBFRF; 04.08.2020
comment
как экстраполировать за пределы домена xy, используя bisplrep и bisplev? - person giammi56; 11.06.2021

Для двухмерной интерполяции griddata являются надежными, локальными и быстрыми. Взгляните на проблема-with-2d-interpolation-in- scipy-non-rectangular-grid на SO.

person denis    schedule 20.10.2010

Возможно, вы захотите посмотреть на следующий метод интерполяции в базовой карте:

mpl_toolkits.basemap.interp http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html

если вам действительно не нужна интерполяция на основе сплайнов.

person reckoner    schedule 10.10.2010
comment
Извините, это не поможет - mpl_toolkits.basemap.interp предназначен только для интерполяции одной прямолинейной сетки в другую. Этому условию моя входная сетка не соответствует. - person Dzz; 11.10.2010