2014-03-03 18 views
5

Ich habe versucht, ein Problem zu lösen. Ich bin überrascht, dass ich im Netz nichts wirklich Nützliches finden konnte.Kovarianz Matrix einer Ellipse

Ich weiß, dass aus den Eigenwerten der Kovarianzmatrix der Ellipse die Haupt- und die Nebenachse der Ellipse berechnet werden können. Wie die folgenden:

a1 = 2*sqrt(e1) 
a2 = 2*sqrt(e2) 

wo a1 und a2 die Haupt- und Nebenachse sind, und e1e2 sind die Eigenwerte der Kovarianzmatrix.

Meine Frage ist: gegebenen Randpunkte (xi,yi) der Bildellipse, wie es möglich ist, die 2 × 2 Kovarianzmatrix dieser Ellipse zu finden?

+0

Sollte die Matrix nicht einfach die Kovarianz aller 'xi'-s' yi'-s sein? – Shai

+0

Ich bin mir nicht sicher! Ich erzeuge Kantenpunkte für einen Kreis mit dem Radius 100. Dann habe ich ein 'p = [xi, yi]' definiert, wobei P eine Kantenpunkte-Matrix 'n x 2' ist. Ich habe den Matlab-Befehl 'cov (P)' verwendet. Ich habe den Radius des Kreises aus der Kovarianzmatrix neu berechnet. Aber das gibt andere Werte als der ursprüngliche Radius. (es gibt 141.140) !! – Omar14

+1

... sollte die Zahl geteilt durch 100 eine Glocke läuten :) –

Antwort

2

nur durch reinen Reverse Engineering (Ich bin nicht vertraut mehr mit diesem Material), kann ich dies tun:

%// Generate circle 
R = 189; 
t = linspace(0, 2*pi, 1000).'; 
x = R*cos(t); 
y = R*sin(t); 

%// Find the radius? 
[~,L] = eig(cov([x,y])); 

%// ...hmm, seems off by a factor of sqrt(2) 
2*sqrt(diag(L))   

%// so it would come out right when I'd include a factor of 1/2 in the sqrt(): 
2*sqrt(diag(L)/2)   

So lassen Sie sich diese Theorie für die allgemeinen Ellipsen testen:

%// Random radii 
a1 = 1000*rand; 
a2 = 1000*rand; 

%// Random rotation matrix 
R = @(a)[ 
    +cos(a) +sin(a); 
    -sin(a) +cos(a)]; 

%// Generate pionts on the ellipse 
t = linspace(0, 2*pi, 1000).'; 
xy = [a1*cos(t) a2*sin(t);] * R(rand); 

%// Find the deviation from the known radii 
%// (taking care of that the ordering may be different) 
[~,L] = eig(cov(xy)); 
min(abs(1-bsxfun(@rdivide, 2*sqrt(diag(L)/2), [a1 a2; a2 a1])),[],2) 

was immer etwas Akzeptables zurückgibt.

Also, scheint zu arbeiten :) Kann jemand verifizieren, dass das in der Tat richtig ist?

+0

Der Faktor 'sqrt (2)' ist, weil die Kovarianzmatrix aus Punkten entlang des Umfangs der Ellipse berechnet wird, keine feste Ellipse. Die OP-Gleichung gilt für die Kovarianzmatrix einer durchgezogenen Ellipse. –

0

Um auf Rodys Antwort zu erweitern, hat die Kovarianzmatrix für eine durchgezogene Ellipse Eigenwerte, die durch lambda_i = r_i^2/4 gegeben sind. Dies führt zur OP-Gleichung r_i = 2*sqrt(lambda_i).

Bei einer (nicht fest) einer Ellipse, wie in dem Fall des OP, sind die Eigenwerte denen des festen Fall double: lambda_i = r_i^2/2, was zu r_i = sqrt(2*lambda_i) (die Rody der gleich 2*sqrt(lambda_i/2)).

Ich konnte nicht direkt eine Referenz dafür finden, aber die Mathematik der Kovarianzmatrix ist identisch mit der der Trägheitsmomente. On Wikipedia können Sie den Fall des "kreisförmigen Reifens" und "festen Scheibe" sehen, die sich durch den gleichen Faktor von 2 unterscheiden.

Hier ist eine Anpassung von Rodys Test, sowohl die festen als auch die nicht-festen Fälle:

% Radius to test with 
r = rand(1,2); 

% Random rotation matrix 
R = @(a)[+cos(a) +sin(a); 
     -sin(a) +cos(a)]; 

% Generate pionts on the ellipse 
N = 1000; 
t = linspace(0, 2*pi, N).'; 
xy = r.*[cos(t),sin(t)] * R(rand); 
% Compute radii, compare to known radii 
L = eig(cov(xy)); 
r_ = sqrt(2*L)'; 
err = max(abs(1 - sort(r_) ./ sort(r))) 

% Generate points in the ellipse (SOLID CASE) 
N = 10000; 
t = 2*pi * rand(N,1); 
xy = r .* sqrt(rand(N,1)) .* [cos(t),sin(t)] * R(rand); 
% Compute radii, compare to known radii 
L = eig(cov(xy)); 
r_ = 2*sqrt(L)'; 
err_solid = max(abs(1 - sort(r_) ./ sort(r))) 

Wenn Sie diesen Code ausführen, werden Sie Fehler von 1e-3 sehen und ~ 6e-3 (für den festen Fall, dass ich viel mehr Punkte zu erzeugen, da das Gebiet mehr Punkte benötigt dicht genug abgetastet werden; je mehr Punkte, desto kleiner der Fehler.