Как интерполировать в повернутую сетку?

Question

У меня есть

[[0, 1, 2, 3, 4],
 [1, 2, 3, 4, 5],
 [2, 3, 4, 5, 6],
 [3, 4, 5, 6, 7],
 [4, 5, 6, 7, 8]]

, и я хочу интерполировать его в повернутую сетку с углом на левом краю.Аналогично:

[[2, ~2, 2],
 [~4, ~4, ~4],
 [6, ~6, 6]]

(я использую ~ для обозначения приблизительных значений.)

(Конечно, мои фактические данные более сложные. Сценарий таков, что я хочу отобразить DEMданные по пикселям на повернутое изображение.)

Вот установка:

import numpy
from scipy import interpolate as interp

grid = numpy.ndarray((5, 5))
for I in range(grid.shape[0]):
    for j in range(grid.shape[1]):
        grid[I, j] = I + j

grid = ndimage.interpolation.shift(
    ndimage.interpolation.rotate(grid, -45, reshape=False),
    -1)

source_x, source_y = numpy.meshgrid(
    numpy.arange(0, 5), numpy.arange(0, 5))
target_x, target_y = numpy.meshgrid(
    numpy.arange(0, 2), numpy.arange(0, 2))

print(interp.griddata(
    numpy.array([source_x.ravel(), source_y.ravel()]).T,
    grid.ravel(),
    target_x, target_y)) 

Это дает мне:

[[2.4467   2.6868 2.4467]
 [4.       4.     4.    ]
 [5.5553   5.3132 5.5553]]

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

Я знаю индексы углов сетки, к которым я хочу интерполировать,То есть у меня

upper_left = 2, 0
upper_right = 0, 2
lower_right = 4, 2
lower_left = 2, 4

Answer 1

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

import numpy as np
from scipy.interpolate import RectBivariateSpline

# input data
data = np.array([[0, 1, 2, 3, 4],
                 [1, 2, 3, 4, 5],
                 [2, 3, 4, 5, 6],
                 [3, 4, 5, 6, 7],
                 [4, 5, 6, 7, 8]])

upper_left = 2, 0
upper_right = 0, 2
lower_right = 2, 4   # note that I swapped this
lower_left = 4, 2    # and this
n_steps = 3, 3


# build interpolator
m, n = data.shape
x, y = np.arange(m), np.arange(n)

interpolator = RectBivariateSpline(x, y, data)

# build grid
ul,ur,ll,lr = map(np.array, (upper_left,upper_right,lower_left,lower_right))
assert np.allclose(ul + lr, ur + ll)    # make sure edges are parallel

x, y = ul[:, None, None] \
       + np.outer(ll-ul, np.linspace(0.0, 1.0, n_steps[0]))[:, :, None] \
       + np.outer(ur-ul, np.linspace(0.0, 1.0, n_steps[1]))[:, None, :]

# intepolate on grid
print(interpolator.ev(x, y))

Печать:

[[2. 2. 2.]
 [4. 4. 4.]
 [6. 6. 6.]]

Answer 2

Хотя это и отвечает на мой вопрос, я надеюсь, что есть более встроенный способ делать вещи.Я все еще ищу лучший подход.

---

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

skimage, обеспечивающее простую функцию для этой цели..

from skimage.transform import AffineTransform, warp
# Set up grid as in question
transform = AffineTransform(rotation=-math.pi / 4,
                            scale=(math.sqrt(2)/2, math.sqrt(2)/2),
                            translations=(0,2))
grid = warp(grid, transform)

В результате получается

[[2. 2. 2. 2. 2.]
 [3. 3. 3. 3. 3.]
 [4. 4. 4. 4. 4.]
 [5. 5. 5. 5. 5.]
 [6. 6. 6. 6. 6.]]

, который затем может быть просто изменен по желанию.

В общем случае, если у нас есть сетка с размерами xи y, и координаты p1, p2, p3, p4 (начиная с верхнего левого угла, по часовой стрелке), к которому мы хотим повернуть, имеем

rotation = math.atan2(p4.x - p1.x, p4.y - p1.y)
scale = (math.sqrt((p2.y - p1.y) ** 2 + (p2.x - p1.x) ** 2) / x,
         math.sqrt((p4.y - p1.y) ** 2 + (p4.x - p1.x) ** 2) / y)
translation = (0, p1.y)