8
$\begingroup$

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?

$\endgroup$
  • 3
    $\begingroup$ Please wait for confirmation from others who observe this behaviour before adding the "bugs" tag $\endgroup$ – dr.blochwave Oct 14 '15 at 7:48
  • $\begingroup$ @blochwave Oh,I get it.thanks for your prompt. $\endgroup$ – yode Oct 14 '15 at 7:57
  • 2
    $\begingroup$ 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 ). $\endgroup$ – george2079 Oct 14 '15 at 14:16
  • $\begingroup$ @george2079 I'm very sorry to make this typo,and I have corrected it just. $\endgroup$ – yode Oct 14 '15 at 16:26
8
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Very nice answer!! $\endgroup$ – yode Oct 15 '15 at 12:42
3
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – yode Oct 14 '15 at 16:42
  • $\begingroup$ @yode will the symbolic RegionDistance really fail with an error, or will it just run for a very long time? $\endgroup$ – MarcoB Oct 14 '15 at 16:55
  • $\begingroup$ @MarcoB Will collapse my Mthematica in my computer.It is running normally in your?? $\endgroup$ – yode Oct 14 '15 at 17:28
  • 1
    $\begingroup$ @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. $\endgroup$ – george2079 Oct 16 '15 at 14:50
  • 2
    $\begingroup$ @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). $\endgroup$ – MarcoB Oct 16 '15 at 16:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.