Message Boards Message Boards

Fit ellipse to a cloud of points

Posted 10 years ago

I got this, but why not fit an ellipse to all points?

Off[NMinimize::eit]; Off[NMinimize::lstol];
DynamicModule[{rng, ptk, szumx, szumy, final}, 
 Panel[Column[{"Theta z przedzia?u od 0 do:", 
    RadioButtonBar[
     Dynamic[rng], {2 Pi -> "2\[Pi]", Pi -> "\[Pi]", 
      Pi/2 -> "\[Pi]/2", Pi/8 -> "\[Pi]/8"}], , 
    "Ilo?? generowanych punktów:", 
    RadioButtonBar[
     Dynamic[ptk], {5 -> "5", 20 -> "20", 50 -> "50", 
      100 -> "100"}], , "Warto?ci wektora przesuni?cia", 
    "Warto?? xi:" RadioButtonBar[
      Dynamic[szumx], {-0.1 -> "-0.1", -.5 -> "-0.5", -.9 -> 
        "-0.9", -3 -> "-3"}], 
    "Warto?? yj:" RadioButtonBar[
      Dynamic[szumy], {0.1 -> "0.1", .5 -> "0.5", .9 -> "0.9", 
       3 -> "3"}], 
    Button["Generuj Elipse", 
     While[True,(*Randomly choose coefficients until acceptable*){a, 
        b, c, d, f, g} = RandomReal[{-10, 10}, 6];
      \[CapitalDelta] = -c d^2 + 2 b d f - a f^2 - b^2 g + a c g;
      j = -b^2 + a c; i = a + c;
      If[\[CapitalDelta] != 0 && j > 0 && \[CapitalDelta]/i < 0, 
       Break[]]];
     ellipse = a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g;
     (*Center of an ellipse in general form is{(c d-b f)/(b^2-
     a c),(a f-b d)/(b^2-a c)}*)
     points = Table[theta = RandomReal[{0, rng}];
       ksol = 
        FindRoot[(ellipse /. {x -> 
              k*Cos[theta] + (c d - b f)/(b^2 - a c), 
             y -> k*Sin[theta] + (a f - b d)/(b^2 - a c)}) == 0, {k, 
          1.}];
       Point[{x, y}] /. {x -> k*Cos[theta] + (c d - b f)/(b^2 - a c), 
          y -> k*Sin[theta] + (a f - b d)/(b^2 - a c)} /. ksol, {ptk}];
     nearpoints = 
      points /. 
       Point[{x_, y_}] :> 
        Point[{x + RandomReal[{szumx, szumy}], 
          y + RandomReal[{szumx, szumy}]}];
     (*ellipse x and y min and max values*)
     yplotrange = 
      Flatten[{y, 
        Sort[{(2*b*d - 2*a*f + 
             Sqrt[(2*b*d - 2*a*f)^2 - 
               4*(b^2 - a*c)*(d^2 - a*g)])/(2*(-b^2 + a*c)), (-2*b*
              d + 2*a*f + 
             Sqrt[(2*b*d - 2*a*f)^2 - 
               4*(b^2 - a*c)*(d^2 - a*g)])/(2*(b^2 - a*c))}]}];
     xplotrange = 
      Flatten[{x, 
        Sort[{(2*c*d - 2*b*f + 

             Sqrt[(-2*c*d + 2*b*f)^2 - 
               4*(b^2 - a*c)*(f^2 - c*g)])/(2*(b^2 - a*c)), (-2*c*d + 
             2*b*f + 
             Sqrt[(-2*c*d + 2*b*f)^2 - 
               4*(b^2 - a*c)*(f^2 - c*g)])/(2*(-b^2 + a*c))}]}];
     (*minimize distance of near points to a new general ellipse*)

     {xs, ys} = Transpose[nearCoords];
     newellipse = aa*x^2 + 2*bb*x*y + cc*y^2 + 2*dd*x + 2*ff*y + gg;
     distance = Plus @@ (newellipse^2 /. {x -> xs, y -> ys});
     {res, coes} = 
      NMinimize[{distance, bb^2 - 4 aa*cc <= 0}, {aa, bb, cc, dd, ff, 
        gg}];
     scaleup = FromDigits[{{1}, Last@RealDigits[1/(gg /. coes)] + 1}];
      esolve = Expand[scaleup*(newellipse /. coes)];
     (*
     {xs,ys}=Transpose[nearCoords];
     newellipse=aa*x^2+2*bb*x*y+cc*y^2+2*dd*x+2*ff*y+gg;
     distance=Plus@@(newellipse^2/.{x->xs,y->ys});
     {res,coes}=NMinimize[distance,{aa,bb,cc,dd,ff,gg}];
     scaleup=FromDigits[{{1},Last@RealDigits[1/(gg/.coes)]+1}];
     esolve=Expand[scaleup*(newellipse/.coes)];
     *)
     final = 
      Show[ContourPlot[{ellipse == 0, esolve == 0}, 
        Evaluate[xplotrange], Evaluate[yplotrange]], Graphics[points],
        Graphics[{Red, nearpoints}], ImageSize -> {500, Automatic}];],
     Dynamic[final]}]]]

