Covariance Matrix of an Ellipse

7
votes

I've been trying to solve a problem. I'm surprised I haven't been able to find anything really useful on the net.

I know that from the eigenvalues of the covariance matrix of the ellipse, the major and the minor axis of the ellipse can be computed. As the following: 

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

where a1 and a2 are the major and minor axis, e1 and e2 are the eigenvalues of covariance matrix.

My question is: given an edge points (xi,yi) of the image ellipse, how it possible to find the 2×2 covariance matrix of that ellipse?

2
matlabimage-processingmatrixcovarianceellipse
Shouldn't the matrix be simply the covariance of all xi-s yi-s? – Shai
I'm not sure!. I generate an edge points for a circle with radius 100. Then I defined a p = [xi,yi], where P is an edge points matrix n x 2. I used the matlab command cov(P). I recomputed the radius of the circle from the Covariance Matrix. But the gives different values from the original radius. (it gives 141,140)!! – Omar14
...that number divided by 100 should ring a bell :) – Rody Oldenhuis

2 Answers

4
votes

Just by pure reverse engineering (I'm not familiar anymore with this material), I can do this: 

%// 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, let's test that theory for general ellipses: 

%// 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)

which always returns something acceptably small.

So, seems to work :) Can anyone verify that this is indeed correct?

3
votes

To expand on Rody's answer, the covariance matrix for a solid ellipse has eigenvalues given by lambda_i = r_i^2/4. This leads to OP's equation of r_i = 2*sqrt(lambda_i).

For a (non-solid) ellipse, as in the OP's case, the eigenvalues are double those of the solid case: lambda_i = r_i^2/2, leading to r_i = sqrt(2*lambda_i) (which is equal to Rody's 2*sqrt(lambda_i/2)).

I could not directly find a reference for this, but the math of the covariance matrix is identical to that of the moments of inertia. On Wikipedia you can see the case of the "circular hoop" and "solid disk", which differ by the same factor of 2.

Here is an adaptation of Rody's test, doing both the solid and the non-solid cases:

% 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)))

If you run this code, you'll see errors of 1e-3 and ~6e-3 (for the solid case I generate many more points, as the area needs more points to be sampled densely enough; the more points, the smaller the error).