Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
How to specify the ellipsoids by equations alone?

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 18 Nov, 2007 05:56:56

Message: 1 of 26

Dears,

When we are talking about circles and spheres we are only
concerened about the radii and centers. And it is easy to
include these parameters in to the equqtion. Coming to
ellipsoids, we need the orientation matrix that will tell
us how the ellipsoid is positioned. Is there a means to
include the orientation information in to the equation of
an ellipsoid?


thank you,

toja

Subject: How to specify the ellipsoids by equations alone?

From: John D'Errico

Date: 18 Nov, 2007 10:03:10

Message: 2 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message <fhok78$obt
$1@fred.mathworks.com>...

> When we are talking about circles and spheres we are only
> concerened about the radii and centers. And it is easy to
> include these parameters in to the equqtion. Coming to
> ellipsoids, we need the orientation matrix that will tell
> us how the ellipsoid is positioned. Is there a means to
> include the orientation information in to the equation of
> an ellipsoid?

The simple definition of an ellipsoid uses
a matrix formulation. Thus for column
vectors x and x0, and 3x3 matrix H, this
equation defines the surface of the ellipse.

  (x - x0)' * H * (x-x0) = 1

Here x0 is a vector that defines the center
of the ellipse.

H is a 3x3 array that defines the entire
shape of the ellipse, both its size and its
orientation in 3-space. (Note that this
formulation also holds in 2-d, where H
will be a 2x2 matrix.) In either event, the
requirement that this equation describes
an ellipse/ellipsoid is for H to be a positive
definite matrix.

I know, your next question is how do you
use this equation to generate points on
the surface of the ellipse. Or perhaps you
will ask how to generate the matrix H.
Since I'm not sure what you will ask, I'll
stop here for now.

HTH,
John

Subject: How to specify the ellipsoids by equations alone?

From: John D'Errico

Date: 18 Nov, 2007 10:03:10

Message: 3 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message <fhok78$obt
$1@fred.mathworks.com>...

> When we are talking about circles and spheres we are only
> concerened about the radii and centers. And it is easy to
> include these parameters in to the equqtion. Coming to
> ellipsoids, we need the orientation matrix that will tell
> us how the ellipsoid is positioned. Is there a means to
> include the orientation information in to the equation of
> an ellipsoid?

The simple definition of an ellipsoid uses
a matrix formulation. Thus for column
vectors x and x0, and 3x3 matrix H, this
equation defines the surface of the ellipse.

  (x - x0)' * H * (x-x0) = 1

Here x0 is a vector that defines the center
of the ellipse.

H is a 3x3 array that defines the entire
shape of the ellipse, both its size and its
orientation in 3-space. (Note that this
formulation also holds in 2-d, where H
will be a 2x2 matrix.) In either event, the
requirement that this equation describes
an ellipse/ellipsoid is for H to be a positive
definite matrix.

I know, your next question is how do you
use this equation to generate points on
the surface of the ellipse. Or perhaps you
will ask how to generate the matrix H.
Since I'm not sure what you will ask, I'll
stop here for now.

HTH,
John

Subject: How to specify the ellipsoids by equations alone?

From: John D'Errico

Date: 18 Nov, 2007 10:03:10

Message: 4 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message <fhok78$obt
$1@fred.mathworks.com>...

> When we are talking about circles and spheres we are only
> concerened about the radii and centers. And it is easy to
> include these parameters in to the equqtion. Coming to
> ellipsoids, we need the orientation matrix that will tell
> us how the ellipsoid is positioned. Is there a means to
> include the orientation information in to the equation of
> an ellipsoid?

The simple definition of an ellipsoid uses
a matrix formulation. Thus for column
vectors x and x0, and 3x3 matrix H, this
equation defines the surface of the ellipse.

  (x - x0)' * H * (x-x0) = 1

Here x0 is a vector that defines the center
of the ellipse.

H is a 3x3 array that defines the entire
shape of the ellipse, both its size and its
orientation in 3-space. (Note that this
formulation also holds in 2-d, where H
will be a 2x2 matrix.) In either event, the
requirement that this equation describes
an ellipse/ellipsoid is for H to be a positive
definite matrix.

I know, your next question is how do you
use this equation to generate points on
the surface of the ellipse. Or perhaps you
will ask how to generate the matrix H.
Since I'm not sure what you will ask, I'll
stop here for now.

HTH,
John

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 18 Nov, 2007 14:06:57

Message: 5 of 26


Hello John,

Thank you for you reply. By the way your comments are
priceless.


The problem is this. I want to determine the intersection
of a line and an ellipsoid. For that I develped the
following code:


% Points P (x,y,z) on a line defined by two points P1
(x1,y1,z1) and P2 (x2,y2,z2) is described by
% P = P1 + u (P2 - P1) or in each coordinate

 x = x1 + u (x2 - x1);
 y = y1 + u (y2 - y1); (1)
z = z1 + u (z2 - z1);

% An ellipsoid centered at P3 (x3,y3,z3) with radii r1,r2
and r3 is described by

(x - x3)2/r1^2 + (y - y3)2/r2^2 + (z - z3)2/r3^2 = 1

% Substituting the equation of the line into the ellipsoid
gives a quadratic equation of the form

% a*u^2 + b*u + c = 0
% where:

% a = r2^2*r3^2*(x2 - x1)^2 + r1^2*r3^2*(y2 - y1)^2 +
r1^2*r2^2*(z2 - z1)^2;
%
% b = 2*( r2^2*r3*(x2 - x1)* (x1 - x3) + r1^2*r3^2*(y2 - y1)
* (y1 - y3) + r1^2*r2^2*(z2 - z1) *(z1 - z3) );
%
% c = r2^2*r3^2*(x1 - x3)^2+ r1^2*r3^2*(y1 - y3)
^2+r1^2*r2^2*(z1 - z3)^2- (r1^2*r2^2*r3^2);

soving for u and substitution in equation 1 gives me the
intersection points. Unfortunately this is only true when
the ellipse has almost no orientation. How could I modify
the equation of the ellipsoid to account for orientation?


that is the challenge?


sincerely,

toja






