Thread Subject:

Vectors on ellipsoid surface

I plotted an ellipsoid [x,y,z]=ellipsoid (standard in Matlab)

Now I have a point C(x,y,z) which lays within this surface.
I am looking for the projection of C on the surface of the ellipsoid,called C', where the normal vector perpendicular to the surface is in line with the vector between C and C'.

Thanks in advance.

> I plotted an ellipsoid [x,y,z]=ellipsoid (standard
> in Matlab)
>
> Now I have a point C(x,y,z) which lays within this
> surface.

But if C lies within the surface, its projection on
the surface is C itself. Or do you mean the point
opposite to C on the surface of the ellipsoid ?

Or do you mean something else, namely that C lies
within the _ellipsoid_, not within the _surface_ of
the ellipsoid ?


> I am looking for the projection of C on the surface
> of the ellipsoid,called C', where the normal vector
> perpendicular to the surface is in line with the
> vector between C and C'.
>
> Thanks in advance.

Best wishes
Torsten.

I had to be more clear indeed, the point C lies within the ellipsoid, so not on the surface.

Point C is a point on the line between P1 and P2, a line which runs through the ellipsoid.

Point C' is the projection of C on the surface of the ellipsoid, along the normal vector.

Best wishes,
Els

Where the normal vector is perpendicular to the Line from P1 to P2 and runs through C and C'.

> I had to be more clear indeed, the point C lies
> within the ellipsoid, so not on the surface.
>
> Point C is a point on the line between P1 and P2, a
> line which runs through the ellipsoid.
>
> Point C' is the projection of C on the surface of the
> ellipsoid, along the normal vector.
>

I think there will be two projection points.
Which one does not matter ?

> Best wishes,
> Els

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.

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

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.

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

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

Sorry make, it:

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

% Given Point (to be projected on ellipse)
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);
tprj = reshape(tprj, 1, []);

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

% Graphical check
% Generate point on ellipse
theta = linspace(0,2*pi);
x = A*[cos(theta); sin(theta)];
clf;
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

Thanks Bruno for all the effort. Though, my ellipsoid is a 3D figure. How difficult is it to adjust your code for 3D?

"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <ht6ob1$20d$1@fred.mathworks.com>...
> Thanks Bruno for all the effort. Though, my ellipsoid is a 3D figure. How difficult is it to adjust your code for 3D?

Oops, my bad I though it's 2D. It's unclear to to me where as the above method can be extended in 3D.

Follow up the Lagrangian method proposed by Tortsen, allow me to synthesize the problem:

Assuming the ellipsoid is E = { x in R^3 : x'*Q*x = 1 }, Q is (3 x 3) sdp matrix.

c is the point to be projected.

We want to minimize
J(x) = | x - c |^2
such that x'*Q*x = 1 (eqt2)

The Lagrangien is

L(x) = J(x) + lambda * (x'*Q*x - 1)

When you derive wrt lambda, you will find (eqt2).
When you derive wrt x, you'll find

(x - c) + lambda * Q*x = 0

or c = x + lambda*Q*x (eqt3)

Q*x is the normal to the ellipsoid, and the relation (eqt3) is a Lagrange Euler condition. It says nothing more than: c is somewhere on the normal vector placed at x.

Now, if we take the dot product of (eqt3) with x, and because x'*Q*x = 1 we find:

<c,x> =|x|^2 + lambda, or:
lambda = <c,x> - |x|^2.

Replace the value of lambda in (eqt3), we get:

c = x + (<c,x> - |x|^2) * Q*x

This can be solved using FSOLVE. I would probably multiply by inv(Q) to make the problem better conditioned:

(<c,x> - |x|^2) * x + P*x + d = 0

where P = inv(Q), and d := -Q\c.

You could also try to solve the above with Newton, Gauss-Newton, Levenberg-Marquardt, etc...

The function which is needed for FSOLVE is very easy to program (you however should organize the parameters so as to have access to c and Q):

function f = myfun(x)
f = (dot(x,c) - dot(x,x))*x + Q\(x-c)
end

Bruno

Here is the implementation of the idea I put earlier:

% Data, generate ellipse centered at origin, { x: x'*Q*x = 1 }
n = 3; % dimension
L=randn(n);
Q = L'*L;

% Given point in R^n (to be projected on ellipse)
c = randn(n,1);

%%
% Engine
% Solve
% (dot(x,c) - dot(x,x))*x + Q\(x-c) = 0
% x'*Q*x = 1
n = size(Q,1);
x = c;
Qinv = inv(Q);
I = eye(n);
maxiter = 100;
tol = 1e-6;
for k=1:maxiter
    d = c-x;
    g = dot(x,d);
    df = g*x - Q\d;
    M = x*(c-2*x)' + g*I + Qinv;
    P = null(x'*Q);
    y = (M*P)\df;
    xold = x;
    x = x-P*y;
    a = x'*Q*x;
    x = x/sqrt(a);
    dx = x-xold;
    if norm(dx)<tol
        break
    end
end


%%
% Graphical check
% Generate many points on ellipse
[U S V] = svd(Q);
s = diag(S);
radii = sqrt(1./s);
A = U*diag(radii);

clf;
if n==2
    theta = linspace(0,2*pi,181);
    ellipse = A*[cos(theta); sin(theta)];
    
    plot(ellipse(1,:),ellipse(2,:));
    axis equal;
    hold on;
    % Plot c projected to the ellipse
    plot(c(1),c(2),'ok');
    plot([x(1) c(1)],[x(2) c(2)],'-r.');
elseif n==3
    N = 64;
    [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N);
    XYZ = U*[X(:) Y(:) Z(:)]';
    X = reshape(XYZ(1,:),[N N]+1);
    Y = reshape(XYZ(2,:),[N N]+1);
    Z = reshape(XYZ(3,:),[N N]+1);
    surf(X,Y,Z);
    axis equal;
    hold on;
    % Plot c projected to the ellipse
    plot3(c(1),c(2),c(3),'ok');
    plot3([x(1) c(1)],[x(2) c(2)],[x(3) c(3)],'-r.');
end

"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <ht5ldf$7tf$1@fred.mathworks.com>...
> I plotted an ellipsoid [x,y,z]=ellipsoid (standard in Matlab)
>
> Now I have a point C(x,y,z) which lays within this surface.
> I am looking for the projection of C on the surface of the ellipsoid,called C', where the normal vector perpendicular to the surface is in line with the vector between C and C'.
>
> Thanks in advance.

  The following is an analytic solution of your problem. It is not an iterative method but gets its solutions directly using matlab's 'roots' function. It assumes an ellipsoid with center at the origin and axes aligned with coordinates axes, satifying the equation

 x^2/a^2 + y^2/b^2 + z^2/c^2 = 1

The point (p,q,r) can be either inside or outside this ellipsoid. Points (x,y,z) on the ellipsoid are found such that (p,q,r) lies on the line through (x,y,z) orthogonal to the surface of the ellipsoid there. In general there will be six solutions to this problem though frequently some of these will be complex-valued. There are always at least two real solutions.

  It should be noted that if two of the three semi-axes a, b, and c, are equal, giving an oblate or prolate spheroid, then there are in general only five solutions, though again some of these may be complex too. It should be noted that in this case if (p,q,r) lies right on the axis of revolution, infinite rings of solutions about this axis become possible. In the code below in this case the algorithm still produces six solutions, but two of them will be equal - that is, equal to within round off error.

  The algorithm below computes the appropriate coefficients of a sextic polynomial in the variable t

 A*t^6 + B*t^5 + C*t^4 + D*t^3 + E*t^2 + F*t + G = 0

and then computes a set of (x,y,z) points on the surface for each root t according to the equations

 x = a^2*p/(a^2+t)
 y = b^2*q/(b^2+t)
 z = c^2*r/(c^2+t)

The significance of the t variable is that it is the common ratio

 (p-x)/(x/p^2) = (q-y)/(y/q^2) = (r-z)/(z/r^2) = t

between components of the line between (p,q,r) and (x,y,z) and the components of a line through (x,y,z) orthogonal to the surface there.

  The code assumes that a, b, c and p, q, r, have been defined.

- - - - - - - - - - - - - - - - - -
 % The program

 % Compute some convenient intermediate values
 p2 = p^2; q2 = q^2; r2 = r^2;
 a2 = a^2; b2 = b^2; c2 = c^2;
 sbc = b2+c2; sca = c2+a2; sab = a2+b2;
 mbc = b2*c2; mca = c2*a2; mab = a2*b2;
 s1 = a2+b2+c2; s2 = mbc+mca+mab; s3 = a2*b2*c2;

 % Now for the polynomial's actual coefficients
 A = 1;
 B = 2*s1;
 C = s1^2+2*s2-p2*a2-q2*b2-r2*c2;
 D = 2*(s1*s2+s3-p2*a2*sbc-q2*b2*sca-r2*c2*sab);
 E = s2^2+2*s3*s1-p^2*a^2*(sbc^2+2*mbc) ...
     -q2*b2*(sca^2+2*mca)-r2*c2*(sab^2+2*mab);
 F = 2*s3*(s2-p2*sbc-q2*sca-r2*sab);
 G = s3*(s3-p2*mbc-q2*mca-r2*mab);

 % Now find the polynomial's six roots
 t = roots([A,B,C,D,E,F,G]);

 % Derive the corresponding sets of coordinates.
 x = a^2*p./(a^2+t);
 y = b^2*q./(b^2+t);
 z = c^2*r./(c^2+t);
 % The resulting three six-element column vectors give the solutions.
- - - - - - - - - - - - - - - - - -

  It should be pointed out that, depending on the position of the (p,q,r) point relative to the ellipsoid, there can be a comparatively large effect from round off errors. For example, if the ellipsoid is nearly a sphere and the point (p,q,r) is nearly at the origin, all points on the surface come close to being solutions, and this would have the effect of magnifying any round off errors in the solutions actually arrived at. This difficulty is inherent in the nature of the stated problem.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hta8lc$e4f$1@fred.mathworks.com>...
>
>
> The significance of the t variable is that it is the common ratio
>
> (p-x)/(x/p^2) = (q-y)/(y/q^2) = (r-z)/(z/r^2) = t
>

That's a way to go Roger. Very good that you though about this. Bravo.

Bruno

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hta8lc$e4f$1@fred.mathworks.com>...
> ........
> The algorithm below computes the appropriate coefficients of a sextic polynomial in the variable t
>
> A*t^6 + B*t^5 + C*t^4 + D*t^3 + E*t^2 + F*t + G = 0
>
> and then computes a set of (x,y,z) points on the surface for each root t according to the equations
>
> x = a^2*p/(a^2+t)
> y = b^2*q/(b^2+t)
> z = c^2*r/(c^2+t)
>
> The significance of the t variable is that it is the common ratio
>
> (p-x)/(x/p^2) = (q-y)/(y/q^2) = (r-z)/(z/r^2) = t
> ........

  I erred in writing the equation explaining the ratio, t. It should read:

"The significance of the t variable is that it is the common ratio

 (p-x)/(x/a^2) = (q-y)/(y/b^2) = (r-z)/(z/c^2) = t "

This error does not affect the code itself, just the explanation.

Roger Stafford

I have extended Roger's idea of parametrization with respect to Lagrange parameter (he call it "t", I call it "lambda") to n-dimensional ellipsoid. The polynomial to be solved is 2*n in order.

Here is the code:

%%%

% Data, generate ellipse centered at origin, { x: x'*Q*x = 1 }
n = 3; % dimension
L = randn(n);
Q = L'*L;

% Given point in R^n (to be projected on ellipse)
c = 1*randn(n,1);

%%
[U S V] = svd(Q);
s = diag(S);
radii = sqrt(1./s);
A = U*diag(radii);

%%
% Engine

% Change variable
% c = U*d
% x = U*y
d = U'*c;
% solve for lambda
% d-y = lambda*S*y
% y'*S*y = 1

% This are polynomial for (s*lambda+1)^2
pd = [s.^2 2*s ones(size(s))];
a = s.*d.^2;
P = 0;
for i=1:n
    % Pi := Prod pd(j) for j#i
    Pi = 1;
    for j=1:n
        if j~=i
            Pi = conv(Pi, pd(j,:));
        end
    end
    % P = sum(Pi)
    P = P + a(i)*Pi;
end
% PC = Prod pd(j) for j#i
Pc = 1;
for j=1:n
    Pc = conv(Pc, pd(j,:));
end
P = [0 0 P]-Pc;
lambda = roots(P);
lambda = reshape(lambda, 1, []);
lambda = lambda(imag(lambda)==0);

% Compute: y = d / (lambda*s + 1)
ls = bsxfun(@times, lambda, s);
y = bsxfun(@rdivide, d, ls+1);
x = U*y;

% Compute l = |x-c|
l = sqrt(sum(bsxfun(@minus, x, c).^2, 1));
% Where the distance is minimal
[lmin imin] = min(l);
lmin

%%
% Graphical check
% Generate many points on ellipse
clf;
if n==2
    theta = linspace(0,2*pi,181);
    ellipse = A*[cos(theta); sin(theta)];
    
    plot(ellipse(1,:),ellipse(2,:));
    axis equal;
    hold on;
    % Plot c projected to the ellipse
    plot(c(1),c(2),'ok');
    for k=1:size(x,2)
        if k==imin
            linespec = '-r.';
        else
            linespec = '-c.';
        end
        plot([x(1,k) c(1)],[x(2,k) c(2)],linespec);
    end
elseif n==3
    N = 64;
    [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N);
    XYZ = U*[X(:) Y(:) Z(:)]';
    min(sqrt(sum(bsxfun(@minus,XYZ,c).^2,1)))
    X = reshape(XYZ(1,:),[N N]+1);
    Y = reshape(XYZ(2,:),[N N]+1);
    Z = reshape(XYZ(3,:),[N N]+1);
    surf(X,Y,Z,'Edgecolor','none','FaceColor','flat');
    colormap gray;
    axis equal;
    hold on;
    plot3(c(1),c(2),c(3),'ok');
    % Plot c projected to the ellipse
    for k=1:size(x,2)
        if k==imin
            linespec = '-r.';
        else
            linespec = '-c.';
        end
        plot3([x(1,k) c(1)],[x(2,k) c(2)],[x(3,k) c(3)],linespec);
    end
else
    fprintf('Sorry I can''t plot\n');
end

% Bruno

There is a detail about the convolution part is somewhat redundant coded, since there is a lot of overlapping of polynomial product. I have not though how to do it in an optimal manner.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <htaqsk$5h1$1@fred.mathworks.com>...
> I have extended Roger's idea of parametrization with respect to Lagrange parameter (he call it "t", I call it "lambda") to n-dimensional ellipsoid. The polynomial to be solved is 2*n in order.
>
> Here is the code:
>
> %%%
>
> % Data, generate ellipse centered at origin, { x: x'*Q*x = 1 }
> n = 3; % dimension
> L = randn(n);
> Q = L'*L;
>
> % Given point in R^n (to be projected on ellipse)
> c = 1*randn(n,1);
>
> %%
> [U S V] = svd(Q);
> s = diag(S);
> radii = sqrt(1./s);
> A = U*diag(radii);
>
> %%
> % Engine
>
> % Change variable
> % c = U*d
> % x = U*y
> d = U'*c;
> % solve for lambda
> % d-y = lambda*S*y
> % y'*S*y = 1
>
> % This are polynomial for (s*lambda+1)^2
> pd = [s.^2 2*s ones(size(s))];
> a = s.*d.^2;
> P = 0;
> for i=1:n
> % Pi := Prod pd(j) for j#i
> Pi = 1;
> for j=1:n
> if j~=i
> Pi = conv(Pi, pd(j,:));
> end
> end
> % P = sum(Pi)
> P = P + a(i)*Pi;
> end
> % PC = Prod pd(j) for j#i
> Pc = 1;
> for j=1:n
> Pc = conv(Pc, pd(j,:));
> end
> P = [0 0 P]-Pc;
> lambda = roots(P);
> lambda = reshape(lambda, 1, []);
> lambda = lambda(imag(lambda)==0);
>
> % Compute: y = d / (lambda*s + 1)
> ls = bsxfun(@times, lambda, s);
> y = bsxfun(@rdivide, d, ls+1);
> x = U*y;
>
> % Compute l = |x-c|
> l = sqrt(sum(bsxfun(@minus, x, c).^2, 1));
> % Where the distance is minimal
> [lmin imin] = min(l);
> lmin
>
> %%
> % Graphical check
> % Generate many points on ellipse
> clf;
> if n==2
> theta = linspace(0,2*pi,181);
> ellipse = A*[cos(theta); sin(theta)];
>
> plot(ellipse(1,:),ellipse(2,:));
> axis equal;
> hold on;
> % Plot c projected to the ellipse
> plot(c(1),c(2),'ok');
> for k=1:size(x,2)
> if k==imin
> linespec = '-r.';
> else
> linespec = '-c.';
> end
> plot([x(1,k) c(1)],[x(2,k) c(2)],linespec);
> end
> elseif n==3
> N = 64;
> [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N);
> XYZ = U*[X(:) Y(:) Z(:)]';
> min(sqrt(sum(bsxfun(@minus,XYZ,c).^2,1)))
> X = reshape(XYZ(1,:),[N N]+1);
> Y = reshape(XYZ(2,:),[N N]+1);
> Z = reshape(XYZ(3,:),[N N]+1);
> surf(X,Y,Z,'Edgecolor','none','FaceColor','flat');
> colormap gray;
> axis equal;
> hold on;
> plot3(c(1),c(2),c(3),'ok');
> % Plot c projected to the ellipse
> for k=1:size(x,2)
> if k==imin
> linespec = '-r.';
> else
> linespec = '-c.';
> end
> plot3([x(1,k) c(1)],[x(2,k) c(2)],[x(3,k) c(3)],linespec);
> end
> else
> fprintf('Sorry I can''t plot\n');
> end
>
> % Bruno
>
> There is a detail about the convolution part is somewhat redundant coded, since there is a lot of overlapping of polynomial product. I have not though how to do it in an optimal manner.
>
> Bruno

Dear Bruno,

Again, thanks for all the effort, this is great. Now I am trying to incorporate it into my own program. Therefore I am trying to understand your code, which is really advanced for me. But what I am not getting is how to give in P1 and P2 (the two points between which C lies), and is the projected C in 90 degrees of this line?

What should I do if my origin is not on [0,0,0], just change [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N); ?

And what should I do if my axis are not straigt but rotated?

Really hope you can help me with this final step.

Best wishes,

Els

I have to say results can be very impressive when a challenging problem catches the interest of CSSM contributors.

"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <htba3v$iq3$1@fred.mathworks.com>...

> But what I am not getting is how to give in P1 and P2 (the two points between which C lies), and is the projected C in 90 degrees of this line?

Sorry? What are P1 and P2? I quote your post #1

[ ...I plotted an ellipsoid [x,y,z]=ellipsoid (standard in Matlab)

Now I have a point C(x,y,z) which lays within this surface.
I am looking for the projection of C on the surface of the ellipsoid,called C', where the normal vector perpendicular to the surface is in line with the vector between C and C'. ]

My understanding is C is given, and C' is the projection on C on the surface of the ellipsoid to be found. In my code, C is called "c" (lower case); and C' is "x(:,imin)".

>
> What should I do if my origin is not on [0,0,0], just change [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N); ?

No, translate it. Simply subtract all your coordinates by the center of the ellipsoid, then the new coordinates are centered about the origin.

>
> And what should I do if my axis are not straigt but rotated?

In my code the axis ellipsoid is not straight, but oriented in the column vectors of matrix U.

I wrote in the comments that I assumed the ellipsoid is E = { x in R^n : x'*Q*x = 1 }, where Q is a (n x n) sdp matrix, this is general ellipsoid *centered* about the origin. The most general form of ellipsoid is E = { x0 + x in R^n : x'*Q*x = 1 }, x0 in R^n and Q is a (n x n) sdp matrix.

If you know
1) the radii: radii(1), radii(2), and radii(3) (three scalars in a vector)
2) the corresponding orientation U1, U2, U3 (three vectors, they must be orthogonal and normed to 1)
3) and the center x0 (a vector) of your ellipsoid

Then define:
U = [U(:,1) U(:,2) U(:3)]
S = diag(1./radii.^2);
Q = U*S*U'

Before calling my code, just subtract x0 to c, then add it back to the result.

c = c-x0
... % Do my code
%
cprojected = x(:,imin)+x0

Bruno

> > But what I am not getting is how to give in P1 and P2 (the two points between which C lies), and is the projected C in 90 degrees of this line?
>
> Sorry? What are P1 and P2? I quote your post #1
>
> My understanding is C is given, and C' is the projection on C on the surface of the ellipsoid to be found. In my code, C is called "c" (lower case); and C' is "x(:,imin)".
>

What you suggest is correct, only c lies on a line which runs from P1 to P2. This is a line through the ellipsoid. C' is the projection of c on the surface of the ellipsoid, where the angle between the line P1(or P2) to c and c to C' is 90 degrees. So I think here has been a misunderstanding.
> >
> > What should I do if my origin is not on [0,0,0], just change [X Y Z] = ellipsoid(0,0,0,radii(1),radii(2),radii(3),N); ?

> If you know
> 1) the radii: radii(1), radii(2), and radii(3) (three scalars in a vector)
> 2) the corresponding orientation U1, U2, U3 (three vectors, they must be orthogonal and normed to 1)
> 3) and the center x0 (a vector) of your ellipsoid
>
> Then define:
> U = [U(:,1) U(:,2) U(:3)]
> S = diag(1./radii.^2);
> Q = U*S*U'