enter image description here

I do not know if I gave a good KKT condition to form an ellipse always in here

{xs, ys} = Transpose[nearCoords];
newellipse = aa*x^2 + 2*bb*x*y + cc*y^2 + 2*dd*x + 2*ff*y + gg;
distance = Plus @@ (newellipse^2 /. {x -> xs, y -> ys});
{res, coes} = 
  NMinimize[{distance, bb^2 - 4 aa*cc <= 0}, {aa, bb, cc, dd, ff, gg}];
scaleup = 
 FromDigits[{{1}, Last@RealDigits[1/(gg /. coes)] + 1}]; esolve = 
 Expand[scaleup*(newellipse /. coes)];
POSTED BY: Radek Drozd
26 Replies
Posted 10 years ago

Ok i did now. Why f must be negative?

POSTED BY: Radek Drozd

Did you clear all the variables?

POSTED BY: Frank Kampas
Posted 10 years ago

what i'm doing wrong?

In[206]:= data = Table[{x, 1/x}, {x, Range[6]}];
val = (a x^2 + b*x*y + c y^2 + d*x + e*y + f)^2 /. 
   Table[Thread[{x, y} -> data[[i]]], {i, Length[data]}];
NSolve[Thread[
  D[Plus @@ val + \[Lambda] (b^2 - 4*a*c + 1), {{a, b, c, d, e, 
      f, \[Lambda]}}] == 0], {a, b, c, d, e, f, \[Lambda]}, Reals]

During evaluation of In[206]:= General::ivar: 7.661932292992361` is not a valid variable. >>

During evaluation of In[206]:= NSolve::ivar: 7.661932292992361` is not a valid variable. >>

Out[208]= NSolve[\!\(
\*SubscriptBox[\(\[PartialD]\), \({{7.661932292992361`, 
      5.640351135084579`, 
      8.174655653861642`, \(-4.7774996460771035`\), \
\(-8.349009029331341`\), \(-7.657080239171673`\), \
\[Lambda]}}\)]\((99361.89451388319`  - 
     217.7210716266105`\ \[Lambda])\)\) == 0, {7.66193, 5.64035, 
  8.17466, -4.7775, -8.34901, -7.65708, \[Lambda]}, Reals]
POSTED BY: Radek Drozd
data = Table[{x, 1/x}, {x, Range[6]}];

val = (a x^2 + b*x*y + c y^2 + d*x + e*y + f)^2 /. 
   Table[Thread[{x, y} -> data[[i]]], {i, Length[data]}];

In[16]:= NSolve[
 Thread[D[Plus @@ val + \[Lambda] (b^2 - 4*a*c + 1), {{a, b, c, d, e, 
      f, \[Lambda]}}] == 0], {a, b, c, d, e, f, \[Lambda]}, Reals]

