我需要一种快速循环遍历 Python 中 Image/Stack 像素的方法
I Need a fast way to loop through pixels of an Image/Stack in Python
我创建了一个 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)
所以在每个像素,它创建一个 window 大小的 kernelxkernelxkernel,找到 window 中像素的中值,并用新的中值替换该像素的值。
我的问题是,它非常慢,我有数千张大图像要处理。必须有一种更快的方法来遍历所有这些像素并仍然能够获得相同的结果。
提前致谢!!
首先,在 python 中循环 3D 矩阵是一个非常非常非常糟糕的主意。为了循环一个大的 3D 矩阵,你最好将 n 维数组降到 Cython or C/C++/Fortran and creating a python extension. However, for this particular case, scipy already contains an implementation of the median filter:
>>> from scipy.ndimage import median_filter
>>> median_filter(my_large_3d_array, radious)
简而言之,在python中没有更快遍历体素的方法(可能numpy iterators would help a bit, but won't increase the performance considerably). If you need to perform more complicated 3D stuff in python, you should consider programming in Cython the loopy interface or, alternatively, using a chunking library such as Dask,它实现了并行 对数组块的操作。
Python 的问题如果 for
循环非常慢,特别是如果它们是嵌套的并且有大数组。因此,没有 标准 pythonic 方法来获得对数组的有效迭代。通常,获得加速的方法是通过向量化操作和 numpy-ticks,但这些都是针对特定问题的,没有通用的 trick,你会学到很多 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
来迭代整个数组,而它是 仅 2ms
对于 Cython版本(x200
加速)。
我创建了一个 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)
所以在每个像素,它创建一个 window 大小的 kernelxkernelxkernel,找到 window 中像素的中值,并用新的中值替换该像素的值。
我的问题是,它非常慢,我有数千张大图像要处理。必须有一种更快的方法来遍历所有这些像素并仍然能够获得相同的结果。
提前致谢!!
首先,在 python 中循环 3D 矩阵是一个非常非常非常糟糕的主意。为了循环一个大的 3D 矩阵,你最好将 n 维数组降到 Cython or C/C++/Fortran and creating a python extension. However, for this particular case, scipy already contains an implementation of the median filter:
>>> from scipy.ndimage import median_filter
>>> median_filter(my_large_3d_array, radious)
简而言之,在python中没有更快遍历体素的方法(可能numpy iterators would help a bit, but won't increase the performance considerably). If you need to perform more complicated 3D stuff in python, you should consider programming in Cython the loopy interface or, alternatively, using a chunking library such as Dask,它实现了并行 对数组块的操作。
Python 的问题如果 for
循环非常慢,特别是如果它们是嵌套的并且有大数组。因此,没有 标准 pythonic 方法来获得对数组的有效迭代。通常,获得加速的方法是通过向量化操作和 numpy-ticks,但这些都是针对特定问题的,没有通用的 trick,你会学到很多 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
来迭代整个数组,而它是 仅 2ms
对于 Cython版本(x200
加速)。