Check if point lies inside error ellipse

32 views (last 30 days)
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

Accepted Answer

KSSV
KSSV on 7 May 2021
REad about inpolygon. This will tell you whether given set of points lies inside the given closed polygon.

More Answers (0)

Categories

Find more on Resizing and Reshaping Matrices in Help Center and File Exchange

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!