Out[16]= {{a -> 0.071953, b -> 0, c -> 3.47449, d -> -0.922589, 
  e -> -6.58444, 
  f -> 3.96278, \[Lambda] -> 0.00264096}, {a -> -0.071953, b -> 0, 
  c -> -3.47449, d -> 0.922589, e -> 6.58444, 
  f -> -3.96278, \[Lambda] -> 0.00264096}}

You want the solution where f is negative.

POSTED BY: Frank Kampas
Posted 10 years ago

I have version 8. Yes now it works.

POSTED BY: Radek Drozd

I think Grad is a new function in Version 10. Try

In[2]:= NSolve[
 Thread[D[x + y + \[Lambda] (x^2 + y^2 - 1), {{x, y, \[Lambda]}}] == 
   0], {x, y, \[Lambda]}]

Out[2]= {{x -> 0.707107, 
  y -> 0.707107, \[Lambda] -> -0.707107}, {x -> -0.707107, 
  y -> -0.707107, \[Lambda] -> 0.707107}}
POSTED BY: Frank Kampas
Posted 10 years ago

i try simulate this in mathematica and this came

In[3]:= NSolve[
 Thread[Grad[x + y + \[Lambda] (x^2 + y^2 - 1), {x, y, \[Lambda]}] == 
   0], {x, y, \[Lambda]}]

During evaluation of In[3]:= NSolve::nsmet: This system cannot be solved with the methods available to NSolve. >>
POSTED BY: Radek Drozd

Here's an example of finding the minimum and maximum points for an objective function = x + y and a single equality constraint x^2 + y^2=1

In[42]:= NSolve[
 Thread[Grad[x + y + \[Lambda] (x^2 + y^2 - 1), {x, y, \[Lambda]}] == 
   0], {x, y, \[Lambda]}]

Out[42]= {{x -> 0.707107, 
  y -> 0.707107, \[Lambda] -> -0.707107}, {x -> -0.707107, 
  y -> -0.707107, \[Lambda] -> 0.707107}}
POSTED BY: Frank Kampas

Well, this really isn't a site for having people write code for your projects or homework. So I'll outline it. Your problem is to minimize a sum of squares involving the quadratic parameter variables {a,b,c,d,e,f} where the quadratic is given by a*x^2+b*x*y+c*y^2+d*x+e*y+f==0. Your constraint is b^2-4*a*c+1==0. Define the function g[a,b,c]=b^2-4*a*c+1. So your objective function (the sum of squares), call it obj[a,b,c] (we ignore that it is a function also of variables {d,e,f}), must satisfy grad(obj) == lambda*grad(g) (where I use grad(...) to denote "take the gradient of ... with respect to the stated set of variables").

You can get seven equations in the seven variables {a,b,c,d,e,f,lambda} from:

g==0 (one eqn) grad obj == lambda*grad g (three eqns, since this we equate vectors of three components) D[obj,d]==0 D[obj,e]==0 D[obj,f]==0

These will be polynomial equations so one can use NSolve to solve the system.

This is about as complete an answer as I am inclined to give. More than that and I'd be doing a substantial part of your project.

POSTED BY: Daniel Lichtblau

The KKT equations extend Lagrange multipliers to inequality constraints by adding a requirement that the multiplier for an inequality constraint be zero if the constraint is not active. This is done by requiring that the product of the multiplier and the constraint value be zero. Therefore using the constraint b^2-4ac=-1 with a Lagrange multiplier is equivalent to the KKT equations, in the special case that there are no inequality constraints.

POSTED BY: Frank Kampas
Posted 10 years ago

And how to introduce it in mathematica? Because i think it's good idea.

POSTED BY: Radek Drozd

My advice is to use a Lagrange multiplier and call it a KKT multiplier. Nobody will be any the wiser. The Lagrange multiplier, used as I noted, will give ellipses.

POSTED BY: Daniel Lichtblau
Posted 10 years ago

