使用高效移動平均濾波器的步幅

Question

我最近在answer to this post中瞭解到strides，並且想知道如何使用它們來計算比我在this post（使用卷積濾波器）中提出的更高效的移動平均濾波器。

這是我到目前爲止。它採用原始數組的視圖，然後將其滾動必要的數量，然後求和內核值以計算平均值。我意識到邊緣處理不正確，但我可以在這之後照顧...有更好更快的方法嗎？其目標是過濾大型浮點陣列，最大尺寸爲5000x5000 x 16層，這是一項任務，其工作相當緩慢。

請注意，我正在尋找8鄰居連接，即3x3濾鏡取9個像素（焦點像素周圍8個）的平均值，並將該值分配給新圖像中的像素。我如何看到該工作

import numpy, scipy 

filtsize = 3 
a = numpy.arange(100).reshape((10,10)) 
b = numpy.lib.stride_tricks.as_strided(a, shape=(a.size,filtsize), strides=(a.itemsize, a.itemsize)) 
for i in range(0, filtsize-1): 
    if i > 0: 
     b += numpy.roll(b, -(pow(filtsize,2)+1)*i, 0) 
filtered = (numpy.sum(b, 1)/pow(filtsize,2)).reshape((a.shape[0],a.shape[1])) 
scipy.misc.imsave("average.jpg", filtered) 

EDIT澄清：

當前代碼：

1. 使用stride_tricks以產生，例如[[0,1,2陣列]，[1,2， 3]，[2,3,4] ...]，它對應於濾波器內核的第一行。
2. 沿着縱軸滾動以獲取內核的中間行[[10,11,12]，[11,12,13]，[13,14,15] ...]並將其添加到數組中我得到了1）
3. 重複以獲得內核[[20,21,22]，[21,22,23]，[22,23,24] ...]的底行。此時，我將每行的總和除以濾鏡中元素的數量，給出每個像素的平均值（移動1行和1列，邊緣有一些奇特點，但我可以以後照顧）。

我期望的是更好地使用stride_tricks來直接獲取9個值或內核元素的總和，對於整個數組，或者有人可以讓我相信另一個更有效的方法。

Answer 1

爲了什麼是值得的，以下是你如何使用「花哨」的大步技巧做到這一點。我昨天會發布這個消息，但是被實際的工作分心了！ :)

@Paul & @eat都有很好的實現使用各種其他方式做到這一點。爲了繼續前面的問題，我想我會發布N維等價物。

但是，您不會顯着擊敗> 1D陣列的scipy.ndimage函數。 （scipy.ndimage.uniform_filter應該擊敗scipy.ndimage.convolve，雖然）

此外，如果您試圖獲得一個多維的移動窗口，您可能會有無意中使您的數組的副本炸燬內存的風險。儘管最初的「滾動」數組只是您原始數組內存的視圖，但複製數組的任何中間步驟都會使您的原始數組的大小大於（即讓我們假設您使用100x100原始陣列...查看（對於（3,3）的過濾器大小）將是98x98x3x3，但使用與原始存儲器相同的內存。但是，任何副本將使用的內存量爲全 98x98x3x3陣列會!!）

基本上，利用瘋狂的跨步技巧是偉大的，當你想在矢量化的ndarray的單軸移動窗口操作。它使得很容易計算諸如移動標準偏差等事情，而且開銷很小。當你想開始沿着多個軸來做這件事的時候，這是可能的，但是你通常會用更專業的功能更好。 （如scipy.ndimage等）

無論如何，這裏是你如何做到這一點：

import numpy as np 

def rolling_window_lastaxis(a, window): 
    """Directly taken from Erik Rigtorp's post to numpy-discussion. 
    <http://www.mail-archive.com/numpy-discussion@scipy.org/msg29450.html>""" 
    if window < 1: 
     raise ValueError, "`window` must be at least 1." 
    if window > a.shape[-1]: 
     raise ValueError, "`window` is too long." 
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) 
    strides = a.strides + (a.strides[-1],) 
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides) 

def rolling_window(a, window): 
    if not hasattr(window, '__iter__'): 
     return rolling_window_lastaxis(a, window) 
    for i, win in enumerate(window): 
     if win > 1: 
      a = a.swapaxes(i, -1) 
      a = rolling_window_lastaxis(a, win) 
      a = a.swapaxes(-2, i) 
    return a 

filtsize = (3, 3) 
a = np.zeros((10,10), dtype=np.float) 
a[5:7,5] = 1 

b = rolling_window(a, filtsize) 
blurred = b.mean(axis=-1).mean(axis=-1) 

所以，我們得到當我們做b = rolling_window(a, filtsize)是一個8x8x3x3陣列，這實際上是一個視圖到相同的內存原來的10x10陣列。我們可以很容易地沿着不同的軸使用不同的過濾器尺寸，或者只沿着N維陣列的選定軸進行操作（即，在4維陣列上的filtsize = (0,3,0,3)會給我們6維視圖）。

然後，我們可以將一個任意函數重複應用到最後一個軸上，以有效地計算移動窗口中的事物。

