8

The data is:

data = {{13, 60}, {36, 35}, {70, 29}, {106, 41}, {118, 72}, {94, 
    94}, {49, 99}, {24, 85}};

When I try to find a fit circle,So my method is work well:

circle = Circle[{x, y}, r];
fitcir = RegionDistance[circle, #] & /@ data // Total // 
  ArgMin[# // N, {x, y, r}] &

(*{68.8023, 63.6189, 43.5349}*)

The effect like:

enter image description here

But the same method for a ellipse will be failure:

ellipse = Circle[{x, y}, {a, b}];
RegionDistance[ellipse , #] & /@ data // Total // 
 ArgMin[# // N, {x, y, a, b}] &

The code will collapse my Mathematica,and furthermore,

RegionDistance[ellipse,{13, 60}]

will be failure to run.

But we can see

In[36]:= ellipse // RegionQ
Out[36]= True

So it is a bug for RegionDistance?

yode
  • 26,686
  • 4
  • 62
  • 167
  • 3
    Please wait for confirmation from others who observe this behaviour before adding the "bugs" tag – dr.blochwave Oct 14 '15 at 07:48
  • @blochwave Oh,I get it.thanks for your prompt. – yode Oct 14 '15 at 07:57
  • 2
    with the above fix under v10.1 it just runs for a long time. I think that's just a consequence of the complexity of the calculation of an analytic distance between a point and a general ellipse. (ie. not a bug IMO ). – george2079 Oct 14 '15 at 14:16
  • @george2079 I'm very sorry to make this typo,and I have corrected it just. – yode Oct 14 '15 at 16:26

2 Answers2

8

I think the most significant speedup can be achieved using machine precision explicitly in the calculation of the distances, rather than arbitrary precision numbers. This is the reason I apply N to data before calculating the distances below.


A tweak is possible in the Nelder-Mead minimization algorithm used by NMinimize (I tested other algorithms and they all seemed slower). It seems effective to provide a reasonable initial guess of the ellipse's parameters. An acceptable estimate can be obtained using the bivariate mean and standard deviation of data. The center of the ellipse is approximated by the Mean of data; the length of the semiaxes is estimated using the standard deviation of data. This of course is inaccurate, but it is good enough to provide a reasonable starting point for the minimizer.

It is also interesting to the note that the post-processing steps (local minimization) taken after the Nelder-Mead algorithm seem crucial to obtaining a good fit in these cases. It is easy to convince oneself of this by turning those off using "PostProcess" -> False as a further option to the Nelder-Mead algorithm. The much slower post-processing is responsible for the slow-down of the counters observed during optimization, and for the final cleanup of the obtained ellipse that appears as a sudden jump in the animated series of intermediate steps shown at the very bottom of this answer.


I then instrumented the minimization code to provide step (s) and evaluation (e) counters updated dynamically during the calculation, so one knows what is going on, and by harvesting the intermediate results using Sow and Reap. The latter was just for fun, to see what stages the optimization traversed.


Here are the code and results of the minimization:

data = {{13, 60}, {36, 35}, {70, 29}, {106, 41}, {118, 72}, {94, 94}, {49, 99}, {24, 85}};

s = 0; e = 0;
Row[{"Evaluation: ", Dynamic[e], ", Steps: ", Dynamic[s]}]

(* Notice application of N to data before calculating distances *)
totdist = Total[RegionDistance[Circle[{x, y}, {a, b}], #] & /@ N[data]];

{{minval, minrules}, {{evalist}, {steplist}}} =
  Reap[
    NMinimize[
      {totdist, a > 0, b > 0},
      {a, b, x, y},
      EvaluationMonitor :> (e++; Sow[{a, b, x, y}, eval]),
      StepMonitor :> (s++; Sow[{a, b, x, y}, step]),
      MaxIterations -> 150,
      Method -> {
         "NelderMead",
         "InitialPoints" -> {Flatten@{StandardDeviation[data], Mean[data]}}
      }
    ],
    {eval, step}
  ];

minrules

Graphics[{
  Red, Circle[{x, y}, {a, b}] /. minrules,
  Black, PointSize[0.02], Point@data, AspectRatio -> 1
}]

running counter

(* Out: {a -> 53.3652, b -> 35.1669, x -> 66.0047, y -> 64.0812} *)

fitted ellipse


One can take a closer look at the intermediate steps in the optimization process that were saved in steplist. Alternatively, evalist can be substituted for steplist if one is interested in each evaluation of the function to be minimized.

ListAnimate[
 Graphics[{
     Red, Circle[{#3, #4}, {#1, #2}],
     Black, PointSize[0.02], Point[data]
    },
    PlotRange -> 
        Transpose @ Through[{Subtract, Plus}[Mean[data], 3 StandardDeviation[data]]],
    Axes -> True
 ] & @@@ steplist,
 AnimationRepetitions -> 1
]

animation of steps


Finally, here is some code to generate a set of "noisier" points from an axis-aligned ellipse:

With[
  {xc = 2, yc = 5, a = 30, b = 22},
  data = (# + RandomReal[{-1, 1}] &) /@
    FindInstance[(x - xc)^2/a^2 + (y - yc)^2/b^2 == 1, {x, y}, Reals, 15][[All, All, 2]]
];

... and here is an example of a best-fit ellipse obtained from such points. This required more iterations (MaxIterations was set to 600), but still only 14.3 s overall:

noisier ellipse

MarcoB
  • 67,153
  • 18
  • 91
  • 189
3

this does not directly address the issue, but this works using RegionDistance in a purely numeric way:

data = {{13, 60}, {36, 35}, {70, 29}, {106, 41}, {118, 72}, {94, 
   94}, {49, 99}, {24, 85}}; 
d[{x_?NumericQ, y_?NumericQ}, a_?NumericQ, b_?NumericQ, data_] := 
 RegionDistance[Circle[{x, y}, {Abs[a], Abs[b]}], #] & /@ data // 
   Total // N

f=NMinimize[d[{x, y}, a, b, data], {x, y, a, b}]

{2.35102, {x -> 65.9529, y -> 64.0925, a -> 53.3151, b -> 35.1669}}

  Show[{RegionPlot[ Circle[{x, y}, {Abs[a], Abs[b]}] /. Last@f ], 
         ListPlot[data]}, AspectRatio -> 1, 
           PlotRange -> {{0, 120}, {0, 120}}]

enter image description here

BTW This is really slow.. see here for some better approaches: Fitting points to tilted, off-center ellipse

george2079
  • 38,913
  • 1
  • 43
  • 110
  • Thanks for your suggestion.Actually I can use FindFit to do it,learning from a friend.But I don't like the FindFit because of its initial value hard to find.And the new function of Geometry Computation is very convenience I thank.your purely numeric way is work.As you say,It's very slow.But all the same RegionDistance[Circle[{x, y}, {Abs[a], Abs[b]}], {13, 60}] will failure to run.Can you give any opinion? – yode Oct 14 '15 at 16:42
  • @yode will the symbolic RegionDistance really fail with an error, or will it just run for a very long time? – MarcoB Oct 14 '15 at 16:55
  • @MarcoB Will collapse my Mthematica in my computer.It is running normally in your?? – yode Oct 14 '15 at 17:28
  • 1
    @yode Do you do understand you are requesting an analytic expression for nearest distance to an arbitrary ellipse? This problem is intractable, or if a solution can be produced it will be a uselessly cumbersome conditional expression involving roots of transcendental expressions. – george2079 Oct 16 '15 at 14:50
  • @george2079 I'm sorry for I don't know how complicated it is actually.I don't know too many problem that I can figure out it with Mathematica.And I think It should be a matter of concern by Wolfram Research.A good tool should help even a common people,which is a Mathematica's pursuit,isn't it?BTW even though kind you and @MarcoB give some solution for this problem so that I can solve it now,I'm still not satisfied with the behavior of RegionDistance in such case. – yode Oct 16 '15 at 15:27
  • 3
    @yode This is not really a shortcoming of RegionDistance: there is no simple and general expression for the distance of a point from an arbitrary ellipse. Take a look at this file (Distance from a Point to an Ellipse, an Ellipsoid, or a Hyperellipsoid). – MarcoB Oct 16 '15 at 16:48
  • @MarcoB Oh,thanks for your PDF to help me learn a lot. – yode Oct 17 '15 at 05:51