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

Question

总的来说，我想做一个过滤器来计算 3D numpy 数组上圆形量的平均值。

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

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

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 过滤器的任何(快速(实现也将不胜感激

Answer 1

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

可以使用

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