为 numpy 3D 数组制作快速自定义过滤器

Question

总的来说，我想制作一个过滤器来计算 3D numpy 数组上的the mean of circular quantities。

我已经查看了 scipy.ndimage.generic_filter，但我无法按照https://ilovesymposia.com/tag/numba 中的描述编译过滤器，显然是由于 Windows 中的 numba 错误。

然后我尝试制作自己的实现来循环遍历数组，并希望之后能够对其进行 jit。它在没有 numba 的情况下运行良好（而且速度很慢），但是 jit 编译失败，我无法解码 TypingError。

numpy 的网格网格不受支持，因此它的行为也必须构建（廉价版本）。

from numba import njit
import numpy as np

@njit
def my_meshgrid(i_, j_,k_):
    #Note: axes 0 and 1 are swapped!
    shape = (len(j_), len(i_), len(k_))
    io = np.empty(shape, dtype=np.int32)
    jo = np.empty(shape, dtype=np.int32)
    ko = np.empty(shape, dtype=np.int32)
    for i in range(len(i_)):
        for j in range(len(j_)):
            for k in range(len(k_)):
                io[j,i,k] = i_[i]
                jo[j,i,k] = j_[j]
                ko[j,i,k] = k_[k]
    return [io,jo, ko]
t3 = my_meshgrid(range(5), range(5,7), range(7,10))
#

@njit
def get_footprint(arr, i , j , k, size=3):
    s = size
    ranges = [range(d-s+1+1,d+s-1) for d in [i,j,k]]
    #Mirror the case where indexes are less than zero
    ind = np.abs(np.meshgrid(*ranges))
    #Mirror the case where indexes are higher than arr.shape:
    for d in range(len(arr.shape)):
        indd = ind[d] - arr.shape[d]
        indd *= -1
        indd = np.abs(indd)
        indd *= -1
        ind[d] = indd
    return arr[ind]


@njit
def mean_angle_filter(degrees, size = 3):
    size = [size]*len(degrees.shape)
    out = np.empty_like(degrees)
    for i in range(degrees.shape[0]):
        for j in range(degrees.shape[1]):
            for k in range(degrees.shape[2]):
                out[i,j,k] = mean_angle(get_footprint(degrees, i,j,k,3))
    return out


@njit
def mean_angle(degrees):
    '''
    https://en.wikipedia.org/wiki/Mean_of_circular_quantities
    '''
    x = np.mean(np.cos(degrees*np.pi/180))
    y = np.mean(np.sin(degrees*np.pi/180))
    return np.arctan2(y,x)*180/np.pi


degrees = np.random.random([20]*3)*90
mean_angle_filter(degrees)

作为 numba 的新手，我很乐意为这个（或类似的）实现获得修复，但在 numpy 中对 mean_angle 过滤器的任何（快速）实现也将不胜感激

Answer 1

您可以大大简化您的代码：

1. 可以使用np.pad() 在边界处进行镜像
2. 使用 SciKit-Image 的 skimage.util.view_as_windows() 可以有效地完成数据的重叠窗口化。
3. 可以使用np.mean(..., axis=x) 计算轴上的平均值，其中x 是int 或tuple 中的tuple，表示您想要表示的轴。

将所有这些放在一起，您将得到一个非常简单的矢量化实现，例如

import numpy as np
import skimage

def mean_angle(degrees, axis=None):
    '''
    https://en.wikipedia.org/wiki/Mean_of_circular_quantities
    '''
    x = np.mean(np.cos(degrees*np.pi/180), axis=axis)
    y = np.mean(np.sin(degrees*np.pi/180), axis=axis)
    return np.arctan2(y,x)*180/np.pi


out = mean_angle(
    skimage.util.view_as_windows(
        np.pad(degrees, (1, 1), mode='symmetric'),
        window_shape=(3, 3, 3)
    ),
    axis=(-1, -2, -3)
)