Determine positions of projected points onto an ellipse

11 ビュー (過去 30 日間)
Salad Box
Salad Box 2022 年 12 月 15 日
コメント済み: Salad Box 2022 年 12 月 20 日
Hi,
I have a set of data points (blue dots below) and I woule like to project them onto an ellipse (with known a, b, x0, y0, and radian or tilt).
I wonder how I can determine the position of each projected points on the ellipse. Any advices/code/peudo code would be lovely. Thank you.

採用された回答

Matt J
Matt J 2022 年 12 月 16 日
編集済み: Matt J 2022 年 12 月 16 日
Use the ellipseprj function below. It requires the download of trustregprob from the file exchange
ab = [2 1]; %ellipse radii
xy0 = [1 1]; %ellipse center
rad=5; R = [cos(rad), -sin(rad); sin(rad),cos(rad)]; %ellipse tilt
xysamp = randn(2,5); %random points to be projected onto the ellipse
xyproj = ellipseprj( R'*diag(1./ab.^2)*R , xy0, xysamp); %the projected points
%%%Visualize
t = linspace(0,2*pi,50)';
xyc = xy0 + ab.*[cos(t), sin(t)]*R; %points for plotting the ellipse
plot(xyc(:,1),xyc(:,2),'b-')
axis equal
grid on
hold on
plot(xysamp(1,:),xysamp(2,:),'ro') %plot random points
plot(xyproj(1,:),xyproj(2,:),'ks') %plot their projections
line([xysamp(1,:);xyproj(1,:)],[xysamp(2,:);xyproj(2,:)],'color','g')
hold off
function Xp=ellipseprj(Q,xc,X0)
%Projects given 2D points onto 2D ellipse
%
% Xp=ellipseprj(Q,xc,X0)
%
%IN:
%
% Q,xc: Ellipse equation matrices. Q is 2x2 and xc is 2x1 such that
% ellipse equation is (y-xc).'*Q*(y-xc)=1
% X0: 2xN matrix of points to be projected
%
%OUT:
%
% Xp: 2xN matrix of projected points
N=size(X0,2);
xc=xc(:);
[Rt,D]=eig(Q);
Rt=real(Rt); D=real(D);
iD=diag(1./diag(D));
iDsqrt=sqrt(iD);
b=-iDsqrt*Rt.'*bsxfun(@minus,xc,X0);
Yp=nan(size(X0));
for i=1:N
Yp(:,i)=trustregprob(iD,b(:,i),1);
end
Xp=bsxfun(@plus, Rt*iDsqrt*Yp,xc);
end
  5 件のコメント
Matt J
Matt J 2022 年 12 月 20 日
@Salad Box, In addition to what John said, it is a pretty basic calculus exercise to compute the tangent to the ellipse at B, and to verify numerically that the tangent is orthogonal to (A-B). I encourage doing that in addition to any visual tests.
Salad Box
Salad Box 2022 年 12 月 20 日
Thank you John and Matt, especially to John, now that explains everything. My axis was distorted I must admit. Well spotted!
Thank you Matt. Your adivce on checking the tangent was very valuable as well providing another way of checking the perpendicularity.
Great answers. Thank you.:)

サインインしてコメントする。

その他の回答 (2 件)

John D'Errico
John D'Errico 2022 年 12 月 16 日
編集済み: John D'Errico 2022 年 12 月 16 日
Simple. Just download my distance2curve utility, found on the file exchange.
For some example data, I'm not sure what you meant by radian as a variable.
t = linspace(0,2*pi,50)';
ab = [2 1];
xy0 = [1 1];
rotmat = @(R) [cos(R), -sin(R);sin(R),cos(R)];
xyc = xy0 + ab.*[cos(t), sin(t)]*rotmat(5);
plot(xyc(:,1),xyc(:,2),'b-')
axis equal
grid on
Now for a given set of points,
xy0 = randn(10,2);
xyproj = distance2curve(xyc,xy0);
hold on
plot(xy0(:,1),xy0(:,2),'ro')
plot(xyproj(:,1),xyproj(:,2),'rs')
line([xy0(:,1),xyproj(:,1)]',[xy0(:,2),xyproj(:,2)]','color','g')
Each point is projected to the nearest point on the ellipse.
distance2curve is found for free download from the file exchange, here:
  4 件のコメント