"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <fhp2ku$4ek$1@fred.mathworks.com>...
> "Toja Debela" <tojadeb@mathworks.com> wrote in message
<fhok78$obt
> $1@fred.mathworks.com>...
>
> > When we are talking about circles and spheres we are
only
> > concerened about the radii and centers. And it is easy
to
> > include these parameters in to the equqtion. Coming to
> > ellipsoids, we need the orientation matrix that will
tell
> > us how the ellipsoid is positioned. Is there a means to
> > include the orientation information in to the equation
of
> > an ellipsoid?
>
> The simple definition of an ellipsoid uses
> a matrix formulation. Thus for column
> vectors x and x0, and 3x3 matrix H, this
> equation defines the surface of the ellipse.
>
> (x - x0)' * H * (x-x0) = 1
>
> Here x0 is a vector that defines the center
> of the ellipse.
>
> H is a 3x3 array that defines the entire
> shape of the ellipse, both its size and its
> orientation in 3-space. (Note that this
> formulation also holds in 2-d, where H
> will be a 2x2 matrix.) In either event, the
> requirement that this equation describes
> an ellipse/ellipsoid is for H to be a positive
> definite matrix.
>
> I know, your next question is how do you
> use this equation to generate points on
> the surface of the ellipse. Or perhaps you
> will ask how to generate the matrix H.
> Since I'm not sure what you will ask, I'll
> stop here for now.
>
> HTH,
> John

Subject: How to specify the ellipsoids by equations alone?

From: Bruno Luong

Date: 18 Nov, 2007 14:20:04

Message: 6 of 26

Parametrize (linearly) your line by a scalar t:

x = M + t.d ( M and d are fixed vectors in R^3)

put it in this equation (H is symetric definite positive
matix).

(x - x0)' * H * (x-x0) = 1

You'll get a second-order polynomial equation with respect to t.

Solve it, and you get your intersection parameters.

Bruno

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 18 Nov, 2007 14:32:38

Message: 7 of 26


Hello, Could you please detail it a bit. I even do not have
clear picture of how you could have come up to x=M+t.d

and how could I get a polynomial out of (x - x0)' * H * (x-
x0) = 1

please clarify it for mr,

Thank you

toja




"Bruno Luong" <brunoluong@yahoo.com> wrote in message
<fhphmk$f27$1@fred.mathworks.com>...
> Parametrize (linearly) your line by a scalar t:
>
> x = M + t.d ( M and d are fixed vectors in R^3)
>
> put it in this equation (H is symetric definite positive
> matix).
>
> (x - x0)' * H * (x-x0) = 1
>
> You'll get a second-order polynomial equation with
respect to t.
>
> Solve it, and you get your intersection parameters.
>
> Bruno

Subject: How to specify the ellipsoids by equations alone?

From: Bruno Luong

Date: 18 Nov, 2007 14:45:28

Message: 8 of 26

> I even do not have
> clear picture of how you could have come up to x=M+t.d

if your line go through 2 points P1 and P2, you can chose M
as P1 and d=P2-P1. t is an unknown scalar.

>
> and how could I get a polynomial out of (x - x0)' * H * (x-
> x0) = 1
>
> please clarify it for mr,
>

Replace x=M+t.d in to the equation and do the algebra ...

(x-x0) = M-x0 + t.d

(x-x0)'*H*(x-x0) = 1

t^2 (d' * H * d) + t * ( 2 . (M-x0)' * H * d) +
(M-x0)'*H*(M-x0) - 1 = 0

This is a polynomial in t to be solved.

Bruno

Subject: How to specify the ellipsoids by equations alone?

From: Bruno Luong

Date: 18 Nov, 2007 14:50:30

Message: 9 of 26

> I even do not have
> clear picture of how you could have come up to x=M+t.d

if your line go through 2 points P1 and P2, you can chose M
as P1 and d=P2-P1. t is an unknown scalar.

>
> and how could I get a polynomial out of (x - x0)' * H * (x-
> x0) = 1
>
> please clarify it for mr,
>

Replace x=M+t.d in to the equation and do the algebra ...

(x-x0) = M-x0 + t.d

(x-x0)'*H*(x-x0) = 1

t^2 (d' * H * d) + t * ( 2 . (M-x0)' * H * d) +
(M-x0)'*H*(M-x0) - 1 = 0

This is a polynomial in t to be solved.

Bruno

Subject: How to specify the ellipsoids by equations alone?

From: John D'Errico

Date: 18 Nov, 2007 15:50:44

Message: 10 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message <fhpgu1$5a4
$1@fred.mathworks.com>...
>
> Hello John,
>
> Thank you for you reply. By the way your comments are
> priceless.
>
>
> The problem is this. I want to determine the intersection
> of a line and an ellipsoid. For that I develped the
> following code:
>
>
> % Points P (x,y,z) on a line defined by two points P1
> (x1,y1,z1) and P2 (x2,y2,z2) is described by
> % P = P1 + u (P2 - P1) or in each coordinate
>
> x = x1 + u (x2 - x1);
> y = y1 + u (y2 - y1); (1)
> z = z1 + u (z2 - z1);


Better is to leave this in vector form.
Parameterized by the scalar value u,
we have

  P(u) = P1 + u*P2

where P is any point on the line, and
P1, P2 are 3x1 vectors in Matlab.

 
> % An ellipsoid centered at P3 (x3,y3,z3) with radii r1,r2
> and r3 is described by
>
> (x - x3)2/r1^2 + (y - y3)2/r2^2 + (z - z3)2/r3^2 = 1
>
> % Substituting the equation of the line into the ellipsoid
> gives a quadratic equation of the form

(snip)

> soving for u and substitution in equation 1 gives me the
> intersection points. Unfortunately this is only true when
> the ellipse has almost no orientation. How could I modify
> the equation of the ellipsoid to account for orientation?
>
>
> that is the challenge?


Return to the formulation I posed before.

  (X - X0)' * H * (X - X0) = 1

If X is a solution to the ellipsoid equation
and is also a point on the line, then we
know that

  (P1 + u*P2 - X0)' * H * (P1 + u*P2 - X0) = 1

