2013-08-25 18 views
6

Ich versuche, einen gleitenden Fensterbetrieb zu vektorisieren.Python - ein gleitendes Fenster vektorisieren

x= vstack((np.array([range(10)]),np.array([range(10)]))) 

x[1,:]=np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]+1],x[1,:]) 

Der n + 1-Wert für jeden aktuellen Wert für Indizes < 5. Aber ich diesen Fehler:

x[1,:]=np.where((x[0,:]<2)&(x[0,:]>0),x[1,x[0,:]+1],x[1,:]) 
IndexError: index (10) out of range (0<=index<9) in dimension 1 

Für die 1-d Fall ein hilfreiches Beispiel entlang der Linien gehen könnte Merkwürdiger würde ich nicht diesen Fehler für den n-1-Wert, die Indizes kleiner als 0 bedeuten würde, es nicht in dem Sinne scheint:

x[1,:]=np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-1],x[1,:]) 

print(x) 

[[0 1 2 3 4 5 6 7 8 9] 
[0 0 1 2 3 5 6 7 8 9]] 

gibt es trotzdem, um dies? ist mein Ansatz völlig falsch? irgendwelche Kommentare würden geschätzt.

EDIT:

Dies ist, was würde ich erreichen möchte ich eine Matrix zu einem numpy Array abflachen, auf dem ich den Mittelwert der 6x6 Nachbarschaft jeder Zelle berechnet werden soll:

matriz = np.array([[1,2,3,4,5], 
    [6,5,4,3,2], 
    [1,1,2,2,3], 
    [3,3,2,2,1], 
    [3,2,1,3,2], 
    [1,2,3,1,2]]) 

# matrix to vector 
vector2 = ndarray.flatten(matriz) 

ncols = int(shape(matriz)[1]) 
nrows = int(shape(matriz)[0]) 

vector = np.zeros(nrows*ncols,dtype='float64') 


# Interior pixels 
if ((i % ncols) != 0 and (i+1) % ncols != 0 and i>ncols and i<ncols*(nrows-1)): 

    vector[i] = np.mean(np.array([vector2[i-ncols-1],vector2[i-ncols],vector2[i-ncols+1],vector2[i-1],vector2[i+1],vector2[i+ncols-1],vector2[i+ncols],vector2[i+ncols+1]])) 
+0

Um zu verdeutlichen, dass Sie "vector2 [i]" nicht in den Mittelwert einschließen wollen oder war das ein Fehler im Code? – Daniel

+0

Ich nicht. Vielen Dank. – JEquihua

+0

Ihr Code berechnet den Durchschnitt einer 3x3-Umgebung jeder Zelle, nicht eine 6x6-Umgebung; war das beabsichtigt? – nneonneo

Antwort

8

Wenn ich das Problem richtig verstehe, möchten Sie den Mittelwert aller Zahlen 1 Schritt um den Index nehmen, den Index vernachlässigen.

ich Ihre Funktion gepatcht haben zu arbeiten, ich glaube, Sie nach so etwas wie diese gingen:

def original(matriz): 

    vector2 = np.ndarray.flatten(matriz) 

    nrows, ncols= matriz.shape 
    vector = np.zeros(nrows*ncols,dtype='float64') 

    # Interior pixels 
    for i in range(vector.shape[0]): 
     if ((i % ncols) != 0 and (i+1) % ncols != 0 and i>ncols and i<ncols*(nrows-1)): 

      vector[i] = np.mean(np.array([vector2[i-ncols-1],vector2[i-ncols],\ 
         vector2[i-ncols+1],vector2[i-1],vector2[i+1],\ 
         vector2[i+ncols-1],vector2[i+ncols],vector2[i+ncols+1]])) 

Ich schrieb dies mit Schneide- und Ansichten mit:

def mean_around(arr): 
    arr=arr.astype(np.float64) 

    out= np.copy(arr[:-2,:-2]) #Top left corner 
    out+= arr[:-2,2:]   #Top right corner 
    out+= arr[:-2,1:-1]   #Top center 
    out+= arr[2:,:-2]   #etc 
    out+= arr[2:,2:] 
    out+= arr[2:,1:-1] 
    out+= arr[1:-1,2:] 
    out+= arr[1:-1,:-2] 

    out/=8.0 #Divide by # of elements to obtain mean 

    cout=np.empty_like(arr) #Create output array 
    cout[1:-1,1:-1]=out  #Fill with out values 
    cout[0,:]=0;cout[-1,:]=0;cout[:,0]=0;cout[:,-1]=0 #Set edges equal to zero 

    return cout 

np.empty_like und anschließend Füllen die Kanten schienen etwas schneller als np.zeros_like. Lassen Sie uns zuerst prüfen, ob sie dasselbe mit IhremArray geben.