Yes, I've read. Whether these two methods can combine with KKT? Because i need that.

POSTED BY: Radek Drozd

"i need generate ellipse alaways from cloud of points"

I am aware of that. Did you read my response from 10 or so hours prior to yours? It mentions two ways for enforcing ellipses.

POSTED BY: Daniel Lichtblau
Posted 10 years ago

and whether the KKT handle with this issue,

enter image description here

it is my application which work without KKT and red elipse trying for matching from clouds of points just do not always get out the ellipse enter image description here

i need generate ellipse alaways from cloud of points

POSTED BY: Radek Drozd

I fail to understand this requirement that KKT multipliers be used. Enforcing the needed nonnegativity condition will be difficult or impossible. Since you want an ellipse and your quadratic parameters can be rescaled, it makes more sense just to insist on the equality b^2-4ac=-1. This can be done with a Lagrange multiplier.

A generally better way to fit ellipses to points is from A. Fitzgibbon, M. Pilu, and R. Fisher. "Direct least square fitting of ellipses". IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(5), 476--480, May 1999. Mathematica implementation code can be fond at the link below. http://forums.wolfram.com/mathgroup/archive/2001/Sep/msg00179.html

POSTED BY: Daniel Lichtblau
Posted 10 years ago

want to make my application was intact, and add in the menu KKT. That's the part I have to take in the menu?

 nearCoords = nearpoints[[All, 1]];
     {xs, ys} = Transpose[nearCoords];
     newellipse = aa*x^2 + 2*bb*x*y + cc*y^2 + 2*dd*x + 2*ff*y + gg;
     distance = Plus @@ (newellipse^2 /. {x -> xs, y -> ys});
     {res, coes} = NMinimize[distance, {aa, bb, cc, dd, ff, gg}];
     scaleup = FromDigits[{{1}, Last@RealDigits[1/(gg /. coes)] + 1}];
     esolve = Expand[scaleup*(newellipse /. coes)];
POSTED BY: Radek Drozd
Posted 10 years ago

I understand it. My professor did not care about the code, but about adding KKT and of course must work to create a few examples. For me it is a little difficult. I understand my application but your code not so much.

POSTED BY: Radek Drozd

I suspect your professor is not going to be satisfied with your using the code I wrote. For one thing, I wrote the code to be as compact as possible, not as understandable as possible. I also wrote it so that the subscripts on the multipliers are the constraints they refer to. This is something that is easily done in Mathematica but very nonstandard. Since there is only one constraint, constructing the KKT equations is straightforward. I suggest you do that yourself for your code.

POSTED BY: Frank Kampas
Posted 10 years ago

ie, I will have to insert it instead of this, yes?

  nearCoords = nearpoints[[All, 1]];
     {xs, ys} = Transpose[nearCoords];
     newellipse = aa*x^2 + 2*bb*x*y + cc*y^2 + 2*dd*x + 2*ff*y + gg;
     distance = Plus @@ (newellipse^2 /. {x -> xs, y -> ys});
     {res, coes} = NMinimize[distance, {aa, bb, cc, dd, ff, gg}];
     scaleup = FromDigits[{{1}, Last@RealDigits[1/(gg /. coes)] + 1}];
     esolve = Expand[scaleup*(newellipse /. coes)];
POSTED BY: Radek Drozd

I modified my code to use NSolve instead of Reduce for this application.

data = Table[{x, 1/x}, {x, Range[6]}];
val = (a x^2 + b*x*y + c y^2 + d*x + e*y + f)^2 /. 
   Table[Thread[{x, y} -> data[[i]]], {i, Length[data]}];

consconvrule = {x_ >= y_ -> y - x, x_ == y_ -> x - y, 
   x_ <= y_ -> x - y, lb_ <= x_ <= ub_ -> (x - lb) (x - ub), 
   ub_ >= x_ >= lb_ -> (x - lb) (x - ub)};