Expand this expression, collecting terms.
You will have a quadratic polynomial with
scalar coefficients in u. Solve for the two
solutions. If real solutions do not exist,
then the line does not intersect your
ellipsoid. If you find real solutions, then
substitute back into your equation of the
line.

  u^2*(P2'*H*P2) + ...
  
  2*u*(P2'*(P1-X0)) + ...

  (P1-X0)'*H*(P1-X0) - 1 = 0

Alternatively, transform the line. Given H
as a 3x3 positive definite matrix, compute
the cholesky decomposition:

  H = L'*L

using the chol function. Then the transformed
line is given by

  Z(u) = L * (P1 + u*P2 - X0)

         = Z1 + u*Z2

where

  Z1 = L*(P1 - X0)

and
  
  Z2 = L*P2

Now solve the transformed problem (note
that I have translated the origin to X0 here.)

  Z'*Z = 1

so that

  (Z1 + u*Z2)' * (Z1 + u*Z2) = 1

Expanding this as we did above, we get

  u^2*(Z2'*Z2) + 2*u*Z2'*Z1 + (Z1'*Z1 - 1) = 0

Again, solve for the roots of this quadratic
polynomial in u. u is the same scalar as it
was before, so we can compute the points
P(u) directly in the untransformed domain.

The question that remains is what is H? For
the simple ellipsoid with axes oriented with
the coordinate axes, H is a diagonal matrix
with diagonal elements 1./[r1^2, r2^2, r3^2].
Thus

  H = diag(1./[r1^2, r2^2, r3^2]);

In that event, we can also write L easily
without even recourse to chol.

  L = diag(1./[r1, r2, r3]);
 
If your ellipsoid has orientation other than
this, it is still easy enough to define H.

John

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 20 Nov, 2007 17:18:22

Message: 11 of 26

Hello John

Thank you again.

If you allow me to extend my question.

%%%

X=P1+u*(P2-P1)

 Return to the formulation You posed before.

 (X - X0)' * H * (X - X0) = 1


 If X is a solution to the ellipsoid equation
 and is also a point on the line, then we
know that
 
   (P1 + u*(P2-P1) - X0)' * H * (P1 + u*(P2-P1) - X0) = 1

%%%

Suppose that I have P1 is the centroid of a set of points
P2,say they are 100. I want to find H that minimizes the
the sum of square of distance differences,sum ((P3-P2)^2).

How would I approach that?











"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <fhpn0k$nv4$1@fred.mathworks.com>...
> "Toja Debela" <tojadeb@mathworks.com> wrote in message
<fhpgu1$5a4
> $1@fred.mathworks.com>...
> >
> > Hello John,
> >
> > Thank you for you reply. By the way your comments are
> > priceless.
> >
> >
> > The problem is this. I want to determine the
intersection
> > of a line and an ellipsoid. For that I develped the
> > following code:
> >
> >
> > % Points P (x,y,z) on a line defined by two points P1
> > (x1,y1,z1) and P2 (x2,y2,z2) is described by
> > % P = P1 + u (P2 - P1) or in each coordinate
> >
> > x = x1 + u (x2 - x1);
> > y = y1 + u (y2 - y1); (1)
> > z = z1 + u (z2 - z1);
>
>
> Better is to leave this in vector form.
> Parameterized by the scalar value u,
> we have
>
> P(u) = P1 + u*P2
>
> where P is any point on the line, and
> P1, P2 are 3x1 vectors in Matlab.
>
>
> > % An ellipsoid centered at P3 (x3,y3,z3) with radii
r1,r2
> > and r3 is described by
> >
> > (x - x3)2/r1^2 + (y - y3)2/r2^2 + (z - z3)2/r3^2 = 1
> >
> > % Substituting the equation of the line into the
ellipsoid
> > gives a quadratic equation of the form
>
> (snip)
>
> > soving for u and substitution in equation 1 gives me
the
> > intersection points. Unfortunately this is only true
when
> > the ellipse has almost no orientation. How could I
modify
> > the equation of the ellipsoid to account for
orientation?
> >
> >
> > that is the challenge?
>
>
> Return to the formulation I posed before.
>
> (X - X0)' * H * (X - X0) = 1
>
> If X is a solution to the ellipsoid equation
> and is also a point on the line, then we
> know that
>
> (P1 + u*P2 - X0)' * H * (P1 + u*P2 - X0) = 1
>
> Expand this expression, collecting terms.
> You will have a quadratic polynomial with
> scalar coefficients in u. Solve for the two
> solutions. If real solutions do not exist,
> then the line does not intersect your
> ellipsoid. If you find real solutions, then
> substitute back into your equation of the
> line.
>
> u^2*(P2'*H*P2) + ...
>
> 2*u*(P2'*(P1-X0)) + ...
>
> (P1-X0)'*H*(P1-X0) - 1 = 0
>
> Alternatively, transform the line. Given H
> as a 3x3 positive definite matrix, compute
> the cholesky decomposition:
>
> H = L'*L
>
> using the chol function. Then the transformed
> line is given by
>
> Z(u) = L * (P1 + u*P2 - X0)
>
> = Z1 + u*Z2
>
> where
>
> Z1 = L*(P1 - X0)
>
> and
>
> Z2 = L*P2
>
> Now solve the transformed problem (note
> that I have translated the origin to X0 here.)
>
> Z'*Z = 1
>
> so that
>
> (Z1 + u*Z2)' * (Z1 + u*Z2) = 1
>
> Expanding this as we did above, we get
>
> u^2*(Z2'*Z2) + 2*u*Z2'*Z1 + (Z1'*Z1 - 1) = 0
>
> Again, solve for the roots of this quadratic
> polynomial in u. u is the same scalar as it
> was before, so we can compute the points
> P(u) directly in the untransformed domain.
>
> The question that remains is what is H? For
> the simple ellipsoid with axes oriented with
> the coordinate axes, H is a diagonal matrix
> with diagonal elements 1./[r1^2, r2^2, r3^2].
> Thus
>
> H = diag(1./[r1^2, r2^2, r3^2]);
>
> In that event, we can also write L easily
> without even recourse to chol.
>
> L = diag(1./[r1, r2, r3]);
>
> If your ellipsoid has orientation other than
> this, it is still easy enough to define H.
>
> John

Subject: How to specify the ellipsoids by equations alone?

From: John D'Errico

Date: 20 Nov, 2007 18:27:54

Message: 12 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message
<fhv4su$49d$1@fred.mathworks.com>...
> Hello John
>
> Thank you again.
>
> If you allow me to extend my question.
>
> %%%
>
> X=P1+u*(P2-P1)
>
> Return to the formulation You posed before.
>
> (X - X0)' * H * (X - X0) = 1
>
>
> If X is a solution to the ellipsoid equation
> and is also a point on the line, then we
> know that
>
> (P1 + u*(P2-P1) - X0)' * H * (P1 + u*(P2-P1) - X0) = 1
>
> %%%
>
> Suppose that I have P1 is the centroid of a set of points
> P2,say they are 100. I want to find H that minimizes the
> the sum of square of distance differences,sum ((P3-P2)^2).
>
> How would I approach that?


I'm not sure what you mean here.

The sum of squared differences from a line
is not dependent on H. Said differently, H does
not define a line. It defines an ellipse. (Do you
wish to compute weighted distances?)

If you want to find that line which minimizes
the sum of squared differences, it is easily
computed from the SVD (singular value
decomposition). This problem is usually known
as total least squares, or perhaps as the errors
in variables problem.

For example, given an array X of data in 3-d.
So X is an nx3 array, each row of which is one
data point.

The line in question can be written as

  P(t) = P1 + P2*t

Here P1 will be the global mean of your data,

  P1 = mean(X,1);

P2 comes from the first right singular vector
of X.

  [U,S,V] = svd(X-repmat(P1,n,1),0);

  P2 = V(:,1)';

HTH,
John

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 20 Nov, 2007 19:14:38

Message: 13 of 26


Basically I am looking for the ellipsoid minimizes the sum
of square distances between my point sets and the ellipsoid
point. It all about optimizing an ellipsoid. If we have the
center point of the set of points(P1), the points them
selves(P2), we can callculate P3 both on lines connecting
P1 and P2 and on the ellipsoid. P1 and P2 are fixed. so,P3
is the function of H ( ellipsoid center fixed).

In our previous posts, we calculated P3 given H, my
question then is: starting with an initial H=H0,trying to
find H such that

  Optimize H such that sum(square(P3-P2)) is minimum



Am I clear? Does the question make sense?


Thank you John

Toja











"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <fhv8va$1rn$1@fred.mathworks.com>...
> "Toja Debela" <tojadeb@mathworks.com> wrote in message b
> <fhv4su$49d$1@fred.mathworks.com>...
> > Hello John
> >
> > Thank you again.
> >
> > If you allow me to extend my question.
> >
> > %%%
> >
> > X=P1+u*(P2-P1)
> >
> > Return to the formulation You posed before.
> >
> > (X - X0)' * H * (X - X0) = 1
> >
> >
> > If X is a solution to the ellipsoid equation
> > and is also a point on the line, then we
> > know that
> >
> > (P1 + u*(P2-P1) - X0)' * H * (P1 + u*(P2-P1) - X0) =
1
> >
> > %%%
> >
> > Suppose that I have P1 is the centroid of a set of
points
> > P2,say they are 100. I want to find H that minimizes
the
> > the sum of square of distance differences,sum ((P3-P2)
^2).
> >
> > How would I approach that?
>
>
> I'm not sure what you mean here.
>
> The sum of squared differences from a line
> is not dependent on H. Said differently, H does
> not define a line. It defines an ellipse. (Do you
> wish to compute weighted distances?)
>
> If you want to find that line which minimizes
> the sum of squared differences, it is easily
> computed from the SVD (singular value
> decomposition). This problem is usually known
> as total least squares, or perhaps as the errors
> in variables problem.
>
> For example, given an array X of data in 3-d.
> So X is an nx3 array, each row of which is one
> data point.
>
> The line in question can be written as
>
> P(t) = P1 + P2*t
>
> Here P1 will be the global mean of your data,
>
> P1 = mean(X,1);
>
> P2 comes from the first right singular vector
> of X.
>
> [U,S,V] = svd(X-repmat(P1,n,1),0);
>
> P2 = V(:,1)';
>
> HTH,
> John
>
>

Subject: How to specify the ellipsoids by equations alone?

From: Bruno Luong

Date: 20 Nov, 2007 19:29:10

Message: 14 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message
<fhvbmu$g9b$1@fred.mathworks.com>...
>
> Basically I am looking for the ellipsoid minimizes the sum
> of square distances between my point sets and the ellipsoid
> point. It all about optimizing an ellipsoid. If we have the
> center point of the set of points(P1), the points them
> selves(P2), we can callculate P3 both on lines connecting
> P1 and P2 and on the ellipsoid. P1 and P2 are fixed. so,P3
> is the function of H ( ellipsoid center fixed).
>
> In our previous posts, we calculated P3 given H, my
> question then is: starting with an initial H=H0,trying to
> find H such that
>
> Optimize H such that sum(square(P3-P2)) is minimum
>
>

Two remarks :
 
1. Most people does the ellipsoide fitting like following:

Find H Minimizing J(H)=sum |d(P2,E(H))|^2

where E(H) is the ellipsoide defined by H,
and d(P2,E) := the distance of P2 to E (i.e., d=min{|P2-X|:
X in E} ).

Your formulation of using P3/P1 seems unusual to me, but why
not?

2. Finding an ellipsoid that fits a set of points (even in a
least-square sense) is a non-linear optimization problem.
Implementing a fitting algorithm by a circle is not quite
straight-forward. There are many scientific papers written
about this topic.

Bruno

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 20 Nov, 2007 20:11:05

Message: 15 of 26


Hi, Bruno

Thank you for your remarks. The use of P1 is to find P3. I
needed P3 because I could not find a way that I can find
what you defined as X. X in E. Do you know litratures on
this subject? Most peopele use minimum vomume enclosing
ellipsoids, which I am not satisfied with.

thanks

toja








"Bruno Luong" <brunoluong@yahoo.com> wrote in message
<fhvci6$14c$1@fred.mathworks.com>...
> "Toja Debela" <tojadeb@mathworks.com> wrote in message
> <fhvbmu$g9b$1@fred.mathworks.com>...
> >
> > Basically I am looking for the ellipsoid minimizes the
sum
> > of square distances between my point sets and the
ellipsoid
> > point. It all about optimizing an ellipsoid. If we have
the
> > center point of the set of points(P1), the points them
> > selves(P2), we can callculate P3 both on lines
connecting
> > P1 and P2 and on the ellipsoid. P1 and P2 are fixed.
so,P3
> > is the function of H ( ellipsoid center fixed).
> >
> > In our previous posts, we calculated P3 given H, my
> > question then is: starting with an initial H=H0,trying
to
> > find H such that
> >
> > Optimize H such that sum(square(P3-P2)) is minimum
> >
> >
>
> Two remarks :
>
> 1. Most people does the ellipsoide fitting like following:
>
> Find H Minimizing J(H)=sum |d(P2,E(H))|^2
>
> where E(H) is the ellipsoide defined by H,
> and d(P2,E) := the distance of P2 to E (i.e., d=min{|P2-
X|:
> X in E} ).
>
> Your formulation of using P3/P1 seems unusual to me, but
why
> not?
>
> 2. Finding an ellipsoid that fits a set of points (even
in a
> least-square sense) is a non-linear optimization problem.
> Implementing a fitting algorithm by a circle is not quite
> straight-forward. There are many scientific papers written
> about this topic.
>
> Bruno

Subject: How to specify the ellipsoids by equations alone?

From: John D'Errico

Date: 20 Nov, 2007 20:38:18

Message: 16 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message
<fhvbmu$g9b$1@fred.mathworks.com>...
>
> Basically I am looking for the ellipsoid minimizes the sum
> of square distances between my point sets and the ellipsoid
> point. It all about optimizing an ellipsoid. If we have the
> center point of the set of points(P1), the points them
> selves(P2), we can callculate P3 both on lines connecting
> P1 and P2 and on the ellipsoid. P1 and P2 are fixed. so,P3
> is the function of H ( ellipsoid center fixed).
>
> In our previous posts, we calculated P3 given H, my
> question then is: starting with an initial H=H0,trying to
> find H such that
>
> Optimize H such that sum(square(P3-P2)) is minimum
>
>
>
> Am I clear? Does the question make sense?

Ok. This is not what I originally thought that you
were looking for.

So you have a set of points P2, and a single point,
P1 in R^3. For each point in P2, given an H that
defines an ellipse, you wish to project the points
in P2 onto the ellipsoid along the ray starting at
P1. Call the set of projected points P3. Now you
wish to find the ellipsoid parameters (H) that
minimize the squared distances between P2 and
P3.

First, start by translating the problem to the origin.
That is, subtract off the point P1 from each point
in P2. In effect, we now only need to worry about
finding H, since P1 is given. (Or do you mean to
find the optimal P1 also?)

I'd use the transformation approach myself at this
point. Thus, given H, compute a cholesky factor
then transform each translated point (P2 - P1).
Then P3 is trivially computed for each point in
the transformed space, and transform back to get
the set P3 in the original domain.

Now its just an optimization problem, on the
6 parameters that make up H. (Note that H is
symmetric, positive definite to define an ellipse.)
Set an optimizer to work on it, perhaps fmincon,
since you have a nonlinear inequality constraint.

John

Subject: How to specify the ellipsoids by equations alone?

From: Bruno Luong

Date: 20 Nov, 2007 20:38:51

Message: 17 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message
<fhvf0p$6qa$1@fred.mathworks.com>...
> Do you know litratures on
> this subject?
>

Toja,

Humm, hard copy archive of my papers is not quite handly.
Not sure where I put them in among the mess in the house.
Tomorrow if I find the reference in the office I'll let you
know.

Bruno

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 21 Nov, 2007 14:03:17

Message: 18 of 26

Hello John

You inspired me and finally managed to find the optimized
ellipsoid except the orientation. But you also hinted that
the optimization problem is performed, on the
6 parameters that make up H.

I used

 H=diag(1./[r1^2,r2^2,r3^2]). Unfortunately, this matrix
does not give the real orientation of the points I want to
fit ellipsoids. It is only a matter of orientation!!!!

The rest is done.... Do you have an idea how I could
include the orientation in the symetric matrix?

Thanks to you and Burno









"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <fhvgjq$t4p$1@fred.mathworks.com>...
> "Toja Debela" <tojadeb@mathworks.com> wrote in message
> <fhvbmu$g9b$1@fred.mathworks.com>...
> >
> > Basically I am looking for the ellipsoid minimizes the
sum
> > of square distances between my point sets and the
ellipsoid
> > point. It all about optimizing an ellipsoid. If we have
the
> > center point of the set of points(P1), the points them
> > selves(P2), we can callculate P3 both on lines
connecting
> > P1 and P2 and on the ellipsoid. P1 and P2 are fixed.
so,P3
> > is the function of H ( ellipsoid center fixed).
> >
> > In our previous posts, we calculated P3 given H, my
> > question then is: starting with an initial H=H0,trying
to
> > find H such that
> >
> > Optimize H such that sum(square(P3-P2)) is minimum
> >
> >
> >
> > Am I clear? Does the question make sense?
>
> Ok. This is not what I originally thought that you
> were looking for.
>
> So you have a set of points P2, and a single point,
> P1 in R^3. For each point in P2, given an H that
> defines an ellipse, you wish to project the points
> in P2 onto the ellipsoid along the ray starting at
> P1. Call the set of projected points P3. Now you
> wish to find the ellipsoid parameters (H) that
> minimize the squared distances between P2 and
> P3.
>
> First, start by translating the problem to the origin.
> That is, subtract off the point P1 from each point
> in P2. In effect, we now only need to worry about
> finding H, since P1 is given. (Or do you mean to
> find the optimal P1 also?)
>
> I'd use the transformation approach myself at this
> point. Thus, given H, compute a cholesky factor
> then transform each translated point (P2 - P1).
> Then P3 is trivially computed for each point in
> the transformed space, and transform back to get
> the set P3 in the original domain.
>
> Now its just an optimization problem, on the
> 6 parameters that make up H. (Note that H is
> symmetric, positive definite to define an ellipse.)
> Set an optimizer to work on it, perhaps fmincon,
> since you have a nonlinear inequality constraint.
>
> John

Subject: How to specify the ellipsoids by equations alone?

From: John D'Errico

Date: 21 Nov, 2007 15:04:27

Message: 19 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message
<fi1dr5$c5g$1@fred.mathworks.com>...
> Hello John
>
> You inspired me and finally managed to find the optimized
> ellipsoid except the orientation. But you also hinted that
> the optimization problem is performed, on the
> 6 parameters that make up H.
>
> I used
>
> H=diag(1./[r1^2,r2^2,r3^2]). Unfortunately, this matrix
> does not give the real orientation of the points I want to
> fit ellipsoids. It is only a matter of orientation!!!!
>
> The rest is done.... Do you have an idea how I could
> include the orientation in the symetric matrix?
>
> Thanks to you and Burno

Good! You are getting there. You might try
this:

H = reshape(params([1 2 3 2 4 5 3 5 6]),3,3);

This matrix will be symmetric, but it may have
negative eigenvalues. The problem is that you
need to constrain the eigenvalues to be
non-negative. So better would be to form H
as

L = [params([1]) 0 0;params(2:3);params(4:6)];
H = L*L';

Note that H is ALWAYS positive definite by this
definition, but also this gives you a direct
factorization of H.

John

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 21 Nov, 2007 17:05:17

Message: 20 of 26

Hello John,

ok, but my problem was not making H positive definite
matrix. Let me state it again:

-I have a set of points, say 200
- I am optimizing H in the least square distance
difference between the points and the ellipsoid. That is
searching the optimal radii r1,r2,r3. But the radii is not
enough to define the ellipsoid. I want to determine the
orientation information from my set of points and include
in H,

Do you have an idea?

Thank you


toja




"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <fi1hdr$1f5$1@fred.mathworks.com>...
> "Toja Debela" <tojadeb@mathworks.com> wrote in message
> <fi1dr5$c5g$1@fred.mathworks.com>...
> > Hello John
> >
> > You inspired me and finally managed to find the
optimized
> > ellipsoid except the orientation. But you also hinted
that
> > the optimization problem is performed, on the
> > 6 parameters that make up H.
> >
> > I used
> >
> > H=diag(1./[r1^2,r2^2,r3^2]). Unfortunately, this
matrix
> > does not give the real orientation of the points I
want to
> > fit ellipsoids. It is only a matter of orientation!!!!
> >
> > The rest is done.... Do you have an idea how I could
> > include the orientation in the symetric matrix?
> >
> > Thanks to you and Burno
>
> Good! You are getting there. You might try
> this:
>
> H = reshape(params([1 2 3 2 4 5 3 5 6]),3,3);
>
> This matrix will be symmetric, but it may have
> negative eigenvalues. The problem is that you
> need to constrain the eigenvalues to be
> non-negative. So better would be to form H
> as
>
> L = [params([1]) 0 0;params(2:3);params(4:6)];
> H = L*L';
>
> Note that H is ALWAYS positive definite by this
> definition, but also this gives you a direct
> factorization of H.
>
> John

Subject: How to specify the ellipsoids by equations alone?

From: Justin Abbott

Date: 21 Nov, 2007 17:52:26

Message: 21 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message
<fi1ogd$e04$1@fred.mathworks.com>...
> Hello John,
>
> ok, but my problem was not making H positive definite
> matrix. Let me state it again:
>
> -I have a set of points, say 200
> - I am optimizing H in the least square distance
> difference between the points and the ellipsoid. That is
> searching the optimal radii r1,r2,r3. But the radii is
not
> enough to define the ellipsoid. I want to determine the
> orientation information from my set of points and
include
> in H,
>
> Do you have an idea?
>
> Thank you
>
>
> toja
>
>
>
>
> "John D'Errico" <woodchips@rochester.rr.com> wrote in
> message <fi1hdr$1f5$1@fred.mathworks.com>...
> > "Toja Debela" <tojadeb@mathworks.com> wrote in message
> > <fi1dr5$c5g$1@fred.mathworks.com>...
> > > Hello John
> > >
> > > You inspired me and finally managed to find the
> optimized
> > > ellipsoid except the orientation. But you also
hinted
> that
> > > the optimization problem is performed, on the
> > > 6 parameters that make up H.
> > >
> > > I used
> > >
> > > H=diag(1./[r1^2,r2^2,r3^2]). Unfortunately, this
> matrix
> > > does not give the real orientation of the points I
> want to
> > > fit ellipsoids. It is only a matter of
orientation!!!!
> > >
> > > The rest is done.... Do you have an idea how I could
> > > include the orientation in the symetric matrix?
> > >
> > > Thanks to you and Burno
> >
> > Good! You are getting there. You might try
> > this:
> >
> > H = reshape(params([1 2 3 2 4 5 3 5 6]),3,3);
> >
> > This matrix will be symmetric, but it may have
> > negative eigenvalues. The problem is that you
> > need to constrain the eigenvalues to be
> > non-negative. So better would be to form H
> > as
> >
> > L = [params([1]) 0 0;params(2:3);params(4:6)];
> > H = L*L';
> >
> > Note that H is ALWAYS positive definite by this
> > definition, but also this gives you a direct
> > factorization of H.
> >
> > John
>

Toja -

I think the fitellipse package on the file exchange may be
of interest to you. It does not do exactly what you want,
but it contains a reference to a paper co-authored by the
late G. Golub with a discussion that I think will really
help.

http://www.mathworks.com/matlabcentral/files/15125/content/
demo/html/ellipsedemo.html

-Justin

Subject: How to specify the ellipsoids by equations alone?

From: John D'Errico

Date: 21 Nov, 2007 18:00:31

Message: 22 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message
<fi1ogd$e04$1@fred.mathworks.com>...
> Hello John,
>
> ok, but my problem was not making H positive definite
> matrix. Let me state it again:
>
> -I have a set of points, say 200
> - I am optimizing H in the least square distance
> difference between the points and the ellipsoid. That is
> searching the optimal radii r1,r2,r3. But the radii is not
> enough to define the ellipsoid. I want to determine the
> orientation information from my set of points and include
> in H,
>
> Do you have an idea?

Yes. Do as I said. Sometimes I even know what
I'm talking about. ;-) Yes, it is a rare occurrence.
It happened one other time before, in the 70's.

In a nutshell, you DO need a positive definite
matrix.

Your objective function for the optimization
must convert the 6 parameters that define
the ellipsoid matrix into a 3x3 matrix. Then
it must compute the errors that you wish to
minimize, and return the total error. Fmincon
will then vary the 6 parameters of that matrix,
optimizing the result. Yes? (Please nod your
head, agreeing with this statement.)

That 3x3 matrix MUST be positive definite
to define an ellipse. (Actually, non-negative
definite is probably adequate, in which case
you would have an ellipse with zero volume.)
If it is not at least non-negative definite, then
you may be finding the best approximation to
your data to a hyperboloid, instead of an ellipsoid.
There are two basic ways to solve this requirement.

1. We can use the ability of fmincon to handle
a nonlinear inequality constraint. The constraint
would be on the most negative eigenvalue of H,
that it must be greater than or equal to zero.

2. If we are clever in how we define H, then
by defining it as the product of two triangular
real matrices, H = L*L', then the matrix MUST
be non-negative definite by construction.

Option 2 is a FAR better choice here, and it is
what I recommended in my last post. Try it.

When fmincon has terminated, recreate the
matrix H from the factors H = L*L'. An eigenvalue
decomposition of H will yield the axis lengths
of the ellipsoid, and the corresponding
eigenvectors define the orientation of each
axis.

As I said, try it. I occasionally do get things
right.

John

Subject: How to specify the ellipsoids by equations alone?

From: Bruno Luong

Date: 22 Nov, 2007 21:36:29

Message: 23 of 26

Toja,

If your data (P2) are not very noisy, there is an idea that
is relatively simple to implement. It's certainly not the
best, far from it.

Let's call (E) the ellipsoide.

If X belong to (E), X verify

(X-X0)'*H*(X-X0) = 1

This can be see as
F(X) = 1, where F is a polynomial of second order, with out
constant term (H is the Hessian of F >0). The polynomial
coefficients {C} that define F depends on X0 and H.

If we solve
F(P2_n) = 1; n=1,...#data

in the least-square sense, with {C} are unknown; then the
problem is linear, and easy to implement (see code bellow).

The big inconvenience of this method is that it is not
robust when applied on noisy data. Better using non-linear
optimization tools for that.

Bruno

%
% Simulate data
%

ndim=2; % <- 2 or 3
L=rand(ndim); % lower triangle
H=L'*L % input H matrix, sym.def.pos.

X0=rand(ndim,1) % input center of the ellipse

sigma=0.02;

[V,D]=eig(H);
W=V*diag(1./sqrt(diag(D)));
ndata=100;
X=randn(ndim,ndata);
normX=sqrt(sum(X.^2,1));
X=X./repmat(normX,ndim,1);
X=X+sigma*randn(size(X)); % add random noise
P2=(repmat(X0,1,ndata)+(W*X))'; % Date points supposeed to
belong to the unknown ellipsoide

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Call solver

[Hest X0est] = solveellipse(P2)

% Plot data and ellipsoide

[V,D]=eig(Hest);
W=V*diag(1./sqrt(diag(D)));
npnt=200;
theta=[0:npnt]/npnt*2*pi;
X=[cos(theta); sin(theta)];
E=(repmat(X0est,1,size(X,2))+(W*X))'; % Ellipse estimated

figure;
if ndim==2
   hold on
   plot(P2(:,1),P2(:,2),'.');
   plot(E(:,1), E(:,2),'r');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Hest X0est] = solveellipse(X)
   
[ndata ndim]=size(X);

%
% Generate all combinaisons of polynomial-order <=2 for ndim
variables
%
order=cell(1,ndim);
order(:)={(0:2)};
ORDER=cell(1,ndim);
[ORDER{:}]=ndgrid(order{:});
ORDER=cellfun(@(grid) grid(:), ORDER, 'UniformOutput',0);
idx=find(sum(cell2mat(ORDER),2)<=2);
ORDER=cellfun(@(order) order(idx), ORDER, 'UniformOutput',0);
ORDER=cell2mat(ORDER);

%
% Remove constant term
%
ORDER=ORDER(2:end,:);

BIGX=repmat(X,[1 1 size(ORDER,1)]);
BIGORDER=repmat(reshape(ORDER',[1 size(ORDER,2)
size(ORDER,1)]),[ndata 1 1]);

M=squeeze(prod(BIGX.^BIGORDER,2));
P=M\ones(ndata,1);

I=repmat(eye(ndim),1,ndim);
I=mat2cell(I,ones(1,ndim),ndim*ones(1,ndim));
J=I';
K=cellfun(@(i,j) (i+j),I,J,'UniformOutput',0);

function loc=findloc(order) % nested function
[dummy loc]=ismember(order,ORDER,'rows');
end

P2_loc=cellfun(@findloc,K);
P2=P(P2_loc);
A=(tril(P2)+triu(P2))/2

I=mat2cell(eye(ndim),ones(1,ndim),ndim);
P1_loc=cellfun(@findloc,I);
P1=P(P1_loc);

X0est=-0.5*(A\P1);
lambda=1/(1+X0est'*A*X0est);
Hest=lambda*A;

end

Subject: How to specify the ellipsoids by equations alone?

From: Bruno Luong

Date: 25 Nov, 2007 21:51:36

Message: 24 of 26

Clean up my code, adding comments.

Bruno

function test_fitellipse

%
% Simulate data
%

ndim=2; % <- 2 or 3 or ... n
L=tril(randn(ndim)); % lower triangle
H=L'*L % input H matrix, sym.def.pos.

X0=randn(ndim,1) % input center of the ellipse

sigma=0.01;

[V,D]=eig(H);
W=V*diag(1./sqrt(diag(D)));
ndata=20;
X=randn(ndim,ndata);
normX=sqrt(sum(X.^2,1));
X=X./repmat(normX,ndim,1);
X=X+sigma*randn(size(X)); % add random noise
P2=(repmat(X0,1,ndata)+(W*X))'; % Date points supposeed to
belong to the unknown ellipsoide

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Call solver

[Hest X0est W] = solveellipse(P2)

% Plot data and ellipsoide

if ndim==2

    npnt=200;
    theta=[0:npnt]/npnt*2*pi;
    X=[cos(theta); sin(theta)];
    E=(repmat(X0est,1,size(X,2))+(W*X))'; % Ellipse estimated

    figure(1);
    clf(1);
    hold(gca,'on')
    plot(P2(:,1),P2(:,2),'.',...
         E(:,1), E(:,2),'r');
    axis(gca, 'equal');
end

end % of test_fitellipse

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [H X0 W] = solveellipse(X)
% function [H X0 W] = solveellipse(X);
%
% INPUT:
% X : (m x n) : m data points in R^n, they are supposed
to belong
% to an ellipsoide (or more generally a second-order
% implicite hyper-surface)
%
% OUTPUT:
% H : Matrix of the Bilinear form associated with the
ellipsoide
% Elippsoide = { X in R^n : (X-X0)' * H * (X-X0) = 1 }
% X0 : (n x 1), center of the ellipsoide
% W : (n x n) square matrix where each column is the axis
of ellipsoide
%

[ndata ndim]=size(X);

if ndim<2
    error('solveellipse: dimension number must be greater
than 1');
end

% This vector will be used in few places for reshapping purpose
uno = ones(1,ndim);

%
% Generate all combinations of polynomial-order <=2 for ndim
variables
%
order=cell(1,ndim);
order(:)={(0:2)};
ORDER=cell(1,ndim); % 1 x ndim cell of vectors (0:2)
[ORDER{:}]=ndgrid(order{:}); % Set {(0:2)}^ndim
ORDER=reshape(cat(ndim+1,ORDER{:}),[],ndim);
ORDER=ORDER(sum(ORDER,2)<=2,:); % second order only

%
% Remove constant term
%
ORDER=ORDER(2:end,:); % similar to ORDER(~any(ORDER,2),:)=[];
npol=size(ORDER,1); % number of basis 2nd order-polynomials

if ndata<npol
    error('solveellipse: Not enough data');
end

    function loc=findloc(order) % nested function
        %
        % Look for the location of the 'order' in the set of
{basis
        % 2nd-order polynomials}
        %
        [dummy loc]=ismember(reshape(order,1,[]),ORDER,'rows');
    end

%
% Both of these are 3d array of (ndata x ndim x npol)
%
BIGX=repmat(X,[1 1 npol]);
BIGORDER=repmat(reshape(ORDER',[1 ndim npol]),[ndata 1 1]);
M=squeeze(prod(BIGX.^BIGORDER,2)); % product of power in
every dimensions
clear BIGX BIGORDER

%
% Solve for polynomial coefficiens that lead to 1 on input
points
%
rhs=ones(ndata,1);
P=M\rhs;
clear M rhs;

%
% Extract the Hessian from the solution
%
I=repmat(eye(ndim),ndim,1);
I=reshape(I,ndim*[1 1 1]);
J=permute(I,[2 1 3]); % swap the first two dimensions
K=squeeze(mat2cell(I+J,uno,uno,ndim));
P2_loc=cellfun(@findloc,K);
P2=P(P2_loc); % second order term
A=(tril(P2)+triu(P2))/2; % devide by 2, accept diagonal terms

%
% Extract the gradient
%
I=mat2cell(eye(ndim),uno,ndim);
P1_loc=cellfun(@findloc,I);
P1=-0.5*P(P1_loc); % -0.5 * first order term

X0=A\P1;
lambda=1/(1+X0'*P1);
H=lambda*A;

if nargout>=3 % Compute main-axis of ellipsoide
    [V D]=eig(H);
    d=diag(D);
    if any(d<0)
        warning('solveellipse:npdH', ...
                'solveellipse: non positive Hessian matrix');
    end
    W = V.*repmat(1./sqrt(d'),ndim,1);
end

end

Subject: How to specify the ellipsoids by equations alone?

From: Toja Debela

Date: 26 Nov, 2007 10:22:26

Message: 25 of 26

Hello Bruno,

Your approach is being investigated. Hest can be used for
the starting point of the optimization. However, X0est is
very far...

I will let you know the comments I have soon.

Bruno,

How could i parametrise variable, X, to solve for the
intersection of two ellipsoids


(X-x1)*H1*(X-x1)=1
(X-x2)*H2*(X-x2)=1

any idea is welcome.


thanx,

toja





"Bruno Luong" <brunoluong@yahoo.com> wrote in message
<ficqp8$i5p$1@fred.mathworks.com>...
> Clean up my code, adding comments.
>
> Bruno
>
> function test_fitellipse
>
> %
> % Simulate data
> %
>
> ndim=2; % <- 2 or 3 or ... n
> L=tril(randn(ndim)); % lower triangle
> H=L'*L % input H matrix, sym.def.pos.
>
> X0=randn(ndim,1) % input center of the ellipse
>
> sigma=0.01;
>
> [V,D]=eig(H);
> W=V*diag(1./sqrt(diag(D)));
> ndata=20;
> X=randn(ndim,ndata);
> normX=sqrt(sum(X.^2,1));
> X=X./repmat(normX,ndim,1);
> X=X+sigma*randn(size(X)); % add random noise
> P2=(repmat(X0,1,ndata)+(W*X))'; % Date points supposeed to
> belong to the unknown ellipsoide
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> % Call solver
>
> [Hest X0est W] = solveellipse(P2)
>
> % Plot data and ellipsoide
>
> if ndim==2
>
> npnt=200;
> theta=[0:npnt]/npnt*2*pi;
> X=[cos(theta); sin(theta)];
> E=(repmat(X0est,1,size(X,2))+(W*X))'; % Ellipse
estimated
>
> figure(1);
> clf(1);
> hold(gca,'on')
> plot(P2(:,1),P2(:,2),'.',...
> E(:,1), E(:,2),'r');
> axis(gca, 'equal');
> end
>
> end % of test_fitellipse
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> function [H X0 W] = solveellipse(X)
> % function [H X0 W] = solveellipse(X);
> %
> % INPUT:
> % X : (m x n) : m data points in R^n, they are supposed
> to belong
> % to an ellipsoide (or more generally a second-
order
> % implicite hyper-surface)
> %
> % OUTPUT:
> % H : Matrix of the Bilinear form associated with the
> ellipsoide
> % Elippsoide = { X in R^n : (X-X0)' * H * (X-X0) =
1 }
> % X0 : (n x 1), center of the ellipsoide
> % W : (n x n) square matrix where each column is the
axis
> of ellipsoide
> %
>
> [ndata ndim]=size(X);
>
> if ndim<2
> error('solveellipse: dimension number must be greater
> than 1');
> end
>
> % This vector will be used in few places for reshapping
purpose
> uno = ones(1,ndim);
>
> %
> % Generate all combinations of polynomial-order <=2 for
ndim
> variables
> %
> order=cell(1,ndim);
> order(:)={(0:2)};
> ORDER=cell(1,ndim); % 1 x ndim cell of vectors (0:2)
> [ORDER{:}]=ndgrid(order{:}); % Set {(0:2)}^ndim
> ORDER=reshape(cat(ndim+1,ORDER{:}),[],ndim);
> ORDER=ORDER(sum(ORDER,2)<=2,:); % second order only
>
> %
> % Remove constant term
> %
> ORDER=ORDER(2:end,:); % similar to ORDER(~any(ORDER,2),:)=
[];
> npol=size(ORDER,1); % number of basis 2nd order-
polynomials
>
> if ndata<npol
> error('solveellipse: Not enough data');
> end
>
> function loc=findloc(order) % nested function
> %
> % Look for the location of the 'order' in the set
of
> {basis
> % 2nd-order polynomials}
> %
> [dummy loc]=ismember(reshape(order,1,
[]),ORDER,'rows');
> end
>
> %
> % Both of these are 3d array of (ndata x ndim x npol)
> %
> BIGX=repmat(X,[1 1 npol]);
> BIGORDER=repmat(reshape(ORDER',[1 ndim npol]),[ndata 1
1]);
> M=squeeze(prod(BIGX.^BIGORDER,2)); % product of power in
> every dimensions
> clear BIGX BIGORDER
>
> %
> % Solve for polynomial coefficiens that lead to 1 on input
> points
> %
> rhs=ones(ndata,1);
> P=M\rhs;
> clear M rhs;
>
> %
> % Extract the Hessian from the solution
> %
> I=repmat(eye(ndim),ndim,1);
> I=reshape(I,ndim*[1 1 1]);
> J=permute(I,[2 1 3]); % swap the first two dimensions
> K=squeeze(mat2cell(I+J,uno,uno,ndim));
> P2_loc=cellfun(@findloc,K);
> P2=P(P2_loc); % second order term
> A=(tril(P2)+triu(P2))/2; % devide by 2, accept diagonal
terms
>
> %
> % Extract the gradient
> %
> I=mat2cell(eye(ndim),uno,ndim);
> P1_loc=cellfun(@findloc,I);
> P1=-0.5*P(P1_loc); % -0.5 * first order term
>
> X0=A\P1;
> lambda=1/(1+X0'*P1);
> H=lambda*A;
>
> if nargout>=3 % Compute main-axis of ellipsoide
> [V D]=eig(H);
> d=diag(D);
> if any(d<0)
> warning('solveellipse:npdH', ...
> 'solveellipse: non positive Hessian
matrix');
> end
> W = V.*repmat(1./sqrt(d'),ndim,1);
> end
>
> end
>
>
>
>

Subject: How to specify the ellipsoids by equations alone?

From: Bruno Luong

Date: 26 Nov, 2007 12:15:46

Message: 26 of 26

"Toja Debela" <tojadeb@mathworks.com> wrote in message
<fie6p2$sn5$1@ginger.mathworks.com>...

>
> How could i parametrise variable, X, to solve for the
> intersection of two ellipsoids
>
>
> (X-x1)*H1*(X-x1)=1
> (X-x2)*H2*(X-x2)=1
>

Find the main axis W2 of H2 (e.g., in my code), and use the
change of coordinate as following

X -> Y := W2^(-1) * (X-x2)

therefore
Y -> X = x2 + W2*Y

X belong to (E2), if and only if Y belong to a unit sphere.

Use azimuth and elevation angles (similar to polar
coordinates) to parametrize unitsphere for Y.

Solve for angles so that X belong to (E1) means solving:

(W2*Y + x2 - x1)' * H1 * (W2*Y + x2 - x1) = 1 .

This is however a non-linear problem.

Bruno

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us