Check if point lies inside error ellipse

12 ビュー (過去 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 件)

カテゴリ

Help Center および File Exchange2-D and 3-D Plots についてさらに検索

製品


リリース

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by