Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Compute eigenvectors of image in python

I'm trying to fit a 2D Gaussian to an image. Noise is very low, so my attempt was to rotate the image such that the two principal axes do not co-vary, figure out the maximum and just compute the standard deviation in both dimensions. Weapon of choice is python.

2d more-or-less gaussian distribution

However I got stuck at finding the eigenvectors of the image - numpy.linalg.py assumes discrete data points. I thought about taking this image to be a probability distribution, sampling a few thousand points and then computing the eigenvectors from that distribution, but I'm sure there must be a way of finding the eigenvectors (ie., semi-major and semi-minor axes of the gaussian ellipse) directly from that image. Any ideas?

Thanks a lot :)

like image 202
Manuel Ebert Avatar asked Jan 25 '12 15:01

Manuel Ebert


People also ask

How eigenvalues and eigenvectors are used in image processing?

An eigenvalue/eigenvector decomposition of the covariance matrix reveals the principal directions of variation between images in the collection. This has applications in image coding, image classification, object recognition, and more.

How do you find eigenvectors and eigenvalues in Numpy?

In NumPy we can compute the eigenvalues and right eigenvectors of a given square array with the help of numpy. linalg. eig(). It will take a square array as a parameter and it will return two values first one is eigenvalues of the array and second is the right eigenvectors of a given square array.


1 Answers

Just a quick note, there are several tools to fit a gaussian to an image. The only thing I can think of off the top of my head is scikits.learn, which isn't completely image-oriented, but I know there are others.

To calculate the eigenvectors of the covariance matrix exactly as you had in mind is very computationally expensive. You have to associate each pixel (or a large-ish random sample) of image with an x,y point.

Basically, you do something like:

    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)

You can instead make use of the fact that it's a regularly-sampled image and compute it's moments (or "intertial axes") instead. This will be considerably faster for large images.

As a quick example, (I'm using a part of one of my previous answers, in case you find it useful...)

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

like image 77
Joe Kington Avatar answered Sep 28 '22 04:09

Joe Kington