print np.allclose(mean_around(matriz),original(matriz)) 
True 

print mean_around(matriz) 
[[ 0.  0.  0.  0.  0. ] 
[ 0.  2.5 2.75 3.125 0. ] 
[ 0.  3.25 2.75 2.375 0. ] 
[ 0.  1.875 2.  2.  0. ] 
[ 0.  2.25 2.25 1.75 0. ] 
[ 0.  0.  0.  0.  0. ]] 

Einige Timings:

a=np.random.rand(500,500) 

print np.allclose(original(a),mean_around(a)) 
True 

%timeit mean_around(a) 
100 loops, best of 3: 4.4 ms per loop 

%timeit original(a) 
1 loops, best of 3: 6.6 s per loop 

Etwa ~ 1500x Speedup.

Sieht aus wie ein guter Ort numba zu verwenden:

def mean_numba(arr): 
    out=np.zeros_like(arr) 
    col,rows=arr.shape 

    for x in xrange(1,col-1): 
     for y in xrange(1,rows-1): 
      out[x,y]=(arr[x-1,y+1]+arr[x-1,y]+arr[x-1,y-1]+arr[x,y+1]+\ 
         arr[x,y-1]+arr[x+1,y+1]+arr[x+1,y]+arr[x+1,y-1])/8. 
    return out 

nmean= autojit(mean_numba) 

können nun gegen alle vorgestellten Methoden vergleichen.

a=np.random.rand(5000,5000) 

%timeit mean_around(a) 
1 loops, best of 3: 729 ms per loop 

%timeit nmean(a) 
10 loops, best of 3: 169 ms per loop 

#CT Zhu's answer 
%timeit it_mean(a) 
1 loops, best of 3: 36.7 s per loop 

#Ali_m's answer 
%timeit fast_local_mean(a,(3,3)) 
1 loops, best of 3: 4.7 s per loop 

#lmjohns3's answer 
%timeit scipy_conv(a) 
1 loops, best of 3: 3.72 s per loop 

A 4-facher Geschwindigkeit mit numba up ist ziemlich nominal anzeigt, dass der numpy Code ungefähr so ​​gut ist wie seine erhalten würde. Ich zog die anderen Codes wie dargestellt, obwohl ich @ CTZhu's Antwort ändern musste, um verschiedene Array-Größen einzuschließen.

+1

Schön. Es ist schneller als meine Version für 'n = 3' um den Faktor zwei, obwohl es für diesen speziellen Fall ziemlich stark abgestimmt ist;). – nneonneo

+0

Ich mag das sehr. Ich bin gerade im Urlaub, aber ich werde das bei meinem speziellen Problem versuchen und zu dir zurückkommen. Ich möchte dies für eine 5000 * 5000-Matrix verwenden und sehen, wie es funktioniert. – JEquihua

+1

@nneonneo 'uniform_filter' war eigentlich die Antwort, die ich in der ersten Iteration dieses Posts verwendet habe, ich bin froh, dass du es vor ein paar Fragen aufbrachtest, es ist unglaublich mächtig und unglaublich schnell. – Daniel

2

Das Problem liegt in x[1,x[0,:]+1], der Index für die 2. Achse: x[0,:]+1 ist [1 2 3 4 5 6 7 8 9 10], in dem Index 10 ist größer als die Dimension von x.

Im Fall x[1,x[0,:]-1] ist der Index der 2. Achse [-1 0 1 2 3 4 5 6 7 8 9] Sie [9 0 1 2 3 4 5 6 7 8] am Ende immer, wie 9 das letzte Element ist und hat einen Index von -1. Der Index des zweiten Elements vom Ende ist -2 und so weiter.

Mit np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-1],x[1,:]) und x[0,:]=[0 1 2 3 4 5 6 7 8 9], was im Wesentlichen los ist, ist, dass die erste Zelle Form genommen wird x[1,:] weil x[0,0] 0 und x[0,:]<5)&(x[0,:]>0 ist False. Die nächsten vier Elemente stammen aus x[1,x[0,:]-1]. Der Rest ist von x[1,:]. Schließlich ist das Ergebnis [0 0 1 2 3 4 5 6 7 8]

Es erscheinen für Gleitfenstersystem von nur 1 Zelle in Ordnung zu sein, aber es wird Sie überraschen mit:

>>> np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-2],x[1,:]) 
array([0, 9, 0, 1, 2, 5, 6, 7, 8, 9]) 

Wenn Sie versuchen, es durch ein Fenster von zwei Zellen zu bewegen .

Für dieses spezifische Problem, wenn wir alles, was in einer Linie halten wollen, dies tun:

