Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

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[
	{#[[1]]^2,#[[1]]*#[[2]],#[[2]]^2,#[[1]],#[[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[[1]]]] . {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[161]:= InputForm[ellipse = ellipseFit[noisyEllipsePts, {x,y}]]
Out[161]//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},
	Contours->{0}, ContourShading->False, PlotPoints->50,
	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

• References:
  • MATLAB to Mathematica 4.1
    • From: Moranresearch@aol.com