代码之家 › 专栏 › 技术社区 › JETM

如何插值到旋转网格中?

numpy python-2.7 python

JETM · 技术社区 · 6 年前

[[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

2 回复 | 直到 6 年前

Paul Panzer 6 年前

这可能不足以满足您的口味,但这里有一个方法,使用您的起点(网格角点)更直接,并适用于样条插值(立方每默认值)。

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.]]

user2201041 user2201041 6 年前

虽然这确实回答了我的问题,但我希望有一种更内在的做事方式。我仍在寻求更好的方法。

这实际上是两个问题:旋转原始网格,然后插值。旋转网格,然后将其转换到正确的左上角可以通过 Affine Transformation .

skimage 提供 easy function 为此目的。

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)