2015-04-19 20 views
7
def GaussianMatrix(X,sigma): 
    row,col=X.shape 
    GassMatrix=np.zeros(shape=(row,row)) 
    X=np.asarray(X) 
    i=0 
    for v_i in X: 
     j=0 
     for v_j in X: 
      GassMatrix[i,j]=Gaussian(v_i.T,v_j.T,sigma) 
      j+=1 
     i+=1 
    return GassMatrix 
def Gaussian(x,z,sigma): 
    return np.exp((-(np.linalg.norm(x-z)**2))/(2*sigma**2)) 

Questo è il mio modo attuale. C'è un modo in cui posso usare l'operazione matrix per fare questo? X è i punti dati.Come calcolare una matrice di kernel gaussiana in modo efficiente in numpy?

risposta

17

Si desidera utilizzare il kernel gaussiano per es. levigatura dell'immagine? Se è così, c'è una funzione gaussian_filter() in SciPy:

In alternativa, questo dovrebbe funzionare:

import numpy as np 
import scipy.stats as st 

def gkern(kernlen=21, nsig=3): 
    """Returns a 2D Gaussian kernel array.""" 

    interval = (2*nsig+1.)/(kernlen) 
    x = np.linspace(-nsig-interval/2., nsig+interval/2., kernlen+1) 
    kern1d = np.diff(st.norm.cdf(x)) 
    kernel_raw = np.sqrt(np.outer(kern1d, kern1d)) 
    kernel = kernel_raw/kernel_raw.sum() 
    return kernel 

ingresso:

import matplotlib.pyplot as plt 
plt.imshow(gkern(21), interpolation='none') 

uscita: enter image description here

+1

Perché si prende la radice quadrata del prodotto esterno (cioè 'kernel_raw = np.sqrt (np.outer (kern1d, kern1d)) ') e non solo moltiplicarli? Mi sento come se mi mancasse qualcosa qui .. – trueter

+0

potresti fornire alcuni dettagli, per favore, su come funziona la tua funzione? Perché hai bisogno di 'np.diff (st.norm.cdf (x))'? –

+0

inoltre, l'implementazione fornisce risultati diversi da quelli di altri nella pagina :( –

3

linalg.norm prende un parametro axis. Con un po 'di sperimentazione Imparai a calcolare la norma per tutte le combinazioni di righe con

np.linalg.norm(x[None,:,:]-x[:,None,:],axis=2) 

Si espande x in una matrice 3D di tutte le differenze, e prende la norma sull'ultima dimensione.

Quindi posso applicare questo al codice aggiungendo il parametro axis al Gaussian:

def Gaussian(x,z,sigma,axis=None): 
    return np.exp((-(np.linalg.norm(x-z, axis=axis)**2))/(2*sigma**2)) 

x=np.arange(12).reshape(3,4) 
GaussianMatrix(x,1) 

produce

array([[ 1.00000000e+00, 1.26641655e-14, 2.57220937e-56], 
     [ 1.26641655e-14, 1.00000000e+00, 1.26641655e-14], 
     [ 2.57220937e-56, 1.26641655e-14, 1.00000000e+00]]) 

Abbinamento:

Gaussian(x[None,:,:],x[:,None,:],1,axis=2) 

array([[ 1.00000000e+00, 1.26641655e-14, 2.57220937e-56], 
     [ 1.26641655e-14, 1.00000000e+00, 1.26641655e-14], 
     [ 2.57220937e-56, 1.26641655e-14, 1.00000000e+00]]) 
+0

Come si specifica il valore per Sigma? – Rojin

15

si può semplicemente gaussian- filtrare un semplice 2D dirac function, il risultato è quindi il filtro f unzione che veniva utilizzato:

import numpy as np 
import scipy.ndimage.filters as fi 

