Determine positions of projected points onto an ellipse

Question

Salad Box 2022-12-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1879377-determine-positions-of-projected-points-onto-an-ellipse

评论： 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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Matt J 2022-12-16

Duplicate of,

https://www.mathworks.com/matlabcentral/answers/1873937-projecting-a-dataset-on-an-ellipse?s_tid=srchtitle

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2022-12-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1879377-determine-positions-of-projected-points-onto-an-ellipse#answer_1129562

编辑：Matt J 2022-12-16

在 MATLAB Online 中打开

Use the ellipseprj function below. It requires the download of trustregprob from the file exchange

https://www.mathworks.com/matlabcentral/fileexchange/53191-quadratic-minimization-with-norm-constraint

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 个评论
显示 3更早的评论隐藏 3更早的评论

John D'Errico 2022-12-20

在 MATLAB Online 中打开

@Salad Box very often it can LOOK like a line is not the line of closest approach, if you have not got the axes correct. You NEED to have the two axes with the same units. Otherwise, things will look distorted. For example

t = linspace(0,2*pi);
x = cos(t);
y = sin(t);
plot(x,y,'b-')

Now you should know absolutely that the curve I drew was a perfect circle. But it doe not look like a circle! It appears to be an ellipse. Next, draw a ray, protruding from the center of the circle at 45 degrees.

hold on

plot([0 1],[0,1],'r-')

hold off

Now look carefully at the plot. If you dod, you will come to the conclusion that the red line is NOT perpendicular to the curve. But since the red line is a radial line emanating fro mthe center of the circle, do you agree the line MUST have been perpendicular?

Now, I'll replot the circle and the radial line. But I'll add an extra command.

plot(x,y,'b-',[0 1],[0,1],'r-')

axis equal

Do you see that now the line is seen to be clearly perpendicular? As well, the circle actually LOOKS like a circle!

You need to be careful when you look at a plot like that, and come to conclusions based on the wrong axis scaling.

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 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.:)

请先登录，再进行评论。

Answer 2

John D'Errico 2022-12-16

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1879377-determine-positions-of-projected-points-onto-an-ellipse#answer_1129222

编辑：John D'Errico 2022-12-16

在 MATLAB Online 中打开

distance2curve.m

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:

https://www.mathworks.com/matlabcentral/fileexchange/34869-distance2curve?s_tid=ta_fx_results

4 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2022-12-19

编辑：Torsten 2022-12-19

在 MATLAB Online 中打开

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 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.

请先登录，再进行评论。

Answer 3

Torsten 2022-12-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1879377-determine-positions-of-projected-points-onto-an-ellipse#answer_1129182

编辑：Torsten 2022-12-16

在 MATLAB Online 中打开

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 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 2022-12-16

编辑：Torsten 2022-12-16

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Determine positions of projected points onto an ellipse

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（2 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无

另请参阅

类别

标签

Community Treasure Hunt

Determine positions of projected points onto an ellipse

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（2 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无