Интерполяция по регулярной сетке в Python

Я изо всех сил пытался внедрить данные для "пустых" пикселей в моей 2D-матрице. В основном, я понимаю (но не глубоко) методы интерполяции, такие как Inverse Distance Weighting, Kriging, Bicubic и т.д. Я точно не знаю начальную точку (либо в постановке проблемы, либо в случае Python).

Определение проблемы: У меня есть матрица MxN (регулярная сетка), в которой каждый пиксель представляет определенное значение измерения (ниже и данные, используемые на этом рисунке здесь). Я хотел бы интерполировать данные для "пространства вопросительных знаков" (пустое пространство, которое также состоит из одинаковых размеров, но пустых пикселей), используя существующие данные, которые у меня есть как синие пиксели.

1) Как я могу интерполировать эти данные. Может ли кто-нибудь дать мне простой пример (например, матрицу 3x3), чтобы понять это ясно?

2) Может ли кто-нибудь указать мне, как выполнять шаги по решению в среде Python?

3) Как сравнить методы интерполяции в смысле точности с помощью Python?

4) Считаете ли вы целесообразным использовать различную интерполяцию в зависимости от плотности данных?

Буду признателен за ваши ответы и предложения.

Ответы

Ответ 1

Какое разумное решение во многом зависит от того, на какие вопросы вы пытаетесь ответить с помощью интерполированных пикселей - caveat emptor: экстраполяция по отсутствующим данным может привести к очень вводящим в заблуждение ответам!

Радиальная базовая функция Интерполяция/сглаживание ядра

С точки зрения практических решений, доступных в Python, одним из способов заполнения этих пикселей было бы использование Scipy реализации интерполяции Radial Basis Function (см. здесь), которая предназначенные для сглаживания/интерполяции рассеянных данных.

Учитывая вашу матрицу M и лежащие в основе 1D координатные массивы r и c (такие, что M.shape == (r.size, c.size)), где отсутствующие записи M установлены на nan, это, похоже, работает довольно хорошо с линейной RBF следующим образом:

import numpy as np
import scipy.interpolate as interpolate

with open('measurement.txt') as fh:
    M = np.vstack(map(float, r.split(' ')) for r in fh.read().splitlines())
r = np.linspace(0, 1, M.shape[0]) 
c = np.linspace(0, 1, M.shape[1])

rr, cc = np.meshgrid(r, c)
vals = ~np.isnan(M)
f = interpolate.Rbf(rr[vals], cc[vals], M[vals], function='linear')
interpolated = f(rr, cc)

Это дает следующую интерполяцию данных, которые вы связали выше, что, хотя и разумно выглядит, показывает, насколько неблагоприятным является отношение недостающих образцов к реальным данным:

Гауссовская регрессия процесса/кригинг

Интерпретация кригинга доступна с помощью регрессии Gaussian Process(), которая сама по себе основана на инструменте DACE Kriging для Matlab) в библиотеке scikit-learn. Это можно вызвать следующим образом:

from sklearn.gaussian_process import GaussianProcess

gp = GaussianProcess(theta0=0.1, thetaL=.001, thetaU=1., nugget=0.01)
gp.fit(X=np.column_stack([rr[vals],cc[vals]]), y=M[vals])
rr_cc_as_cols = np.column_stack([rr.flatten(), cc.flatten()])
interpolated = gp.predict(rr_cc_as_cols).reshape(M.shape)

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

Inpainting

В качестве окончательного в стороне, полностью визуально мотивированное решение будет использовать функциональность OpenCV inpainting, хотя это предполагает 8-битные массивы (0 - 255) и не имеет простой математической интерпретации.