       Re: MATLAB to Mathematica 4.1

• To: mathgroup at smc.vnet.net
• Subject: [mg30790] Re: [mg30694] MATLAB to Mathematica 4.1
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 19 Sep 2001 00:16:29 -0400 (EDT)
• References: <200109080622.CAA25718@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```Moranresearch at aol.com wrote:
>
> Could someone please help translate this six line of MATLAB code into Mathematica? I would like to use this method to fit ellipses to experimental data.
> Source code for the algorithm is supplied and a beautiful demonstration is available on http://vision.dai.ed.ac.uk/maurizp/ElliFitDemo/demo.htmlThank you. John R Moran PhD, MD
>
> % x,y are list of coordinates
> function a = fit_ellipse(x,y)
>
> % Build design matrix
> D = [x.*x x.*y y.*y x y ones(size(x))];
>
> % Build scatter matrix
> S = D'*D;
>
> % Build 6x6 constarint matrix
> C(6,6) = 0; C(1,3) = 2; C(2,2) = -1; C(3,1) = 2;
>
> % Solve eigensystem
> [gevec,geval] = eig(inv(S)**C);
>
> % Find positive eigenvalue
> [PosR,PosC] = find(geval > 0 & ~ isinf(geval));
>
> % Extract eigenvector corresponding to positive eigenvalue
> a = gevec(:,PosC);
>
> Source code for the algorithm is supplied and a beautiful demonstration is available on http://vision.dai.ed.ac.uk/maurizp/ElliFitDemo/demo.html

It can be done as follows.

conicParameters[pts_] := Map[
{#[]^2,#[]*#[],#[]^2,#[],#[],1}&, pts]

scatterMat[pts_] := With[{sixtuple=conicParameters[pts]},
Transpose[sixtuple].sixtuple]
constrainedMat = Table[0, {6}, {6}];
constrainedMat[[1,3]] = constrainedMat[[3,1]] = 2; constrainedMat[[2,2]]
= -1;

ellipseFit[data_, {x_,y_}] := Module[
{gevals, gevecs, smat=scatterMat[data], posval, soln},
{gevals, gevecs} = Eigensystem[Inverse[smat] . constrainedMat];
posval = Ordering[gevals, -1];
gevecs[[posval[]]] . {x^2,x*y,y^2,x,y,1}
]

Here is a simple example. I'll start with several random x-values in the
interval from -1 to 1, then compute y-values to place them on the upper
part of the unit circle. We then take a random rotate-and-sheer (just a
random 2x2 matrix), and a random displacement in 2D. This gives us
points on some ellipse. Finally we add noise so they do not lie quite on
an ellipse anymore.

ellipseMat = Table[Random[], {2}, {2}]
translation = Table[Random[], {2}]
initXVals = Table[2*Random[]-1, {20}]
unitCircPts = Map[{#,Sqrt[1-#^2]}&, initXVals]
ellipsePts =  Map[ellipseMat.#+translation&, unitCircPts]
noisyEllipsePts = Map[#+.01*{Random[],Random[]}&, ellipsePts]

Now we compute the Fitzgibbon-Pilu-Fisher best fitting ellipse for the
particular run I did. Results are of course dependent on the random
values used in a given test.

In:= InputForm[ellipse = ellipseFit[noisyEllipsePts, {x,y}]]
Out//InputForm=
-0.018554479398057374 + 0.46270036775190626*x - 0.6208625862435061*x^2 -
0.17771140084160036*y + 0.5435317417625108*x*y - 0.2703435845125049*y^2

The remaining code below is to give a colorful plot of points and
computed ellipse.

colorPoint[pt_] := {RGBColor[Random[],Random[],Random[]], pt}

pts = ListPlot[noisyEllipsePts, DisplayFunction->Identity];
pts = (pts /. pt_Point :> colorPoint[pt]) /.
Graphics[{a___},b_] :> Graphics[{PointSize[.02],a},b];

p1 = ContourPlot[ellipse, {x,-1,1}, {y,-1,1},
DisplayFunction->Identity];

This will give a nice picture.

Show[Graphics[p1] /. GrayLevel[_]->RGBColor[.1,.4,1], pts,
DisplayFunction->\$DisplayFunction]

This is a good example of constrained optimization where the constraint
as well as objective function is quadratic. Which means that the
Lagrange multiplier lambda becomes a (generalized) eigenvalue lambda. I
guess this demonstrates some sort of internal self-consistency of the
Greek alphabet.

Daniel Lichtblau
Wolfram Research

```

• Prev by Date: RE: Combinations
• Next by Date: RE: Re: Summing list subsets
• Previous by thread: MATLAB to Mathematica 4.1
• Next by thread: Extract coefs. and exponents from a list . Package problem