在numpy/scipy中優化實現一個旋轉蒙板

Question

這是我第一次嘗試在numpy中使用步幅，它相對於不同濾鏡的簡單迭代而言的確提高了速度，但它仍然非常慢（並且感覺像是有至少有一兩件事情是完全冗餘或低效的）。

所以我的問題是：有沒有更好的方法來執行此操作或調整我的代碼，使其顯着更快？

該算法對每個像素執行9個不同濾鏡的局部評估，並選擇具有最小標準偏差的一個（我試圖實現Nagau和Matsuyma（1980）「複雜區域照片的結構分析」，如所述在圖像分析書中）。其結果是既平滑和邊緣銳化圖像（相當冷靜，如果你問我！）

import numpy as np 
from scipy import ndimage 
from numpy.lib import stride_tricks 

def get_rotating_kernels(): 

    kernels = list() 

    protokernel = np.arange(9).reshape(3, 3) 

    for k in xrange(9): 

     ax1, ax2 = np.where(protokernel==k) 
     kernel = np.zeros((5,5), dtype=bool) 
     kernel[ax1: ax1+3, ax2: ax2+3] = 1 
     kernels.append(kernel) 

    return kernels 


def get_rotation_smooth(im, **kwargs): 

    kernels = np.array([k.ravel() for k in get_rotating_kernels()], 
       dtype=bool) 

    def rotation_matrix(section): 

     multi_s = stride_tricks.as_strided(section, shape=(9,25), 
      strides=(0, section.itemsize)) 

     rot_filters = multi_s[kernels].reshape(9,9) 

     return rot_filters[rot_filters.std(1).argmin(),:].mean() 

    return ndimage.filters.generic_filter(im, rotation_matrix, size=5, **kwargs) 

from scipy import lena 
im = lena() 
im2 = get_rotation_smooth(im) 

（只是一個註釋，那麼get_rotating_kernel還沒有真正得到優化，因爲幾乎沒有時間花費在那裏反正）

在我的上網本，它花了126s和莉娜畢竟是一個相當小的形象。

編輯：

我建議改變rot_filters.std(1)到rot_filters.var(1)節省相當多的平方根，它在5秒的順序剃掉的東西。

Answer 1

我相信你會有一個艱難的時刻使用Python + scipy優化顯著。但是，我能夠通過使用as_strided直接生成rot_filters而不是通過布爾索引來做出小改進。這是基於一個非常簡單的n維windows函數。 （在我意識到scipy中存在2d卷積函數之前，我寫它來解決this problem。）以下代碼在我的機器上提供了適度的10％加速;請參閱下面的它是如何工作的解釋：

import numpy as np 
from scipy import ndimage 
from numpy.lib import stride_tricks 

# pass in `as_strided` as a default arg to save a global lookup 
def rotation_matrix2(section, _as_strided=stride_tricks.as_strided): 
    section = section.reshape(5, 5) # sqrt(section.size), sqrt(section.size) 
    windows_shape = (3, 3, 3, 3)  # 5 - 3 + 1, 5 - 3 + 1, 3, 3 
    windows_strides = section.strides + section.strides 
    windows = _as_strided(section, windows_shape, windows_strides) 
    rot_filters = windows.reshape(9, 9) 
    return rot_filters[rot_filters.std(1).argmin(),:].mean() 

def get_rotation_smooth(im, _rm=rotation_matrix2, **kwargs): 
    return ndimage.filters.generic_filter(im, _rm, size=5, **kwargs) 

if __name__ == '__main__': 
    import matplotlib.pyplot as plt 
    from scipy.misc import lena 
    im = lena() 
    im2 = get_rotation_smooth(im) 
    #plt.gray()  # Uncomment these lines for 
    #plt.imshow(im2) # demo purposes. 
    #plt.show() 

上述功能rotation_matrix2等同於以下兩個函數（它們共同其實比你原來的功能慢一點，因爲windows是更廣義的）。這與您的原始代碼完全相同 - 將9個3x3窗口創建爲5x5陣列，然後將它們重塑爲9x9陣列進行處理。

def windows(a, w, _as_strided=stride_tricks.as_strided): 
    windows_shape = tuple(sa - sw + 1 for sa, sw in zip(a.shape, w)) 
    windows_shape += w 
    windows_strides = a.strides + a.strides 
    return _as_strided(a, windows_shape, windows_strides) 

def rotation_matrix1(section, _windows=windows): 
    rot_filters = windows(section.reshape(5, 5), (3, 3)).reshape(9, 9) 
    return rot_filters[rot_filters.std(1).argmin(),:].mean() 

windows作品與任何尺寸的陣列，只要在窗口具有相同的維數。下面是它如何工作的故障：

windows_shape = tuple(sa - sw + 1 for sa, sw in zip(a.shape, w)) 

我們可以認爲windows陣列作爲正d陣列的正d陣列。外n-d陣列的形狀由較大陣列內窗口的自由度決定;在每個維度中，窗口可以佔用的位置數量等於較大數組的長度減去窗口長度加一。在這種情況下，我們有一個3x3的窗口放入一個5x5的陣列，所以外面的二維陣列是一個3x3的陣列。

windows_shape += w 

內部n-d數組的形狀與窗口本身的形狀相同。在我們的例子中，這又是一個3x3陣列。

現在的大步。我們必須爲外部n-d數組和內部n-d數組定義步幅。但事實證明它們是一樣的！畢竟，窗口移動通過更大的數組，就像單個索引在數組中移動一樣，對吧？

windows_strides = a.strides + a.strides 

現在，我們有我們需要創建Windows的所有信息：

return _as_strided(a, windows_shape, windows_strides) 

Answer 2

對於複雜的每像素+鄰域操作，您可以考慮使用cython來提高性能。它允許將代碼有效地編寫成for-loops，並且接近於Python語法，後面將其轉換爲C代碼。

爲了尋找靈感，你可以在例如在scikit圖像的代碼來看看：

https://github.com/scikit-image/scikit-image/blob/master/skimage/filter/_denoise.pyx