>>> for i in [1, 2, 3, 4, 5, 6]: 
    print hstack((np.where(x[1,x[0,:]-i]<x[0, -i], x[1,x[0,:]-i], 0)[:5], x[0,5:])) 

[0 0 1 2 3 5 6 7 8 9] 
[0 0 0 1 2 5 6 7 8 9] 
[0 0 0 0 1 5 6 7 8 9] 
[0 0 0 0 0 5 6 7 8 9] 
[0 0 0 0 0 5 6 7 8 9] 
[0 0 0 0 0 5 6 7 8 9] 

Edit: Jetzt ist Ihre ursprüngliche Frage besser verstehe ich, im Grunde Sie eine 2D nehmen wollen Array und berechnen N * N-Zelle Durchschnitt um jede Zelle. Das ist ziemlich üblich. Zuerst möchten Sie wahrscheinlich N auf ungerade Zahlen beschränken, andernfalls ist ein 2 * 2 Durchschnitt um eine Zelle schwierig zu definieren. Angenommen, wir wollen 3 * 3 Durchschnitt:

#In this example, the shape is (10,10) 
>>> a1=\ 
array([[3, 7, 0, 9, 0, 8, 1, 4, 3, 3], 
    [5, 6, 5, 2, 9, 2, 3, 5, 2, 9], 
    [0, 9, 8, 5, 3, 1, 8, 1, 9, 4], 
    [7, 4, 0, 0, 9, 3, 3, 3, 5, 4], 
    [3, 1, 2, 4, 8, 8, 2, 1, 9, 6], 
    [0, 0, 3, 9, 3, 0, 9, 1, 3, 3], 
    [1, 2, 7, 4, 6, 6, 2, 6, 2, 1], 
    [3, 9, 8, 5, 0, 3, 1, 4, 0, 5], 
    [0, 3, 1, 4, 9, 9, 7, 5, 4, 5], 
    [4, 3, 8, 7, 8, 6, 8, 1, 1, 8]]) 
#move your original array 'a1' around, use range(-2,2) for 5*5 average and so on 
>>> movea1=[a1[np.clip(np.arange(10)+i, 0, 9)][:,np.clip(np.arange(10)+j, 0, 9)] for i, j in itertools.product(*[range(-1,2),]*2)] 
#then just take the average 
>>> averagea1=np.mean(np.array(movea1), axis=0) 
#trim the result array, because the cells among the edges do not have 3*3 average 
>>> averagea1[1:10-1, 1:10-1] 
array([[ 4.77777778, 5.66666667, 4.55555556, 4.33333333, 3.88888889, 
    3.66666667, 4.  , 4.44444444], 
    [ 4.88888889, 4.33333333, 4.55555556, 3.77777778, 4.55555556, 
    3.22222222, 4.33333333, 4.66666667], 
    [ 3.77777778, 3.66666667, 4.33333333, 4.55555556, 5.  , 
    3.33333333, 4.55555556, 4.66666667], 
    [ 2.22222222, 2.55555556, 4.22222222, 4.88888889, 5.  , 
    3.33333333, 4.  , 3.88888889], 
    [ 2.11111111, 3.55555556, 5.11111111, 5.33333333, 4.88888889, 
    3.88888889, 3.88888889, 3.55555556], 
    [ 3.66666667, 5.22222222, 5.  , 4.  , 3.33333333, 
    3.55555556, 3.11111111, 2.77777778], 
    [ 3.77777778, 4.77777778, 4.88888889, 5.11111111, 4.77777778, 
    4.77777778, 3.44444444, 3.55555556], 
    [ 4.33333333, 5.33333333, 5.55555556, 5.66666667, 5.66666667, 
    4.88888889, 3.44444444, 3.66666667]]) 

Ich glaube, Sie brauchen nicht zu glätten Sie 2D-Array, die Verwirrung verursacht. Wenn Sie die Kantenelemente anders handhaben möchten, als sie einfach wegzuschneiden, sollten Sie maskierte Arrays erstellen, indem Sie np.ma im Schritt "Move your original array around" verwenden.

+0

Warum funktioniert es nicht umgekehrt, 10 ist wieder das erste Element? oder wie kann ich machen was ich will? – JEquihua

+0

Oh, im Gegensatz zu Matlab beginnt der Index von Python bei 0. Wenn Sie also "int" verwenden, ist der maximale Index für einen Vektor der Länge 10 9 und wenn Sie x [10] versuchen, erhalten Sie einen 'indexError'. Für 'x = [0 1 2 3 4 5 6 7 8 9]', um 9 zu bekommen, reicht 'x [-1]' oder 'x [9]', aber 'x [10]' wird nicht. –

+0

Ich werde meine Frage bearbeiten, um zu zeigen, was ich wirklich erreichen möchte. Ich wollte einfach keine lange Frage, aber hier geht es. Wie ich denke, du verstehst mich ein bisschen falsch. – JEquihua

