Calculer les vecteurs propres de l'image en python

Question

J'essaie d'adapter une gaussienne 2D à une image.Le bruit est très faible, j'ai donc tenté de faire pivoter l'image de telle sorte que les deux axes principaux ne varient pas conjointement, de déterminer le maximum et de calculer simplement l'écart type dans les deux dimensions.L'arme de choix est le python.

Cependant, je suis resté bloqué pour trouver les vecteurs propres de l'image - numpy.linalg.py suppose des points de données discrets.J'ai pensé à prendre cette image comme une distribution de probabilité, à échantillonner quelques milliers de points, puis à calculer les vecteurs propres à partir de cette distribution, mais je suis sûr qu'il doit y avoir un moyen de trouver les vecteurs propres (c'est-à-dire semi-majeurs et semi-majeurs). axes mineurs de l'ellipse gaussienne) directement à partir de cette image.Des idées?

Merci beaucoup :)

Answer 1

Juste un petit mot, il existe plusieurs outils pour adapter une gaussienne à une image.La seule chose à laquelle je peux penser par tête est scikits.apprendre, qui n'est pas complètement orienté image, mais je sais qu'il y en a d'autres.

Calculer les vecteurs propres de la matrice de covariance exactement comme vous l'aviez en tête est très coûteux en termes de calcul.Vous devez associer chaque pixel (ou un échantillon aléatoire de grande taille) de l'image à un point x, y.

En gros, tu fais quelque chose comme :

    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)

Vous pouvez plutôt utiliser le fait qu'il s'agit d'une image régulièrement échantillonnée et calculer ses moments (ou "axes intertiaux") à la place.Ce sera considérablement plus rapide pour les grandes images.

À titre d'exemple rapide, (j'utilise une partie de l'un de mes réponses précédentes, au cas où cela vous serait utile...)

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()

Answer 2

Ajuster une gaussienne de manière robuste peut être délicat.Il y avait un article amusant sur ce sujet dans le magazine IEEE Signal Processing :

Hongwei Guo, "A Simple Algorithme for Adapte A Gaussien Fonction" IEEE Signal Processing Magazine, septembre 2011, pp.134--137

Je donne ici une implémentation du cas 1D :

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

(Faites défiler vers le bas pour voir les ajustements résultants)

Answer 3

Avez-vous essayé l'analyse en composantes principales (ACP) ?Peut-être le Forfait MDP pourrait faire le travail avec un minimum d'effort.