Is Q even necessary? After this code, I let the engine start. But I have the feeling somethings going wrong. What I have now is this:

C=[1.0150 -1.5380 -6.0150]


cx=1; cy=-5;cz=1; % Center of my ellipsoid


ax=12; ay=12; az=20; % radii of my ellipsoid


radii(1)=ax; radii(2)=ay; radii(3)=az;



U(:,1) = [0.1;0.2;1]; U(:,2) = [0.5;0.3;0.7]; U(:,3) = [0.6;0.6;0.4] % example of

possible axis of my ellipsoid


% Data, generate ellipse centered at origin, { x: x'*Q*x = 1 }

c = C'-[cx; cy; cz]

% c = 1*randn(n,1)


U = [U(:,1) U(:,2) U(:,3)]

S = diag(1./radii.^2)

s = diag(S)

% Q = U*s*U'

A = U*diag(radii);

"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <htbpdd$e9i$1@fred.mathworks.com>...

>
> U(:,1) = [0.1;0.2;1]; U(:,2) = [0.5;0.3;0.7]; U(:,3) = [0.6;0.6;0.4] % example of
>
> possible axis of my ellipsoid
>

Axis of an ellipsoid MUST be orthogonal. You can't just put a random orientation like you did.

Bruno

Ok, you are correct. The orthogonal axis problem I solved now.

Still my question about the P1 and P2 remains. Where can I give in these points? So that a line between P1 and P2 comes into existance, where c lies on. And C' is a projection of the c point perpendicular to this line.

Really hope you can help me to solve this final problem.

"Els " <y.e.t.reeuwijk@student.utwente.nl> wrote in message <htdetc$c98$1@fred.mathworks.com>...
> Ok, you are correct. The orthogonal axis problem I solved now.
>
> Still my question about the P1 and P2 remains. Where can I give in these points? So that a line between P1 and P2 comes into existance, where c lies on. And C' is a projection of the c point perpendicular to this line.

Project the ellipsoid boundary on the plane perpendicular to P1-P2, then use 2D version of the algorithm.

Bruno

Dear Bruno,

Really did my best to understand your code and your explenations, but I am really loosing it now. I just don't get what you mean and how I should do this.

My input values are
P1
P2
[xr,yr,zr] % which are the surface points of the ellipsoid, which is rotated along the z axis. The z axis is along vector t.
[cx,cy,cz] % the center of the ellipsoid
[ax,ay,az] % length of the axis in the unrotated x,y,z direction
C % point on line P1-P2

Output should be c' on the surface of my ellipsoid.

I have the feeling that we are so close to a solution, and really hope you could help me with a code. That would be great. Or just a starter, because you allready have done so much.

Best wishes,

Els

You want C' on (E) and CC' orthogonal to P1P2, right? Let us denote by (P) the plane passing through C and perpendicular to vector P1P2. This plane must contain C' (because you want CC' to be orthogonal to P1P2) and C as well.

