2009-11-20 34 views
17

Bir grafik modülü için, bir kısmı rasgele şekiller kümesinin minimum sınırlayıcı elipsini hesaplamak olan bir ödev verdim. Elips eksende hizalanmış olmak zorunda değildir.Sınırlayıcı elips

Bu, AWT şekillerini kullanarak java (euch) içinde çalışıyor, bu nedenle nesnelerin çevrelemesini/kesişimini denetlemek için sağlanan tüm araçları kullanabiliyorum.

+4

Trampet ... ve soru nedir? – ChssPly76

+1

Bu, 3.44'te soruları yazmanın getirdiği şey! Gecenin bu saatinde ödev yaptığımı mı düşünüyorsun ve yarın için bile değil mi? Üniversite bana ne yaptı? ;) – Martin

+0

wow ... sizler havalı şeyler yapıyorsunuz. Belli ki, eksik olanı kaçırmadıkça, bu önemsiz ... – mjv

cevap

33

Minimum Volume Enclosing Ellipsoid ya da 2D kasanızda minimum alanı arıyorsunuz. Bu optimizasyon problemi dışbükeydir ve verimli bir şekilde çözülebilir. İçerdiğim bağlantıdaki MATLAB kodunu gözden geçirin - uygulama önemsizdir ve matris dönüşümünden daha karmaşık bir şey gerektirmez.

Matematikle ilgilenen herkes this document numaralı telefonu okumalıdır.

Ayrıca, elipsin çizilmesi de basittir - bu here bulunabilir, ancak elips üzerinde noktalar oluşturmak için MATLAB-spesifik bir işleve ihtiyacınız olacaktır. Eğer kanonik forma denklemi dönüştürmek nasıl görmek için bu kodu kullanabilirsiniz

matrix form http://mathurl.com/yz7flxe.png

Ama algoritma matris formunda elipsin denklemini verir beri,

Singular Value Decomposition (SVD) adresinden canonical http://mathurl.com/y86tlbw.png

. Ve sonra canonical form kullanarak elips çizmek oldukça kolaydır.

İşte MATLAB kodunun sonucu 10 rastgele 2B nokta (mavi). PCA gibi results

diğer yöntemler elips dış noktaları varyans bir göstergesi olduğu ayrışma (eigen/tekil değer) elde edilen elips minimum sınırlayıcı elips olacağını garanti etmez.

DÜZENLEME:

kimse belgesini okursanız Yani, 2D olarak bu konuda gitmek için iki yol vardır: Burada optimum algoritmanın sözde kod var - Suboptimal algoritma açıkça belgede açıklanmıştır:

Optimal algoritma:

Input: A 2x10 matrix P storing 10 2D points 
     and tolerance = tolerance for error. 
Output: The equation of the ellipse in the matrix form, 
     i.e. a 2x2 matrix A and a 2x1 vector C representing 
     the center of the ellipse. 

// Dimension of the points 
d = 2; 
// Number of points 
N = 10; 

// Add a row of 1s to the 2xN matrix P - so Q is 3xN now. 
Q = [P;ones(1,N)] 

// Initialize 
count = 1; 
err = 1; 
//u is an Nx1 vector where each element is 1/N 
u = (1/N) * ones(N,1)  

// Khachiyan Algorithm 
while err > tolerance 
{ 
    // Matrix multiplication: 
    // diag(u) : if u is a vector, places the elements of u 
    // in the diagonal of an NxN matrix of zeros 
    X = Q*diag(u)*Q'; // Q' - transpose of Q  

    // inv(X) returns the matrix inverse of X 
    // diag(M) when M is a matrix returns the diagonal vector of M 
    M = diag(Q' * inv(X) * Q); // Q' - transpose of Q 

    // Find the value and location of the maximum element in the vector M 
    maximum = max(M); 
    j = find_maximum_value_location(M); 

    // Calculate the step size for the ascent 
    step_size = (maximum - d -1)/((d+1)*(maximum-1)); 

    // Calculate the new_u: 
    // Take the vector u, and multiply all the elements in it by (1-step_size) 
    new_u = (1 - step_size)*u ; 

    // Increment the jth element of new_u by step_size 
    new_u(j) = new_u(j) + step_size; 

    // Store the error by taking finding the square root of the SSD 
    // between new_u and u 
    // The SSD or sum-of-square-differences, takes two vectors 
    // of the same size, creates a new vector by finding the 
    // difference between corresponding elements, squaring 
    // each difference and adding them all together. 

    // So if the vectors were: a = [1 2 3] and b = [5 4 6], then: 
    // SSD = (1-5)^2 + (2-4)^2 + (3-6)^2; 
    // And the norm(a-b) = sqrt(SSD); 
    err = norm(new_u - u); 

    // Increment count and replace u 
    count = count + 1; 
    u = new_u; 
} 

// Put the elements of the vector u into the diagonal of a matrix 
// U with the rest of the elements as 0 
U = diag(u); 

// Compute the A-matrix 
A = (1/d) * inv(P * U * P' - (P * u)*(P*u)'); 

// And the center, 
c = P * u; 
+2

Doğrusal cebirlerde güveniyoruz!Bunu paylaştığın için teşekkür ederim Jacob. Bir şekilde çok daha karmaşık bir çözüm bekliyordum. ama ahh! Algoritma değil cebir düşünüyordum. +1, keşke +2 yapabilseydim, başka bir şey için gidenleri desteklemeliyim "a == b ve a = b" soruları arasındaki fark nedir? Teşekkür ederim. – mjv

+1

Lol, çok teşekkürler! Bu gerçekten büyük bir tesadüf, dün kendi araştırmam için buna bir çözüm buldum! Arkasındaki matematik anlamak oldukça zor ama müthiş kısmı uygulama önemsiz. – Jacob

+1

Bunu açıklayabiliyor musunuz, bu daha çok java gibi, gerçekten matlab'a geldiğinde hiçbir fikrim yok yani burada kayboldum :( – Martin

İlgili konular