MATLAB Answers

Check if point lies inside error ellipse

11 ビュー (過去 30 日間)
Oliver Hancock
Oliver Hancock 2021 年 5 月 7 日
回答済み: 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

採用された回答

KSSV
KSSV 2021 年 5 月 7 日
REad about inpolygon. This will tell you whether given set of points lies inside the given closed polygon.

その他の回答 (0 件)

カテゴリ

Community Treasure Hunt

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

Start Hunting!

Translated by