center, major and minor axis of ellipsoid

30 views (last 30 days)
buthayna alnaalwa
buthayna alnaalwa on 12 Apr 2019
Commented: Walter Roberson on 18 Mar 2021
Here is my code to plot the ellipsoids od 50,80,95,99% confidence intervals from my data (x1_bladder, y1_rectum), and it works, but i don not know how to plot the major and minor axis of this ellipsoid, can anyone help me? Thanks in advance
% my Data
Ith_MinT = [x1_bladder];
Can_MinT = [y1_rectum];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Calculate moments of the data sets %%%
% Observed
Ith_mean = mean(Ith_MinT); % x mean
Can_mean= mean(Can_MinT); % y mean
CV = cov(Ith_MinT,Can_MinT); % the angle of ellipse is determined by the cov of data
[Evec, Eval]=eig(CV); % Eigen values and vectors of covariance matrix
%%% Plot observed multivariate contours %%
% Observed data
xCenter = Ith_mean; % ellipses centered at sample averages, center at x,y data mean
yCenter = Can_mean;
theta = 0 : 0.01 : 2*pi; % angles used for plotting ellipses 0:360 degree
% compute angle for rotation of ellipse
% rotation angle will be angle between x axis and first eigenvector
x_vec= [1 ; 0]; % vector along x-axis
cosrotation =dot(x_vec,Evec(:,1))/(norm(x_vec)*norm(Evec(:,1)));
rotation =pi/2-acos(cosrotation); % rotation angle
R = [sin(rotation) cos(rotation); ...
-cos(rotation) sin(rotation)]; % create a rotation matrix
% create chi squared vector
chisq = [1.368 3.2188 4.605 5.991 9.210]; % percentiles of chi^2 dist df=2% the scale of the ellipsoid depends on confidence interval chosen
% chisq = [1.368 3.2188 4.605 5.991];
% size ellipses for each quantile
for i = 1:length(chisq)
% calculate the radius of the ellipse
xRadius(i)=(chisq(i)*Eval(1,1))^.5; % primary
yRadius(i)=(chisq(i)*Eval(2,2))^.5; % secondary
% lines for plotting ellipse
x{i} = xRadius(i)* cos(theta);
y{i} = yRadius(i) * sin(theta);
% rotate ellipse
rotated_Coords{i} = R*[x{i} ; y{i}];
% center ellipse
x_plot{i}=rotated_Coords{i}(1,:)'+xCenter;
y_plot{i}=rotated_Coords{i}(2,:)'+yCenter;
end
% Set up plot
figure
set(gca,'Title',text('String','Probability Ellipses of M1 full weighted of patient plan',...
'FontAngle', 'italic', 'FontWeight', 'bold'),...
'xlabel',text('String', '$\mathbf{M1 full weighted bladder}$', 'Interpreter', 'latex'),...
'ylabel',text('String', '$\mathbf{M1 full weighted rectum}$', 'Interpreter', 'latex'))
hold on
% Plot data points
plot(Ith_MinT,Can_MinT,'o');% scattered data plot
% Plot contours
for j = 1:length(chisq)
plot(x_plot{j},y_plot{j})
end
legend('Data points','50%', '80%', '90%', '95%', '99%')
hold on

Answers (2)

KSSV
KSSV on 12 Apr 2019
% my Data
Ith_MinT = rand(1,100) ;
Can_MinT = rand(1,100) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Calculate moments of the data sets %%%
% Observed
Ith_mean = mean(Ith_MinT); % x mean
Can_mean= mean(Can_MinT); % y mean
CV = cov(Ith_MinT,Can_MinT); % the angle of ellipse is determined by the cov of data
[Evec, Eval]=eig(CV); % Eigen values and vectors of covariance matrix
%%% Plot observed multivariate contours %%
% Observed data
xCenter = Ith_mean; % ellipses centered at sample averages, center at x,y data mean
yCenter = Can_mean;
theta = 0 : 0.01 : 2*pi; % angles used for plotting ellipses 0:360 degree
% compute angle for rotation of ellipse
% rotation angle will be angle between x axis and first eigenvector
x_vec= [1 ; 0]; % vector along x-axis
cosrotation =dot(x_vec,Evec(:,1))/(norm(x_vec)*norm(Evec(:,1)));
rotation =pi/2-acos(cosrotation); % rotation angle
R = [sin(rotation) cos(rotation); ...
-cos(rotation) sin(rotation)]; % create a rotation matrix
% create chi squared vector
chisq = [1.368 3.2188 4.605 5.991 9.210]; % percentiles of chi^2 dist df=2% the scale of the ellipsoid depends on confidence interval chosen
% chisq = [1.368 3.2188 4.605 5.991];
% size ellipses for each quantile
for i = 1:length(chisq)
% calculate the radius of the ellipse
xRadius(i)=(chisq(i)*Eval(1,1))^.5; % primary
yRadius(i)=(chisq(i)*Eval(2,2))^.5; % secondary
% lines for plotting ellipse
x{i} = xRadius(i)* cos(theta);
y{i} = yRadius(i) * sin(theta);
% rotate ellipse
rotated_Coords{i} = R*[x{i} ; y{i}];
% center ellipse
x_plot{i}=rotated_Coords{i}(1,:)'+xCenter;
y_plot{i}=rotated_Coords{i}(2,:)'+yCenter;
end
% Set up plot
figure
set(gca,'Title',text('String','Probability Ellipses of M1 full weighted of patient plan',...
'FontAngle', 'italic', 'FontWeight', 'bold'),...
'xlabel',text('String', '$\mathbf{M1 full weighted bladder}$', 'Interpreter', 'latex'),...
'ylabel',text('String', '$\mathbf{M1 full weighted rectum}$', 'Interpreter', 'latex'))
hold on
% Plot data points
plot(Ith_MinT,Can_MinT,'o');% scattered data plot
% Plot contours
for j = 1:length(chisq)
P = [x_plot{j} y_plot{j}] ;
plot(x_plot{j},y_plot{j})
% Get center
C = mean(P) ;
plot(C(1),C(2),'*')
% Get distances
d = sqrt(sum((C-[P(:,1) P(:,2)]).^2,2)) ;
M1 = [C(1)-min(d),C(2) ; C(1)+min(d),C(2)] ;
plot(M1(:,1),M1(:,2),'r')
M2 = [C(1),C(2)-max(d) ; C(1) ,C(2)+max(d)] ;
plot(M2(:,1),M2(:,2),'r')
end
legend('Data points','50%', '80%', '90%', '95%', '99%')
hold on
  1 Comment
buthayna alnaalwa
buthayna alnaalwa on 12 Apr 2019
thanks alot for your help!
but this gives me axis without rotation.would you please see the attached grpah.this is the results i gort now !

Sign in to comment.


Waqar Khan
Waqar Khan on 18 Mar 2021
Hello, I need to make ellipse using two points such as X=0.098, and Y=0.765. how can i find center from these X and Y and after that draw an ellipse. Please need help.
  3 Comments

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!