在尝试使用这里的几个优秀解决方案后,我遇到了包含NaN的数据问题。无论是uniform_filter还是np.cumsum()解决方案都会导致NaN通过输出数组传播而不仅仅是被忽略。
我的解决方案基本上只是将@Jaime答案中的窗口求和函数与卷积交换,这对NaN是稳健的。
def windowed_sum(arr: np.ndarray, radius: int) -> np.ndarray:
"""radius=1 means the pixel itself and the 8 surrounding pixels"""
kernel = np.ones((radius * 2 + 1, radius * 2 + 1), dtype=int)
return convolve(arr, kernel, mode="constant", cval=0.0)
def windowed_var(arr: np.ndarray, radius: int) -> np.ndarray:
"""Note: this returns smaller in size than the input array (by radius)"""
diameter = radius * 2 + 1
win_sum = windowed_sum(arr, radius)[radius:-radius, radius:-radius]
win_sum_2 = windowed_sum(arr * arr, radius)[radius:-radius, radius:-radius]
return (win_sum_2 - win_sum * win_sum / diameter / diameter) / diameter / diameter
def windowed_std(arr: np.ndarray, radius: int) -> np.ndarray:
output = np.full_like(arr, np.nan, dtype=np.float64)
var_arr = windowed_var(arr, radius)
std_arr = np.sqrt(var_arr)
output[radius:-radius, radius:-radius] = std_arr
return output
这个执行速度比uniform_filter
慢一点,但仍然比许多其他方法(如数组堆叠、迭代等)要快得多。
>>> data = np.random.random((1024, 1024))
>>> %timeit windowed_std(data, 4)
158 ms ± 695 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
与执行相同大小数据的uniform_filter相比,该算法仅需约36毫秒
包含一些NaN值:
data = np.arange(100, dtype=np.float64).reshape(10, 10)
data[3:4, 3:4] = np.nan
windowed_std(data, 1)
array([[ nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],
[ nan, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, nan],
[ nan, 8.21, nan, nan, nan, 8.21, 8.21, 8.21, 8.21, nan],
[ nan, 8.21, nan, nan, nan, 8.21, 8.21, 8.21, 8.21, nan],
[ nan, 8.21, nan, nan, nan, 8.21, 8.21, 8.21, 8.21, nan],
[ nan, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, nan],
[ nan, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, nan],
[ nan, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, nan],
[ nan, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, 8.21, nan],
[ nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]])