Effacer les filtres
Effacer les filtres

Check if point lies inside error ellipse

11 vues (au cours des 30 derniers jours)
Oliver Hancock
Oliver Hancock le 7 Mai 2021
Réponse apportée : KSSV le 7 Mai 2021
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

Réponse acceptée

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

Plus de réponses (0)

Catégories

En savoir plus sur 2-D and 3-D Plots dans Help Center et File Exchange

Produits


Version

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by