Calcular los vectores propios de la imagen en python

Estoy tratando de ajustar un gaussiano 2D a una imagen. El ruido es muy bajo, por lo que mi intento fue rotar la imagen de modo que los dos ejes principales no covariaran, calcular el máximo y calcular la desviación estándar en ambas dimensiones. El arma de elección es Python.Calcular los vectores propios de la imagen en python

2d more-or-less gaussian distribution

Sin embargo, me quedé atrapado en encontrar los vectores propios de la imagen - numpy.linalg.py asume puntos de datos discretos. Pensé en tomar esta imagen como una distribución de probabilidad, probando algunos miles de puntos y luego calculando los autovectores de esa distribución, pero estoy seguro de que debe haber una forma de encontrar los vectores propios (es decir, semi-mayor y semi- ejes menores de la elipse gaussiana) directamente desde esa imagen. ¿Algunas ideas?

Muchas gracias :)

Fuente

2012-01-25 Manuel Ebert

Solo una nota rápida, hay varias herramientas para ajustar un gaussiano a una imagen. Lo único que se me ocurre en la cabeza es scikits.learn, que no está completamente orientado a la imagen, pero sé que hay otros.

Para calcular los vectores propios de la matriz de covarianza exactamente como tenía en mente es muy costoso desde el punto de vista informático. Debe asociar cada píxel (o una muestra aleatoria grande) de una imagen con un punto x, y.

Básicamente, haces algo como:

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)

su lugar, puede hacer uso del hecho de que es una imagen muestreada periódicamente y calcular sus momentos (o "ejes inercial") en su lugar. Esto será considerablemente más rápido para imágenes grandes.

Como un ejemplo rápido, (estoy usando una parte de uno de mis previous answers, en caso de que les sea útil ...)

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