def gkern2(kernlen=21, nsig=3): 
    """Returns a 2D Gaussian kernel array.""" 

    # create nxn zeros 
    inp = np.zeros((kernlen, kernlen)) 
    # set element at the middle to one, a dirac delta 
    inp[kernlen//2, kernlen//2] = 1 
    # gaussian-smooth the dirac, resulting in a gaussian filter mask 
    return fi.gaussian_filter(inp, nsig) 
+0

Grande, esattamente quello di cui avevo bisogno, grazie mille! – Minato

7

Io stesso usato la risposta accettata per la mia elaborazione delle immagini, ma lo trovo (e le altre risposte) troppo dipendente da altri moduli. Inoltre, la risposta accettata a volte produce kernel con un numero di voci pari a zero alla fine.

Pertanto, ecco la mia soluzione compatta:

import numpy as np 


def gkern(l=5, sig=1.): 
    """ 
    creates gaussian kernel with side length l and a sigma of sig 
    """ 

    ax = np.arange(-l // 2 + 1., l // 2 + 1.) 
    xx, yy = np.meshgrid(ax, ax) 

    kernel = np.exp(-(xx**2 + yy**2)/(2. * sig**2)) 

    return kernel/np.sum(kernel) 
2

Un 2D matrice kernel gaussiano può essere calcolata con la radiodiffusione NumPy,

def gaussian_kernel(size=21, sigma=3): 
    """Returns a 2D Gaussian kernel. 
    Parameters 
    ---------- 
    size : float, the kernel size (will be square) 

    sigma : float, the sigma Gaussian parameter 

    Returns 
    ------- 
    out : array, shape = (size, size) 
     an array with the centered gaussian kernel 
    """ 
    x = np.linspace(- (size // 2), size // 2) 
    x /= np.sqrt(2)*sigma 
    x2 = x**2 
    kernel = np.exp(- x2[:, None] - x2[None, :]) 
    return kernel/kernel.sum() 

Per le piccole dimensioni del kernel questo dovrebbe essere ragionevolmente veloce.

Nota: ciò semplifica la modifica del parametro sigma rispetto alla risposta accettata.

+0

Penso che intendessi 'np.linspace (- (size // 2), size // 2) 'Altrimenti, l'intervallo è un po 'più lungo sul lato sinistro di zero (perché' (-21) // 2 = -11', considerando '21 // 2 = 10'. –

+0

@ CiprianTomoiagă Grazie, risolto. – rth

2

Sto cercando di migliorare su FuzzyDuck's answer qui. Penso che questo approccio sia più breve e più facile da capire. Qui sto usando signal.scipy.gaussian per ottenere il kernel gaussiano 2D.

import numpy as np 
from scipy import signal 

def gkern(kernlen=21, std=3): 
    """Returns a 2D Gaussian kernel array.""" 
    gkern1d = signal.gaussian(kernlen, std=std).reshape(kernlen, 1) 
    gkern2d = np.outer(gkern1d, gkern1d) 
    return gkern2d 

Tracciando utilizzando matplotlib.pyplot:

import matplotlib.pyplot as plt 
plt.imshow(gkern(21), interpolation='none') 

Gaussian kernel plotted using matplotlib

2

Costruire sulla risposta di Teddy Hartanto. Puoi semplicemente calcolare le tue funzioni Gaussiane unidimensionali e quindi usare np.outer per calcolare quello bidimensionale. Modo molto veloce ed efficiente.

Con il codice qui sotto è anche possibile utilizzare diversi Sigma per ogni dimensione

import numpy as np 
def generate_gaussian_mask(shape, sigma, sigma_y=None): 
    if sigma_y==None: 
     sigma_y=sigma 
    rows, cols = shape 

    def get_gaussian_fct(size, sigma): 
     fct_gaus_x = np.linspace(0,size,size) 
     fct_gaus_x = fct_gaus_x-size/2 
     fct_gaus_x = fct_gaus_x**2 
     fct_gaus_x = fct_gaus_x/(2*sigma**2) 
     fct_gaus_x = np.exp(-fct_gaus_x) 
     return fct_gaus_x 

    mask = np.outer(get_gaussian_fct(rows,sigma), get_gaussian_fct(cols,sigma_y)) 
    return mask 
Problemi correlati