2012-01-25 19 views
5

私は2Dガウスを画像にフィットさせようとしています。騒音は非常に低いので、2つの主軸が共に変化しないように画像を回転させ、最大値を計算して両方の次元の標準偏差を計算するだけでした。選択の武器は、Pythonです。Pythonで画像の固有ベクトルを計算する

2d more-or-less gaussian distribution

しかし、私は、画像の固有ベクトルを見つけることで行き詰まって - numpy.linalg.pyは、離散的なデータポイントを想定しています。私はこの画像を確率分布にして数千点をサンプリングし、その分布から固有ベクトルを計算することを考えましたが、固有ベクトル(すなわち、半主成分および半主成分)を見つける方法がなければならないと確信していますガウシアン楕円の短軸)をその画像から直接得る。何か案は?

どうもありがとう:)

答えて

17

画像にガウス分布を合わせるためのいくつかのツールがあります。私が頭の上から考えることができるのは、完全に画像指向ではないscikits.learnですが、他にもあります。

共分散行列の固有ベクトルを正確に計算するには、計算コストが非常に高くなります。画像の各ピクセル(または大規模なランダムサンプル)をx、yポイントに関連付ける必要があります。代わりに、あなたが代わりにそれを定期的にサンプリングされた画像だという事実を利用すると、それは瞬間(または「intertial軸」)だ計算することができ

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) 

基本的に、あなたのような何かを行います。これは、大きな画像ではかなり高速になります。簡単な例として、

、(私は私のprevious answersの1の一部を使用しています、場合にあなたはそれが便利...)ガウシアンフィッティング

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

あなたは(PCA)主成分分析をしてみてくださいましたか?たぶん、MDP packageは最小限の労力で仕事をすることができます。

3

堅牢にすることができトリッキー。 IEEE信号処理マガジンのこのトピックに関する楽しい記事がありました:。

Hongwei郭、IEEE 信号処理マガジン、2011年9月、頁「ガウス関数を装着するための単純なアルゴリズムは、」134--137

私はここ1Dケースの実装を与える:

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

(結果の適合を確認するためにスクロールダウン)

関連する問題