This plane (P) is 2D affine space. It intersects the 3D ellipsoid (E) in a ellipse contour. Let us call the contour (F).

So let's summarize:

There is a 2D-plane (P) that contains a point C and a ellipse contour (F). Everything happens in this plane. If you take the projection of C on (F) and call it C', it verifies everything you want. No? If yes, the problem is reduces to 2D problem.

How can you project on (P) and reduce mathematically to two coordinates?

By selecting an orthonormal basis of (P), like this:

B = null((P1(:)-P2(:))').

B is 3 x 2, each columns is a vector // to (P). B is orthonormal.

The projection any point in R^3 in (P) can be carried out by the conversion formulas:
  x -> xi = B'*(x - C) (eqt1);
and the reverse R^2 -> R^3
  xi -> x = C + B*xi (eqt2).

Assuming the ellipsoid equation (E) is { x such that x'*A*x + b'*x + c = 0 } (eqt3)

Injecting (eqt2) in (eqt3), we get:

xi'*B'*A*B*xi + 2*C'*A*B*xi + C'*A*C + b'*B*xi + b'*C + c = 0.

This is nothing than the equation of the ellipse (F), the intersection of (E) and (P) in 2D xi-coordinates:

{ xi : xi'*A2*xi + b2'*xi + c2 = 0 } where

A2 = B'*A*B,
b2 = 2*B'*A*C + B'*b,
c2 = C'*A*C + b'*C + c.

The xi-coordinates of C is (0,0)' (from eqt1).

Now we need to project C = (0,0)' on (F) using the method we have showed you.

That gives the xi-coordinates of C'. Convert in 3D coordinates using (eqt2) and you are done.

If something is not clear you should ask a specific question. I can't guess what is clear what is not or where you get stuck.

Bruno

I clean up the code and better structure it. The convolution part is now carried out in an optimal manner (not at all important for 2D or 3D, but it's much nicer that way). The submission on FEX is here:

http://www.mathworks.com/matlabcentral/fileexchange/27711-euclidian-projection-on-ellipsoid

Bruno

Dear Bruno,

This is amazing, i must agree with Mark. Such a 'simple' question evaluates to a great contribution. I do not know how to thank you. My thesis on kinematics of ligaments in the knee will improve really much with this.

I added this code to EllGeo2Alg

radii = reshape(radii, 1, []);
W = bsxfun(@rdivide, U, radii);
A = W*W';
b = -2*A*x0;
c = x0'*A*x0-1;

B=null((P1(:)-P2(:))') ;
C=[0,0]';
A2=B'*A*B;
b2=2*B'*A*C + B'*b;
c2=C'*A*C + b'*C + C;

So that the output of the function is now A2,b2 and C2.

So by doing the next code, I thought I would get C'.
-------------------------------------------------------------------
[A2 b2 c] = EllGeo2Alg(radii, U, x0,P1,P2)
x = EllPrj(c, radii, U);

% Compute l = |x-c|
l = sqrt(sum(bsxfun(@minus, x, c).^2, 1));
% Where the distance is minimal
[lmin imin] = min(l);
----------------------------------------------------------------------

Only it allready goes wrong in b2=2*B'*A*C + B'*b;
Should I transform C(0,0) to 3D before I do this? Or is my whole order of doing things wrong?

By transforming C(0,0) to 3D I should use eq 2 (C+ B*xi), where I do know B and xi, but what parameter should I choose for C? And how can your file in the file exchanger help me with this?

Hope these questions are specific enough.

Again, thanks for all the effort and help. It is getting much clearer now.

Best wishes,

Els

Here are few mistakes I catch:

>
> I added this code to EllGeo2Alg
>
> radii = reshape(radii, 1, []);
> W = bsxfun(@rdivide, U, radii);
> A = W*W';
> b = -2*A*x0;
> c = x0'*A*x0-1;
>
> B=null((P1(:)-P2(:))') ;
> C=[0,0]'; % <- THIS IS WRONG, you should use 3D coordinates
> A2=B'*A*B;
> b2=2*B'*A*C + B'*b;
> c2=C'*A*C + b'*C + C; % <- This is wrong you should use the lower-case c (of A, b c)

It should goes like this (not check):

radii = reshape(radii, 1, []);
W = bsxfun(@rdivide, U, radii);
A = W*W';
b = -2*A*x0;
c = x0'*A*x0-1;

B=null((P1(:)-P2(:))') ;

A2 = B'*A*B;
b2 = 2*B'*A*C+B'*b;
c2 = C'*A*C + b'*C + c;
[radii U x0] = EllAlg2Geo(A2, b2, c2);

C2 = [0; 0]; % xi-coordinates of C
Cp2 = EllPrj(C2, radii, U, x0); % C'
l = sqrt(sum(Cp2.^2, 1));
% Where the distance is minimal
[lmin imin] = min(l);
Cp2 = Cp2(:,imin);

CP = C + B*Cp2; % <- This is the 3D coordinates of the projection C'

% Bruno