# Behaviour of RegionDistance with Ellipse?

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: 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:= ellipse // RegionQ
Out= True


So it is a bug for RegionDistance?

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

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 -> {
"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
}] (* Out: {a -> 53.3652, b -> 35.1669, x -> 66.0047, y -> 64.0812} *) 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
] 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: – yode
Oct 15, 2015 at 12:42

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}}] BTW This is really slow.. see here for some better approaches: Fitting points to tilted, off-center ellipse

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