Check if point lies inside error ellipse

Question

Oliver Hancock 2021-5-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/823870-check-if-point-lies-inside-error-ellipse

回答： KSSV 2021-5-7

I need to display a message saying "geometric centre (gc) is in 99%/95%/90% confidence ellipse" i.e. use if and elseif to determine which ellipse the point gc lies in (if any)

clear all
I = imread('AM\cropped\09_58_00.jpg');
figure
image(I)
AMfigure = xlsread("AllDataCollated.xlsx","Sheet2");
x = AMfigure(:,1);
yBad = AMfigure(:,2);
y = size(I, 1) - yBad ;
t =  AMfigure(:,5);
hold on
for k=0
    
[rows, columns, numberOfColorChannels] = size(I);
ind = ((k*40)+1):40*(k+1); % indices of data to plot. When k=0 ind equals 1 to 40. When k=1 ind equals 41 to 80  
plot(x(ind),y(ind),'w.','MarkerSize', 7)
xMean = mean(x(ind));
yMean = mean(y(ind));
plot(xMean, yMean, 'g+', 'LineWidth', 1, 'MarkerSize', 6);
disp("mean of data:        x=" + xMean + "  y=" + yMean);
data = [x(ind) y(ind)];
s = std(data);
disp(s)
    
% Calculate the eigenvectors and eigenvalues
covariance = cov(data);
[eigenvec, eigenval ] = eig(covariance);
% Get the index of the largest eigenvector
[largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));
largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);
% Get the largest eigenvalue
largest_eigenval = max(max(eigenval));
% Get the smallest eigenvector and eigenvalue
if(largest_eigenvec_ind_c == 1)
    smallest_eigenval = max(eigenval(:,2));
    smallest_eigenvec = eigenvec(:,2);
else
    smallest_eigenval = max(eigenval(:,1));
    smallest_eigenvec = eigenvec(1,:);
end
% Calculate the angle between the x-axis and the largest eigenvector
angle = atan2(largest_eigenvec(2), largest_eigenvec(1));
% This angle is between -pi and pi.
% Let's shift it such that the angle is between 0 and 2pi
if(angle < 0)
    angle = angle + 2*pi;
end
% Get the coordinates of the data mean
avg = mean(data);
% Get the 90% confidence interval error ellipse
chisquare_val90 = 2.1459;
theta_grid = linspace(0,2*pi);
phi = angle; 
X090=avg(1);
Y090=avg(2);
a=chisquare_val90*sqrt(largest_eigenval);
b=chisquare_val90*sqrt(smallest_eigenval);
% the ellipse in x and y coordinates 
ellipse_x_r  = a*cos( theta_grid );
ellipse_y_r  = b*sin( theta_grid );
%Define a rotation matrix
R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];
%let's rotate the ellipse to some angle phi
r_ellipse90 = [ellipse_x_r;ellipse_y_r]' * R;
% Draw the error ellipse
plot(r_ellipse90(:,1) + X090,r_ellipse90(:,2) + Y090,'g-')
hold on;
% Get the 95% confidence interval error ellipse
chisquare_val95 = 2.4477;
theta_grid = linspace(0,2*pi);
phi = angle; 
X095=avg(1);
Y095=avg(2);
a=chisquare_val95*sqrt(largest_eigenval);
b=chisquare_val95*sqrt(smallest_eigenval);
% the ellipse in x and y coordinates 
ellipse_x_r  = a*cos( theta_grid );
ellipse_y_r  = b*sin( theta_grid );
%Define a rotation matrix
R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];
%let's rotate the ellipse to some angle phi
r_ellipse95 = [ellipse_x_r;ellipse_y_r]' * R;
% Draw the error ellipse
plot(r_ellipse95(:,1) + X095,r_ellipse95(:,2) + Y095,'y-')
hold on;
% Get the 99% confidence interval error ellipse
chisquare_val99 = 3.0348;
theta_grid = linspace(0,2*pi);
phi = angle; 
X099=avg(1);
Y099=avg(2);
a=chisquare_val99*sqrt(largest_eigenval);
b=chisquare_val99*sqrt(smallest_eigenval);
% the ellipse in x and y coordinates 
ellipse_x_r  = a*cos( theta_grid );
ellipse_y_r  = b*sin( theta_grid );
%Define a rotation matrix
R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];
%let's rotate the ellipse to some angle phi
r_ellipse99 = [ellipse_x_r;ellipse_y_r]' * R;
% Draw the error ellipse
plot(r_ellipse99(:,1) + X099,r_ellipse99(:,2) + Y099,'r-')
hold on;
figure2 = xlsread("Ellipse centres","sheet1");
x = figure2(:,6);
y = figure2(:,7);
c = ((k+1):(k+1));
plot(x(c),y(c),'cx','markersize',7)
gcX = x(c);
gcY = y(c);
gc = [x(c), y(c)];
disp("geometric centre at: x=" + x(c) + " y=" + y(c));
time = ((k*40)+1);
title({t(time),"(HHMMSS)"})
filename = "analysis of time " + t(time);
disp(filename)
 
    
end