4

Es klingt, als ob Sie versuchen, eine 2D-Faltung zu berechnen. Wenn Sie in der Lage sind scipy zu verwenden, würde ich vorschlagen, versuchen scipy.signal.convolve2d:

matriz = np.random.randn(10, 10) 

# to average a 3x3 neighborhood 
kernel = np.ones((3, 3), float) 

# to compute the mean, divide by size of neighborhood 
kernel /= kernel.sum() 

average = scipy.signal.convolve2d(matriz, kernel) 

Der Grund dafür ist der Mittelwert aller 3x3 Nachbarschaften berechnet zu sehen ist, wenn Sie „abrollen“ convolve2d in seine Bestandteile Schleifen. Effektiv (und ignorieren, was an den Rändern der Quelle und Kernel-Arrays geschieht), ist es Computing:

X, Y = kernel.shape 
for i in range(matriz.shape[0]): 
    for j in range(matriz.shape[1]): 
     for ii in range(X): 
      for jj in range(Y): 
       average[i, j] += kernel[ii, jj] * matriz[i+ii, j+jj] 

Also, wenn jeder Wert in Ihrem Kernel 1/(1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1) == 1/9, können Sie den angezeigten Code umschreiben als:

for i in range(matriz.shape[0]): 
    for j in range(matriz.shape[1]): 
     average[i, j] = 1./9 * matriz[i:i+X, j:j+Y].sum() 

Welche genau das gleiche wie die Berechnung der Durchschnitt der Werte in matriz über einen 3x3-Bereich ist, ab i, j.

Ein Vorteil dieser Vorgehensweise besteht darin, dass Sie die Gewichtung Ihrer Umgebung leicht ändern können, indem Sie die Werte in Ihrem Kernel entsprechend einstellen. So zum Beispiel, wenn Sie den Mittelwert in jeder Nachbarschaft doppelt so viel Gewicht wie die anderen geben wollten, könnten Sie Ihren Kernel wie diese bauen:

kernel = np.ones((3, 3), float) 
kernel[1, 1] = 2. 
kernel /= kernel.sum() 

und der Faltungscode würde das gleiche bleiben, aber die Berechnung würde eine andere Art von Durchschnitt (eine "mittengewichtete") ergeben. Hier gibt es viele Möglichkeiten; Hoffentlich bietet dies eine schöne Abstraktion für die Aufgabe, die du machst.

3

Es gibt einfach eine Funktion in der Scipy-Standardbibliothek, die den Mittelwert über gleitende Fenster extrem schnell berechnet. Es heißt uniform_filter. Sie können es verwenden, um Ihre Mean-of-Nachbarschaftsfunktion zu implementieren, wie folgt:

from scipy.ndimage.filters import uniform_filter 
def neighbourhood_average(arr, win=3): 
    sums = uniform_filter(arr, win, mode='constant') * (win*win) 
    return ((sums - arr)/(win*win - 1)) 

Dies gibt einen Array X wo X[i,j] ist der Durchschnitt aller Nachbarn von i,j in arr ohne i,j selbst. Beachten Sie, dass die erste und die letzte Spalte sowie die erste und die letzte Zeile Randbedingungen unterliegen und daher für Ihre Anwendung ungültig sein können (Sie können mode= verwenden, um die Begrenzungsregel bei Bedarf zu steuern).

Da uniform_filter eine hocheffiziente lineare Zeit Algorithmus in geraden C (linear nur in der Größe von arr) implementiert verwendet, sollte es leicht andere Lösungen übertreffen, vor allem, wenn win groß ist.

+0

Sehr interessant.Welche Bedingungen unterliegen den Grenzen? Ich denke, ich möchte die üblichen Bedingungen, aber ich habe das nicht in meiner Frage gepostet. Wie ist das ausgeschlossen (i, j) selbst? Würde es Ihnen etwas ausmachen, den Code ein wenig zu erklären? – JEquihua

+0

'uniform_filter' zentriert das Fenster standardmäßig auf jedes '(i, j)', so dass es z.B. ein 3x3-Fenster "(i-1: i + 2, j-1: j + 2)". Für Werte, die außerhalb des ursprünglichen Arrays liegen, verwendet "uniform_filter" einen Füllwert, der durch "mode" bestimmt wird. Wenn Sie sich nicht um unvollständige Fenster kümmern, können Sie einfach die erste und letzte Zeile und die erste und letzte Spalte löschen oder auf Null setzen. – nneonneo

+1

Es schließt "(i, j)" wegen des '- arr' Bits aus, das den ursprünglichen Wert aus der Fenstersumme entfernt. – nneonneo