From: Els on
You are correct. The one with the shortest distance to the surface is the one I am interested in. But if that is too much trouble, I can calculate that myself. If I just had the two points, that would make my day.

Thanks again.
From: Torsten Hennig on
> Where the normal vector is perpendicular to the Line
> from P1 to P2 and runs through C and C'.

Ok, so you try to find the point C' on the surface
of the ellipsoid nearest to C.
If c = (cx,cy,cz) and C'=(csx,csy,csz), the problem
is

min: (cx-csx)^2 + (cy-csy)^2 + (cz-csz)^2
s.t.
(csx-xc)^2/rx^2 + (csy-yc)^2/ry^2 + (csz-zc)^2/rz^2 = 1.

Introducing a Lagrange parameter lambda, the solution
is given by a stationary point of the Lagrange function

L(lambda,csx,csy,csz) =
(cx-csx)^2 + (cy-csy)^2 + (cz-csz)^2 - lambda*
( (csx-xc)^2/rx^2 + (csy-yc)^2/ry^2 + (csz-zc)^2/rz^2 - 1)

So calculate first-order derivatives of L with
respect to csx,csy,csz and lambda, set them to zero
and solve the resulting system of 4 equations with
fsolve, e.g.

Best wishes
Torsten.
From: Els on
Dear Torsten,

This sounds great, the solution is near. But I have never calculated a first-order derivative in Matlab before. And the xc,yc and zc, are they the same as cx,cy and cz?

Thanks again.
From: Torsten Hennig on
> Dear Torsten,
>
> This sounds great, the solution is near. But I have
> never calculated a first-order derivative in Matlab
> before. And the xc,yc and zc, are they the same as
> cx,cy and cz?
>

No, (xc,yc,zc) is the centre of the ellipsoid.

> Thanks again.

L(lambda,csx,csy,csz) =
(cx-csx)^2 + (cy-csy)^2 + (cz-csz)^2 - lambda*
( (csx-xc)^2/rx^2 + (csy-yc)^2/ry^2 + (csz-zc)^2/rz^2 - 1)

So calculate first-order derivatives of L with
respect to csx,csy,csz and lambda, set them to zero
and solve the resulting system of 4 equations with
fsolve, e.g.

The first-order derivatives are given by

dL/d(lambda) =
-((csx-xc)^2/rx^2 + (csy-yc)^2/ry^2 + (csz-zc)^2/rz^2 - 1)

dL/d(csx) =
-2*(cx-csx)-2*lambda*(csx-xc)/rx^2

dL/d(csy) =
-2*(cy-csy)-2*lambda*(csy-yc)/ry^2

dL/d(csz) =
-2*(cz-csz)-2*lambda*(csz-zc)/rz^2

Now solve
dL/d(lambda) = 0
dL/d(csx) = 0
dL/d(csy) = 0
dL/d(csz) = 0
for lambda, csx,csy and csz using fsolve.

Best wishes
Torsten.
From: Bruno Luong on
This problem of projection on ellipse in 2D can be solved without optimization toolbox, here is the demo:

% Data, generate ellipse centered at origin, x'*Q*x = c
L=randn(2);
Q = L'*L; % (2 x 2) sdp matrix
c = rand;

% Given point (to be projected on ellipse), it can be inside or outside
P = randn(2,1);

% Engine
[U S V] = svd(Q);
s = diag(S);
A = U*diag(sqrt(c./s));
R = [0 1;
-1 0];
B = R*Q*A;
alpha = -P.'*R*Q*A;
rho = sqrt(s(1)/s(2));
beta = c*(rho-1/rho);

Poly = beta*[0 -2 0 2 0] + ...
alpha(1)*[-1 0 0 0 1] + ...
alpha(2)*[0 2 0 2 0];
tprj = roots(Poly);
tprj = tprj(imag(tprj)==0); % remove imaginary solution
tprj = reshape(tprj, 1, []);

% Candidate of the projections
xprj = A*[(1-tprj.^2); 2*tprj];
xprj = bsxfun(@rdivide, xprj, (1+tprj.^2));

% Graphical check
clf
% Generate point on ellipse
theta = linspace(0,2*pi);
t = tan(theta/2);
x = A*[(1-t.^2); 2*t];
x = bsxfun(@rdivide, x, (1+t.^2));
plot(x(1,:),x(2,:));
axis equal;
hold on;
% Plot P projected to the ellipse
plot(P(1),P(2),'or');
for k=1:size(xprj,2)
plot([xprj(1,k) P(1)],[xprj(2,k) P(2)],'-r');
end

% Bruno