2012-01-25 10 views
5

Próbuję dopasować obraz 2D w skali Gaussa do obrazu. Hałas jest bardzo niski, więc starałem się obrócić obraz w taki sposób, aby dwie główne osie nie zmieniały się wzajemnie, obliczyć maksimum i obliczyć standardowe odchylenie w obu wymiarach. Broń z wyboru to python.Wyliczenie wektorów własnych obrazu w pytonie

2d more-or-less gaussian distribution

jednak utknąłem na znalezienie wektorów własnych obrazu - numpy.linalg.py zakłada dyskretnych punktów danych. Zastanawiałem się nad tym, czy zdjęcie to jest rozkładem prawdopodobieństwa, próbując kilka tysięcy punktów, a następnie obliczając wektory własne z tego rozkładu, ale jestem pewien, że musi istnieć sposób na znalezienie wektorów własnych (tj. Pół-dur i pół- małe osie elipsy gaussowskiej) bezpośrednio z tego obrazu. Jakieś pomysły?

Thanks a lot :)

Odpowiedz

17

Po krótkiej notatce dostępnych jest kilka narzędzi do dopasowania obrazu gaussowskiego. Jedyne, co mogę sobie wyobrazić, to: scikits.learn, które nie jest całkowicie zorientowane na obraz, ale wiem, że są inne.

Aby obliczyć wektory własne macierzy kowariancji dokładnie tak, jak było to na myśli, jest bardzo kosztowna pod względem obliczeniowym. Musisz powiązać każdy piksel (lub dużą, losową próbkę) obrazu z punktem x, y.

Zasadniczo, można zrobić coś takiego:

import numpy as np 
    # grid is your image data, here... 
    grid = np.random.random((10,10)) 

    nrows, ncols = grid.shape 
    i,j = np.mgrid[:nrows, :ncols] 
    coords = np.vstack((i.reshape(-1), j.reshape(-1), grid.reshape(-1))).T 
    cov = np.cov(coords) 
    eigvals, eigvecs = np.linalg.eigh(cov) 

Można zamiast skorzystania z faktu, że jest to regularnie próbą obraz i obliczyć to chwile (lub „intertial osiach”) zamiast. Będzie to znacznie szybsze w przypadku dużych obrazów.

Jako szybki przykład (używam część jednego z moich previous answers, w przypadku okazać się przydatna ...)

import numpy as np 
import matplotlib.pyplot as plt 

def main(): 
    data = generate_data() 
    xbar, ybar, cov = intertial_axis(data) 

    fig, ax = plt.subplots() 
    ax.imshow(data) 
    plot_bars(xbar, ybar, cov, ax) 
    plt.show() 

def generate_data(): 
    data = np.zeros((200, 200), dtype=np.float) 
    cov = np.array([[200, 100], [100, 200]]) 
    ij = np.random.multivariate_normal((100,100), cov, int(1e5)) 
    for i,j in ij: 
     data[int(i), int(j)] += 1 
    return data 

def raw_moment(data, iord, jord): 
    nrows, ncols = data.shape 
    y, x = np.mgrid[:nrows, :ncols] 
    data = data * x**iord * y**jord 
    return data.sum() 

def intertial_axis(data): 
    """Calculate the x-mean, y-mean, and cov matrix of an image.""" 
    data_sum = data.sum() 
    m10 = raw_moment(data, 1, 0) 
    m01 = raw_moment(data, 0, 1) 
    x_bar = m10/data_sum 
    y_bar = m01/data_sum 
    u11 = (raw_moment(data, 1, 1) - x_bar * m01)/data_sum 
    u20 = (raw_moment(data, 2, 0) - x_bar * m10)/data_sum 
    u02 = (raw_moment(data, 0, 2) - y_bar * m01)/data_sum 
    cov = np.array([[u20, u11], [u11, u02]]) 
    return x_bar, y_bar, cov 

def plot_bars(x_bar, y_bar, cov, ax): 
    """Plot bars with a length of 2 stddev along the principal axes.""" 
    def make_lines(eigvals, eigvecs, mean, i): 
     """Make lines a length of 2 stddev.""" 
     std = np.sqrt(eigvals[i]) 
     vec = 2 * std * eigvecs[:,i]/np.hypot(*eigvecs[:,i]) 
     x, y = np.vstack((mean-vec, mean, mean+vec)).T 
     return x, y 
    mean = np.array([x_bar, y_bar]) 
    eigvals, eigvecs = np.linalg.eigh(cov) 
    ax.plot(*make_lines(eigvals, eigvecs, mean, 0), marker='o', color='white') 
    ax.plot(*make_lines(eigvals, eigvecs, mean, -1), marker='o', color='red') 
    ax.axis('image') 

if __name__ == '__main__': 
    main() 

enter image description here

1

Próbowałaś Analiza głównych składowych (PCA)? Może MDP package może wykonać pracę przy minimalnym wysiłku.

3

Montaż Gaussa solidnie może być zdradliwy. Było zabawa artykuł na ten temat w IEEE Signal Processing Magazine.

Hongwei Guo, "prosty algorytm do montażu Gaussa funkcji" IEEE Signal Processing Magazine, wrzesień 2011, str 134--137

daję implementację 1D przypadku tutaj:

http://scipy-central.org/item/28/2/fitting-a-gaussian-to-noisy-data-points

(Przewiń w dół, aby zobaczyć wynikających pasuje)