使用 numpy 旋转 2d 子数组而没有混叠效果

Question

我只想将二维数组中的正值像素围绕中心点旋转一定程度。数据表示来自羽流扩散模型的气溶胶浓度，烟囱位置是旋转的原点。

我想根据风向旋转这个散布模式。

首先计算沿 x 轴的风向的浓度，然后使用围绕我的阵列中心点（烟囱位置）的二维线性旋转将所有点的浓度转换为它们的旋转位置> 0。旋转公式的输入 X、Y 是像素索引。

我的问题是输出是别名的，因为整数变成了浮点数。为了获得整数，我向上或向下舍入了输出。然而，这会产生空单元，随着角度的增加，空单元的数量会越来越多。

谁能帮我找到解决问题的方法？如果可能的话，我想使用 numpy 或最少的软件包来解决这个问题...

我的脚本中处理计算浓度和将像素旋转 50°N 的部分如下。谢谢你的帮助。

def linear2D_rotation(xcoord,ycoord,azimuth_degrees):
    radians = (90 - azimuth_degrees) * (np.pi / 180) # in radians
    xcoord_rotated = (xcoord * np.cos(radians)) - (ycoord * np.sin(radians))
    ycoord_rotated = (xcoord * np.sin(radians)) + (ycoord * np.cos(radians))
    return xcoord_rotated,ycoord_rotated

u_orient = 50 # wind orientation in degres from North
kernel = np.zeros((NpixelY, NpixelX))  # initialize matrix
Yc = int((NpixelY - 1) / 2)  # position of central pixel
Xc = int((NpixelX - 1) / 2)  # position of central pixel

nk = 0
for Y in list(range(0,NpixelX)):
    for X in list(range(0,NpixelY)):
        # compute concentrations only in positive x-direction
        if (X-Xc)>0:
            # nnumber of pixels to origin point (chimney)
            dx = ((X-Xc)+1)
        dy = ((Y-Yc)+1)
        # distance of point to origin (chimney)
        DX = dx*pixel_size_X
        DY = dy*pixel_size_Y
        # compute diffusivity coefficients
        Sy, Sz = calcul_diffusivity_coeff(DX, stability_class)
        # concentration at ground level below the centerline of the plume
        C = (Q / (2 * np.pi * u * Sy * Sz)) * \
            np.exp(-(DY / (2 * Sy)) ** 2) * \
            (np.exp(-((Z - H) / (2 * Sz)) ** 2) + np.exp(-((Z + H) / (2 * Sz)) ** 2))  # at point away from center line
        C = C * 1e9  # convert MBq to Bq

        # rotate only if concentration value at pixel is positive
        if C > 1e-12:
            X_rot, Y_rot = linear2D_rotation(xcoord=dx, ycoord=dy,azimuth_degrees=u_orient)
            X2 = int(round(Xc+X_rot))
            Y2 = int(round(Yc-Y_rot)) # Y increases downwards
            # pixels that fall out of bounds -> ignore
            if (X2 > (NpixelX - 1)) or (X2 < 0) or (Y2 > (NpixelY - 1)):
                continue
            else:
                # replace new pixel position in kernel array
                kernel[Y2, X2] = C

待旋转的原数组

旋转 40°N 的阵列显示数据丢失

Answer 1

您的问题描述不是 100% 清楚，但这里有一些建议：

1.) 不要重新发明轮子。对于像旋转像素这样的事情，有标准的解决方案。使用它们！在这种情况下

scipy.ndimage.affine_transform 执行旋转
指定旋转的齐次坐标矩阵
最近邻插值（下面代码中的参数order=0）。

2.) 不要在不必要的地方循环。通过不处理非正像素获得的速度与循环损失的速度没有任何关系。在手写 python 代码赶上它们之前，编译函数可以绕过很多冗余零。

3.) 不要指望一对一映射像素的解决方案，因为事实是有些点不是最近邻点，有些点与多个其他点最近邻。考虑到这一点，您可能需要考虑更高阶、更平滑的插值。

将您的解决方案与标准工具解决方案进行比较，我们发现后者

更快地给出可比较的结果，并且没有那些孔伪影。

代码（没有绘图）。请注意，我必须转置和 flipud 才能对齐结果：

import numpy as np
from scipy import ndimage as sim
from scipy import stats

def mock_data(n, Theta=50, put_neg=True):
    y, x = np.ogrid[-20:20:1j*n, -9:3:1j*n, ]
    raster = stats.norm.pdf(y)*stats.norm.pdf(x)
    if put_neg:
        y, x = np.ogrid[-5:5:1j*n, -3:9:1j*n, ]
        raster -= stats.norm.pdf(y)*stats.norm.pdf(x)
        raster -= (stats.norm.pdf(y)*stats.norm.pdf(x)).T
    return {'C': raster * 1e-9, 'Theta': Theta}

def rotmat(Theta, offset=None):
    theta = np.radians(Theta)
    c, s = np.cos(theta), np.sin(theta)
    if offset is None:
        return np.array([[c, -s] [s, c]])
    R = np.array([[c, -s, 0], [s, c,0], [0,0,1]])
    to, fro = np.identity(3), np.identity(3)
    offset = np.asanyarray(offset)
    to[:2, 2] = offset
    fro[:2, 2] = -offset
    return to @ R @ fro

def f_pp(C, Theta):
    m, n = C.shape
    clipped = np.maximum(0, 1e9 * data['C'])
    clipped[:, :n//2] = 0
    M = rotmat(Theta, ((m-1)/2, (n-1)/2))
    return sim.affine_transform(clipped, M, order = 0)

def linear2D_rotation(xcoord,ycoord,azimuth_degrees):
    radians = (90 - azimuth_degrees) * (np.pi / 180) # in radians
    xcoord_rotated = (xcoord * np.cos(radians)) - (ycoord * np.sin(radians))
    ycoord_rotated = (xcoord * np.sin(radians)) + (ycoord * np.cos(radians))
    return xcoord_rotated,ycoord_rotated

def f_OP(C, Theta):
    kernel = np.zeros_like(C)
    m, n = C.shape
    for Y in range(m):
        for X in range(n):
            if X > n//2:
                c = C[Y, X] * 1e9
                if c > 1e-12:
                    dx = X - n//2 + 1
                    dy = Y - m//2 + 1
                    X_rot, Y_rot = linear2D_rotation(xcoord=dx, ycoord=dy,azimuth_degrees=Theta)
                    X2 = int(round(n//2+X_rot))
                    Y2 = int(round(m//2-Y_rot)) # Y increases downwards
                    # pixels that fall out of bounds -> ignore
                    if (X2 > (n - 1)) or (X2 < 0) or (Y2 > (m - 1)):
                        continue
                    else:
                        # replace new pixel position in kernel array
                        kernel[Y2, X2] = c
    return kernel

n = 100
data = mock_data(n, 70)

使用 numpy 旋转 2d 子数组而没有混叠效果

Rotating a 2d sub-array using numpy without aliasing effects

numpy

image-processing

rotation

python-2.4

multidimensional-array