2012-01-25 19 views
5

में छवि के ईगेंवेक्टर की गणना करें मैं एक छवि में 2 डी गाऊशियन फिट करने की कोशिश कर रहा हूं। शोर बहुत कम है, इसलिए मेरा प्रयास छवि को घुमाने के लिए था कि दो प्रमुख अक्षों में सह-भिन्नता नहीं है, अधिकतम पता लगाएं और दोनों आयामों में मानक विचलन की गणना करें। पसंद का हथियार अजगर है।पाइथन

2d more-or-less gaussian distribution

हालांकि मैं छवि के eigenvectors ढूँढने में अटक गया - numpy.linalg.py असतत डेटा बिंदुओं मान लिया गया है। मैंने इस छवि को एक संभाव्यता वितरण के रूप में लेने के बारे में सोचा, कुछ हज़ार अंक नमूना और फिर उस वितरण से ईजिनवेक्टरों की गणना करना, लेकिन मुझे यकीन है कि ईगेंवेक्टर (यानी, अर्ध-प्रमुख और अर्द्ध- गॉसियन एलीपसे के नाबालिग अक्ष) सीधे उस छवि से। कोई विचार?

धन्यवाद एक बहुत :)

उत्तर

17

बस एक त्वरित नोट, एक छवि के लिए एक गाऊशियन फिट करने के लिए कई उपकरण हैं। मेरे सिर के शीर्ष से बाहर निकलने वाली एकमात्र चीज scikits.learn है, जो पूरी तरह से छवि उन्मुख नहीं है, लेकिन मुझे पता है कि अन्य भी हैं।

कॉन्वर्सिस मैट्रिक्स के ईजिनवेक्टरों की गणना करने के लिए बिल्कुल आपके मन में बहुत ही कम्प्यूटेशनल महंगा है। आपको x, y बिंदु के साथ छवि के प्रत्येक पिक्सेल (या एक बड़े-आश यादृच्छिक नमूना) को जोड़ना होगा।

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) 

आप के बजाय तथ्य यह है कि यह एक नियमित रूप से नमूना छवि है और यह बजाय क्षणों (या "intertial कुल्हाड़ियों") है की गणना का उपयोग कर सकते:

मूल रूप से, आप की तरह कुछ है। बड़ी छवियों के लिए यह काफी तेज होगा।

एक त्वरित उदाहरण के रूप में, (मैं अपने previous answers में से एक का एक हिस्सा उपयोग कर रहा हूँ, मामले में आप इसे उपयोगी पाते ...)

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

आप की कोशिश की प्रधानाचार्य घटक विश्लेषण (पीसीए)? शायद MDP package कम से कम प्रयास के साथ काम कर सकता है।

3

एक गाऊसी फिटिंग मजबूती के साथ किया जा सकता है मुश्किल। आईईईई सिग्नल प्रोसेसिंग पत्रिका में इस विषय पर एक मजेदार लेख थी:।

Hongwei गुओ, आईईईई सिग्नल प्रोसेसिंग पत्रिका, सितंबर 2011, पीपी "एक गाऊसी समारोह फिटिंग लिए एक सरल एल्गोरिथ्म" 134--137

मैं -1 डी मामले के एक कार्यान्वयन यहाँ दे:

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

(जिसके परिणामस्वरूप फिट देखने के लिए नीचे स्क्रॉल करें)

+०१२३५१६४१०६१

 संबंधित मुद्दे

  • कोई संबंधित समस्या नहीं^_^