以下是使用 sobel 过滤器避免 for 循环的方法:
import numpy as np
from scipy.ndimage import sobel
def sobHypot_vec(rec):
r = np.abs(rec)
x = sobel(r, 0, mode='constant')
y = sobel(r, 1, mode='constant')
h = np.hypot(x, y)
h = np.apply_over_axes(np.mean, h, [0,1])
return h.argmax()
我不确定 sobel 过滤器在您的应用程序中是否特别必要,如果没有您特定的 20 层“图像”,这很难测试,但您可以尝试使用np.gradient
而不是运行两次 sobel。优点是gradient
在三个维度上运行。您可以忽略第三个组件,只假设前两个组件。这似乎很浪费,但在我的测试中实际上仍然更快。
对于各种随机生成的图像,这给出了与您的代码相同的答案,但要对其进行r = np.random.rand(1024,1024,20) + np.random.rand(1024,1024,20)*1j
测试以确保,并且可能会摆弄dx, dy
np.gradient
def grad_max(rec):
g = np.gradient(np.abs(rec))[:2] # ignore derivative in third dimension
h = np.hypot(*g)
h = np.apply_over_axes(np.mean, h, [0,1]) # mean along first and second dimension
return h.argmax()
使用此代码进行计时:
def sobHypot_clean(rec):
rs = rec.shape
hype = np.ones(rs)
r = np.abs(rec)
for i in xrange(rs[-1]):
ri = r[...,i]
x = sobel(ri, 0, mode='constant')
y = sobel(ri, 1, mode='constant')
hype[...,i] = np.hypot(x,y).mean()
return hype.argmax()
定时:
In [1]: r = np.random.rand(1024,1024,20) + np.random.rand(1024,1024,20)*1j
# Original Post
In [2]: timeit sobHypot(r)
1 loops, best of 3: 9.85 s per loop
#cleaned up a bit:
In [3]: timeit sobHypot_clean(r)
1 loops, best of 3: 7.64 s per loop
# vectorized:
In [4]: timeit sobHypot_vec(r)
1 loops, best of 3: 5.98 s per loop
# using np.gradient:
In [5]: timeit grad_max(r)
1 loops, best of 3: 4.12 s per loop
请在您自己的图像上测试这些函数中的任何一个,以确保它们提供所需的输出,因为不同类型的数组的反应可能与我所做的简单随机测试不同。