Torsten
Torsten 2022 年 12 月 19 日
編集済み: Torsten 2022 年 12 月 19 日
Among the t-values you supply, the code chooses the point of xyc that gives the shortest distance to your point.
So choosing a finer t-grid gives a better result.
Change
t = linspace(0,2*pi,50)';
to
t = linspace(0,2*pi,5000)';
e.g.
Salad Box
Salad Box 2022 年 12 月 20 日
Thank you for your answer Torsten. I did what you have suggested to modify t to have 5000 intervals. It hasn't changed results. The projected point B still has the same x and y values as when t has 50 intervals. I am not sure whether t is the key to this.

サインインしてコメントする。


Torsten
Torsten 2022 年 12 月 16 日
編集済み: Torsten 2022 年 12 月 16 日
a = 2;
b = 1; % main axes
x0 = 1;
y0 = 1; % translation
alpha = 5; % rotation angle of ellipse in degrees
px = randn(1,10);
py = randn(1,10); % Sample points to be projected
syms theta xp yp xp_given yp_given
% Generate ellipse equation dependent on theta
x = a*cos(theta);
y = b*sin(theta);
% Solve for theta that minimizes the distance
f = (x-xp)^2 + (y-yp)^2;
df = diff(f,theta);
sol_theta = solve(df==0,theta,'MaxDegree',4)
sol_theta = 
% Transform given points (xp_given,yp_given) to coordinate system of the centered ellipse
% and project on it
P_slash = [cosd(-alpha) -sind(-alpha);sind(-alpha) cosd(-alpha)]*[xp_given - x0 ; yp_given - y0];
sol_theta = subs(sol_theta,[xp yp],[P_slash(1) P_slash(2)])
sol_theta = 
% Transform projected points back to given ellipse
Px_proj = cosd(alpha)*subs(x,theta,sol_theta) - sind(alpha)*subs(y,theta,sol_theta) + x0;
Py_proj = sind(alpha)*subs(x,theta,sol_theta) + cosd(alpha)*subs(y,theta,sol_theta) + y0;
% Numerical example
for i = 1:numel(px)
px_proj = double(subs(Px_proj,[xp_given yp_given],[px(i) py(i)]));
py_proj = double(subs(Py_proj,[xp_given yp_given],[px(i) py(i)]));
idx = abs(imag(px_proj))<1e-6 & abs(imag(py_proj))<1e-6;
px_proj = px_proj(idx);
py_proj = py_proj(idx);
d = sqrt((px(i)-px_proj).^2 + (py(i)-py_proj).^2);
[dist(i),index] = min(d);
pxp(i) = px_proj(index);
pyp(i) = py_proj(index);
end
%Graphical presentation
thetanum = 0:2*pi/100:2*pi;
xx = a*cos(thetanum);
yy = b*sin(thetanum);
xxx = xx*cosd(alpha) - yy*sind(alpha);
yyy = xx*sind(alpha) + yy*cosd(alpha);
xxx = xxx + x0;
yyy = yyy + y0;
plot(xxx,yyy)
hold on
arrayfun(@(i)plot([px(i), real(pxp(i))],[py(i), real(pyp(i))]),1:numel(px))
axis equal
grid on
  2 件のコメント
Salad Box
Salad Box 2022 年 12 月 16 日
編集済み: Salad Box 2022 年 12 月 16 日
Thank you for your answer. I am just wondering what theta, x_given and y_given are.
should I replace x_give and y_given with my own data points?
I assume x and y are real numbers, the thing is that my data set contains 85 data points, and in the case of the ellipse because I used theta = 0:0.01:2*pi, it actually generated 629 points to draw the ellipse.
So I am wondering how I can use
f = (x-x_given)^2 + (y-y_given)^2;
if the number of data points are not even matching (x_given and y_given are 85; x and y are 629)?
Torsten
Torsten 2022 年 12 月 16 日
編集済み: Torsten 2022 年 12 月 16 日
px and py should be your data points.
You must admit 0 <= theta < 2*pi to determine the optimal projection point on the ellipse, but the graphical part is only for illustration - it's not a necessary part of the code.
Inputs are - as said -
a = 2;
b = 1; % main axes
x0 = 1;
y0 = 1; % translation
alpha = 5; % rotation angle of ellipse in degrees
px = randn(1,10);
py = randn(1,10); % Sample points to be projected
and outputs are pxp and pyp - the x-and y-coordinates of the points on the ellipse nearest to the (px /py) - and maybe (if it's of interest) the array "dist" which contains the Euclidean distance between the given data points and their counterparts on the ellipse.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeFunctions についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by