KKTN[obj_, cons_List, vars_List] := 
 Module[{convcons, gradobj, gradcons, \[CapitalLambda]}, 
  convcons = (cons /. consconvrule);
  {gradobj, gradcons} = D[{obj, convcons}, {vars}];
  \[CapitalLambda] = Subscript[\[Lambda], #] & /@ cons;
  NSolve[Flatten[{Thread[gradobj == \[CapitalLambda].gradcons], 
     Thread[\[CapitalLambda]*convcons == 0] /. 
      Subscript[\[Lambda], Equal[a_, b_]] -> 0, 
     cons, {objval == obj}}], Join[{objval}, vars, \[CapitalLambda]], 
   Reals]]

resn = KKTN[Plus @@ val, {b^2 - 4*a*c <= -10^-3}, {a, b, c, d, e, f}]
{{a -> -0.00227535, b -> 0, c -> -0.109873, d -> 0.0291748, 
  e -> 0.208218, f -> -0.125314, 
  Subscript[\[Lambda], b^2 - 4 a c <= -(1/1000)] -> -0.00264096, 
  objval -> 2.64096*10^-6}}

Show[ContourPlot[
  Evaluate[a x^2 + b*x*y + c y^2 + d*x + e*y + f /. resn[[1]]] == 
   0, {x, 0, 16}, {y, -4, 4}], ListPlot[data]] gives an ellipse which passes through the points.
POSTED BY: Frank Kampas
Posted 10 years ago

because it requires me to do my professor

POSTED BY: Radek Drozd

Why do you want to use the KKT conditions when it's simpler to use NMinimize with the constraint b^2- 4ac <= - epsilon, where epsilon is a small number, like 10^-3 ?

POSTED BY: Frank Kampas
Posted 10 years ago

I know that it was your code. Here is not successful this add and resolve? So now I do not know how to do it in order to add KKT to my application with the condition that the ellipse always emerged from the cloud of points

POSTED BY: Radek Drozd

That's code I wrote to construct the Kuhn-Karush-Tucker conditions and solve them with Reduce. It's somewhat nonstandard as it finds maxima and minima simultaneously. It's designed for small problems which can be handled symbolically.

POSTED BY: Frank Kampas
Posted 10 years ago

if someone help me to understand this I managed to add it to my application?

 consconvrule = {
    x_ >= y_ -> y - x,
    x_ == y_ -> x - y ,
    x_ <= y_ -> x - y,
    lb_ <= x_ <= ub_ -> (x - lb) (x - ub),
    ub_ >= x_ >= lb_ -> (x - lb) (x - ub)
    };

 KKT[obj_, cons_List, vars_List] :=
 Module[
  {convcons, gradobj, gradcons, ?},
  convcons = (cons /. consconvrule);
  {gradobj, gradcons} = D[{obj, convcons}, {vars}];
  ? = Subscript[?, #] & /@ cons;
  LogicalExpand @ Reduce[
    Flatten[{
      Thread[gradobj == ?.gradcons],
      Thread[?*convcons == 0] /.
       Subscript[?, Equal[a_, b_]] -> 0,
      cons,
      {objval == obj}
      }],
    Join[{objval}, vars, ?],
    Reals,
    Backsubstitution -> True
    ]
  ]

to this application

Off[NMinimize::eit]; Off[NMinimize::lstol];
DynamicModule[{rng, ptk, szumx, szumy, final}, 
 Panel[Column[{"Theta z przedzia?u od 0 do:", 
    RadioButtonBar[
     Dynamic[rng], {2 Pi -> "2\[Pi]", Pi -> "\[Pi]", 
      Pi/2 -> "\[Pi]/2", Pi/8 -> "\[Pi]/8"}], , 
    "Ilo?? generowanych punktów:", 
    RadioButtonBar[
     Dynamic[ptk], {5 -> "5", 20 -> "20", 50 -> "50", 
      100 -> "100"}], , "Warto?ci wektora przesuni?cia", 
    "Warto?? xi:" RadioButtonBar[
      Dynamic[szumx], {-0.1 -> "-0.1", -.5 -> "-0.5", -.9 -> 
        "-0.9", -3 -> "-3"}], 
    "Warto?? yj:" RadioButtonBar[
      Dynamic[szumy], {0.1 -> "0.1", .5 -> "0.5", .9 -> "0.9", 
       3 -> "3"}], 
    Button["Generuj Elipse", 
     While[True,(*Randomly choose coefficients until acceptable*){a, 
        b, c, d, f, g} = RandomReal[{-10, 10}, 6];
      \[CapitalDelta] = -c d^2 + 2 b d f - a f^2 - b^2 g + a c g;
      j = -b^2 + a c; i = a + c;
      If[\[CapitalDelta] != 0 && j > 0 && \[CapitalDelta]/i < 0, 
       Break[]]];
     ellipse = a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g;
     (*Center of an ellipse in general form is{(c d-b f)/(b^2-
     a c),(a f-b d)/(b^2-a c)}*)
     points = Table[theta = RandomReal[{0, rng}];
       ksol = 
        FindRoot[(ellipse /. {x -> 
              k*Cos[theta] + (c d - b f)/(b^2 - a c), 
             y -> k*Sin[theta] + (a f - b d)/(b^2 - a c)}) == 0, {k, 
          1.}];
       Point[{x, y}] /. {x -> k*Cos[theta] + (c d - b f)/(b^2 - a c), 
          y -> k*Sin[theta] + (a f - b d)/(b^2 - a c)} /. ksol, {ptk}];
     nearpoints = 
      points /. 
       Point[{x_, y_}] :> 
        Point[{x + RandomReal[{szumx, szumy}], 
          y + RandomReal[{szumx, szumy}]}];
     (*ellipse x and y min and max values*)
     yplotrange = 
      Flatten[{y, 
        Sort[{(2*b*d - 2*a*f + 
             Sqrt[(2*b*d - 2*a*f)^2 - 
               4*(b^2 - a*c)*(d^2 - a*g)])/(2*(-b^2 + a*c)), (-2*b*
              d + 2*a*f + 
             Sqrt[(2*b*d - 2*a*f)^2 - 
               4*(b^2 - a*c)*(d^2 - a*g)])/(2*(b^2 - a*c))}]}];
     xplotrange = 
      Flatten[{x, 
        Sort[{(2*c*d - 2*b*f + 
             Sqrt[(-2*c*d + 2*b*f)^2 - 
               4*(b^2 - a*c)*(f^2 - c*g)])/(2*(b^2 - a*c)), (-2*c*d + 
             2*b*f + 
             Sqrt[(-2*c*d + 2*b*f)^2 - 
               4*(b^2 - a*c)*(f^2 - c*g)])/(2*(-b^2 + a*c))}]}];
     (*minimize distance of near points to a new general ellipse*)
     nearCoords = nearpoints[[All, 1]];
     {xs, ys} = Transpose[nearCoords];
     newellipse = aa*x^2 + 2*bb*x*y + cc*y^2 + 2*dd*x + 2*ff*y + gg;
     distance = Plus @@ (newellipse^2 /. {x -> xs, y -> ys});
     {res, coes} = NMinimize[distance, {aa, bb, cc, dd, ff, gg}];
     scaleup = FromDigits[{{1}, Last@RealDigits[1/(gg /. coes)] + 1}];
     esolve = Expand[scaleup*(newellipse /. coes)];
     final = 
      Show[ContourPlot[{ellipse == 0, esolve == 0}, 
        Evaluate[xplotrange], Evaluate[yplotrange]], Graphics[points],
        Graphics[{Red, nearpoints}], ImageSize -> {500, Automatic}];],
     Dynamic@final
    (*,Dynamic@Grid[Join[{{"points","nearby points"}},
    Transpose[{points,nearpoints}][[All,All,1]]],Frame->All,
    Alignment->Left]*)}]]]
POSTED BY: Radek Drozd

Group Abstract Group Abstract