为numpy 3D阵列制作快速定制过滤器

Question

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

我已经调查过scipy.ndimage.generic_过滤器但是我不能像https://ilovesymposia.com/tag/numba中描述的那样编译过滤器，显然是因为windows中的一个numba错误。在

然后，我尝试创建自己的实现，在数组中循环，并希望以后能够jit它。它在没有numba的情况下运行得很好（也很慢），但是jit编译失败了，我无法解码typengerror。在

numpy的meshgrid不受支持，所以它的行为也必须构建（廉价版本）。在

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 filter也会很感激

Answer 1

您可以大大简化代码：

1. 边界上的镜像可以使用^{}
2. 使用SciKit Image的^{}可以有效地完成数据的重叠窗口。在
3. 计算轴上的平均值可以使用np.mean(..., axis=x)来完成，其中x是int或{}的int，表示你想要的轴。在

把这些放在一起，你就得到了一个非常简单的矢量化实现，比如

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