python - 我需要一种快速的方法来循环遍历 Python 中图像/堆栈的像素

Question

我创建了一个 3D 中值滤波器，它确实有效，如下所示:

def Median_Filter_3D(image,kernel):

window = np.zeros(shape=(kernel,kernel,kernel), dtype = np.uint8)
n = (kernel-1)/2    #Deals with Image border
imgout = np.empty_like(image)
w,h,l = image.shape()

%%在每个像素上开始循环

for y in np.arange(0,(w-n*2),1):
    for x in np.arange(0,(h-n*2),1):
        for z in np.arange(0,(l-n*2),1):
            window[:,:,:] = image[x:x+kernel,y:y+kernel,z:z+kernel]
            med = np.median(window)
            imgout[x+n,y+n,z+n] = med 
return(imgout)

因此在每个像素处，它创建一个大小为 kernelxkernelxkernel 的窗口，找到窗口中像素的中值，并用新的中值替换该像素的值。

我的问题是，它非常慢，我有数千张大图像要处理。必须有一种更快的方法来遍历所有这些像素并仍然能够获得相同的结果。

提前致谢!

Answer 1

首先，在 python 中循环 3D 矩阵是一个非常非常非常糟糕的主意。为了循环一个大的 3D 矩阵，你最好转到 Cython。或 C/C++/Fortran 并创建 python 扩展。但是，对于这种特殊情况，scipy 已经包含了 median filter 的实现。对于 n 维数组:

>>> from scipy.ndimage import median_filter
>>> median_filter(my_large_3d_array, radious)

简而言之，在 python 中没有更快遍历体素的方法(也许 numpy iterators 会有所帮助，但不会显着提高性能)。如果你需要在 python 中执行更复杂的 3D 东西，你应该考虑在 Cython 中编程 loopy 接口(interface)，或者使用分块库，如 Dask ，它为数组 block 实现并行操作。

---

Python 的问题是 for 循环非常慢，特别是当它们是嵌套的并且有大数组时。因此，没有标准的 pythonic 方法来获得对数组的有效迭代。通常，获得加速的方法是通过矢量化操作和 numpy-ticks，但这些都是针对特定问题的，没有通用的技巧，您将在 SO 中学到很多 numpy 技巧.

作为一种通用方法，如果您确实需要遍历数组，您可以在 Cython 中编写代码。 . Cython 是 Python 的类 C 扩展。您使用 Python 语法编写代码，但指定变量类型(如在 C 中，使用 int 或 float)。然后该代码会自动编译为 C 并可以从 python 调用。一个简单的例子:

示例 Python 循环函数:

import numpy as np

def iter_A(A):
    B = np.empty(A.shape, dtype=np.float64)

    for i in range(A.shape[0]):
        for j in range(A.shape[1]):
            B[i, j] = A[i, j] * 2
    return B

我知道上面的代码有点多余，可以写成B = A * 2，但它的目的只是为了说明 python 循环非常慢。

函数的 Cython 版本:

import numpy as np
cimport numpy as np

def iter_A_cy(double[:, ::1] A):
    cdef Py_ssize_t H = A.shape[0], W = A.shape[1]
    cdef double[:, ::1] B = np.empty((H, W), dtype=np.float64)
    cdef Py_ssize_t i, j

    for i in range(H):
        for j in range(W):
            B[i, j] = A[i, j] * 2

    return np.asarray(B)

两种实现的测试速度:

>>> import numpy as np
>>> A = np.random.randn(1000, 1000)
>>> %timeit iter_A(A)
1 loop, best of 3: 399 ms per loop
>>> %timeit iter_A_cy(A)
100 loops, best of 3: 2.11 ms per loop

注意:您不能按原样运行 Cython 函数。你需要把它放在一个单独的文件中，然后先编译它(或者在 IPython Notebook 中使用 %%cython magic)。

它表明原始 python 版本用 400ms 迭代整个数组，而 Cython 版本仅 2ms(x200 加速)。