但是，因爲我們在mean（或std或其他）的每一步中都存儲了比我們原始數組大得多的臨時數組，所以這不是高效的存儲器！它的速度也不會很快。

爲ndimage等效就是：

blurred = scipy.ndimage.uniform_filter(a, filtsize, output=a) 

這將處理各種邊界條件，做到就地「模糊」，而不需要數組的臨時副本，並非常快。大步技巧是一個函數一起一個軸適用於移動窗口的好方法，但他們沒有這樣做沿多個軸線的好方法，通常....

只是我的$ 0.02時，在任何速度...

Answer 2

我有信心需要解決的一件事是您的視圖陣列b。

它有來自未分配的內存的一些項目，所以你會得到崩潰。

鑑於新的算法的描述，需要固定的第一件事是事實，你邁進的a分配外：

bshape = (a.size-filtsize+1, filtsize) 
bstrides = (a.itemsize, a.itemsize) 
b = numpy.lib.stride_tricks.as_strided(a, shape=bshape, strides=bstrides) 

---

更新

因爲我仍然不太清楚方法，似乎有更簡單的方法來解決這個問題，我只是想把它放在這裏：

A = numpy.arange(100).reshape((10,10)) 

shifts = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)] 
B = A[1:-1, 1:-1].copy() 
for dx,dy in shifts: 
    xstop = -1+dx or None 
    ystop = -1+dy or None 
    B += A[1+dx:xstop, 1+dy:ystop] 
B /= 9 

......這看起來很直接。唯一無關的操作是它只分配並填充B一次。無論如何，所有的添加，分割和索引都必須完成。如果你正在做16個樂隊，如果你的意圖是保存一個圖像，你仍然只需要分配B一次。即使這沒有幫助，它可能會澄清爲什麼我不瞭解問題，或者至少作爲衡量其他方法加速的基準。這將運行在2.6秒在我的筆記本float64的，其中0.5是B

Answer 3

創立了5K X 5K陣列上讓我們來看看：

這不是很清楚的形式你的問題，但是現在我假設你會喜歡顯着提高這種平均。

import numpy as np 
from numpy.lib import stride_tricks as st 

def mf(A, k_shape= (3, 3)): 
    m= A.shape[0]- 2 
    n= A.shape[1]- 2 
    strides= A.strides+ A.strides 
    new_shape= (m, n, k_shape[0], k_shape[1]) 
    A= st.as_strided(A, shape= new_shape, strides= strides) 
    return np.sum(np.sum(A, -1), -1)/ np.prod(k_shape) 

if __name__ == '__main__': 
    A= np.arange(100).reshape((10, 10)) 
    print mf(A) 

現在，你會期待什麼樣的性能改進？

更新：
首先，一個警告：在它的代碼的當前狀態不正常適應「內核」的形狀。然而，這不是我現在主要關心的問題（無論如何，這個想法已經存在，如何正確適應）。

我剛纔選擇了4D A的新造型直觀，這對我來說真的有意義思考一個2D「內核」中心爲中心，以原始2D A.

可是對每個網格位置4D整形可能實際上並不是最好的。我認爲這裏真正的問題是總結的表現。人們應該能夠找到（4D A的）「最佳順序」，以充分利用你的機器緩存架構。然而，對於那些與你的機器緩存「協作」的「小」陣列和那些不支持（至少不那麼直截了當）的大陣列，這個順序可能不一樣。

更新2：
這裏是mf略加修改。顯然，最好先重塑一個3D數組，然後不要總結點積（這樣做的好處是，內核可以是任意的）。然而，它仍然比Pauls更新功能慢3倍（在我的機器上）。

def mf(A): 
    k_shape= (3, 3) 
    k= np.prod(k_shape) 
    m= A.shape[0]- 2 
    n= A.shape[1]- 2 
    strides= A.strides* 2 
    new_shape= (m, n)+ k_shape 
    A= st.as_strided(A, shape= new_shape, strides= strides) 
    w= np.ones(k)/ k 
    return np.dot(A.reshape((m, n, -1)), w) 

Answer 4

我還沒有與Python熟悉，寫出來的代碼，而是要加快回旋兩個最好的方法是要麼單獨的過濾器或使用傅立葉變換。

分離濾波器：卷積是O（M * N），其中M和N分別是圖像和濾波器中的像素數。由於使用3×3內核進行平均濾波相當於先用3×1內核和1×3內核進行濾波，則可以通過連續卷積來獲得(3+3)/(3*3) =〜30％的速度提升，其中兩個1 -d內核（隨着內核變大，這顯然會變得更好）。當然，您仍然可以在這裏使用步幅技巧。

傅立葉變換：conv(A,B)相當於ifft(fft(A)*fft(B))，即在直接空間卷積變成在傅立葉空間中，在那裏A是你的圖像和B是您的濾波器相乘。由於傅里葉變換的（元素方式）乘法要求A和B的大小相同，所以B是一個size(A)的數組，其內核位於圖像的中心位置，其他位置爲零。要在陣列的中心放置一個3x3的內核，您可能需要將A填充爲奇數大小。根據傅立葉變換的實現，這可能比卷積快得多（並且如果多次應用相同的濾波器，則可以預先計算fft(B)，節省另外30％的計算時間）。