Separation of two ellipses defined by coordinates

Question

I have two ellipses that are defined by points (coordinates in 2D space). The center is not at the absolute center (otherwise it was impossible to get it).

What I don't know is which coordinates belong to the first ellipse and which coordinates belong to the second. The goal is to get the coordinates of each ellipse (almost circle) separately. I tried separating them by distance from the center, but since they are not ideal circles, it is not possible to get each one in its entirety.

Is a good way to solve the equation of a circle and increase its diameter and see if there are coordinates in the circle that gain frequency and separate them based on the frequencies?

Thanks for help

Text file with coordinates and image with visualisation

Thank you very much everyone for your help

Answer 1

The figure you attached seems to suggest that you are asking to differentiate between two circles.

But in looking at the actual data provided, there seems to be 4 circles involved, all amost concentric and at different radial position bands from each other.

You can use a curve fit to do a simple separation between the small circle and the larger one. Here is the proposed procedure:

1. Place the origin at the estimated center and adjust the coordinates of all the points to the new origin.

2. Convert the points into polar coordinate function $r(\theta)$ with $$r = \sqrt{x^2+y^2}$$ and $$\theta = {\rm atan2}(x,\;y)$$ _{Some programming environments use atan2(y,x), so check documentation.}

   I understand that the origin does not coincide with the center of the circles/ellipses exactly but with a small error.

3. Do a curve fit to the polar data to the following form $$ r(\theta) \approx r_0 + a \cos \theta + b \sin \theta$$

   This is a 1st order harmonic analysis and the coefficients $(a,b)$ corresponds to a better estimate of the center, and the value $r_0$ is the average circle radius for all the data. We can use this value to differentiate between the large circle and the small circle.

4. Adjust the $(x,y)$ points to the new center and calculate the centered radius and azimuth angle similarly to step 2 above with $$R = \sqrt{(x-a)^2+(y-b)^2}$$ and $$\Theta = {\rm atan2}(x-a,\;y-b)$$

5. Now compare the values of $R$ for each data point, to $r_0$, and if it is larger then the point belongs to the larger circle, otherwise to the smaller circle.

   

   An example of the above procedure with the data values given is shown above.

Link to Google Sheet that does the above

if you look at a histogram of radial distances from the center, you clearly see 3 populations which you can use to classify the data points.

Answer 2

Are your inputs point coordinates? Can you rely on your ellipses not to intersect? If both are true, you could try to approach this as follows:

1. Fit a single ellipse to all your points. The result should lie somewhere between the two ellipses you want to describe.

2. Plug each point into the equation of the ellipsis and observe the error value. That error will have one sign for points in the one side and another sign for points on the other side. Use that information to obtain a set of points that very likely come from the outer ellipse, and another set of points very likely from the inner. I'd probably take the $\tfrac13$ largest (positive) values and the $\tfrac13$ smallest (negative) values, leaving the middle third as “undecided”, but that's just gut feeling.

3. Fit one ellipse to each set of points.

To me an ellipse in this context is actually any conic section. So start with the equation $ax^2+by^2+cxy+dx+ey+f=0$. For an almost circular ellipse you can probably set $a=1$ which will avoid issues stemming from the fact that multiples of this equation describe the same ellipse.

So “fitting an ellipse” would mean setting $a=1$ then determining $b,c,d,e,f$ such that the left hand side is close to zero for all points. Practically speaking, you want to minimize the sums of squares of these left hand sides, making the whole problem a linear least squares optimization problem. You can reduce that to a linear system of five equations.

The “error value” from step 2 would again be the result of plugging a point into the left hand side of the equation. If the point were fully on the ellipse, you'd get a zero. Since your ellipse at that stage is somewhere between the two ellipses you actually want, for the points on the one ellipse you should get positive values, while for those on the other ellipse you would get negative values. (With $a>0$ and the ellipse almost circular you should get positive values outside the ellipse and negative inside, but which side is which doesn't matter that much.)

If your two ellipses intersect, the single ellipse fitted to all the points might do a poor job in separating the points that belong to one ellipse from the points that belong to the other. At that stage some different approach would be needed. You might use the fact that 5 points define an ellipse, trying to find 5 points from the same ellipse and then define a second ellipse based on the points which are furthest away from the ellipse you already found. But getting 5 points from a single ellipse might be a non-trivial task.

Answer 3

Follows a MATHEMATICA script which extracts the more external circle, proceeding until no circle points remains. The points are supposed to be under the label tab. The procedure, after a convenient scaling and reduction to the barycenter, fits the smallest circle which contains the remaining points. At each fit, the points associated to the found circle are eliminated from the remaining data, until no data.

$$ C\to (x-a)^2+(y-b)^2-r^2=0 $$

so given a data set $\{p_k,\ k =1\cdots,n\}$ we cast the following minimization problem

$$ \min_{a,b,r}r\ \ \text{s.t.}\ \ \ \cases{a^2+b^2\le \delta_1\\ 0\le r\le r_{max}\\ \|p_k-(a,b)\|^2\le r^2} $$

Clear[Circlefit]
Circlefit[subtab_, tol_] := Module[{a, b, r, ret = {}, perm = {}, subn, restrs, sol, k, stk, parms, delta = 0.5, rmax = 2},
  subn = Length[subtab];
  restrs = Table[(subtab[[k]] - {a, b}).(subtab[[k]] - {a, b}) <= r^2, {k, 1, subn}];
  sol = Quiet@NMinimize[Join[{r, a^2 + b^2 < delta, 0 < r < rmax}, restrs], {a, b, r}, Method -> "DifferentialEvolution"];
  res[p_] := Abs[Norm[p - {a, b}] - Abs[r]] /. sol[[2]];
  For[k = 1, k <= subn, k++, stk = subtab[[k]]; 
   If[res[stk] < tol, AppendTo[ret, stk], AppendTo[perm, stk]]];
  parms = {a, b, r} /. sol[[2]];
  Return[{ret, perm, parms}]
  ]


data = Import["/home/Downloads/centre_and_coordinates_of_ellipses.txt", "Table"];
data0 = Take[data, {4, Length[data]}];
f = StringReplace[#, "," -> ""] &;
data01 = f /@ data0;
f = StringReplace[#, "(" -> ""] &;
data02 = f /@ data01;
f = StringReplace[#, ")" -> ""] &;
data03 = f /@ data02;
strtoreal[str_] := Internal`StringToDouble@str
tab = Table[{strtoreal[data03[[i, 1]]], strtoreal[data03[[i, 2]]]}, {i, 1, Length[data03]}];

n = Length[tab];
cg = Total[tab]/n;
scale = 1/2000;
tabcg = Table[(tab[[k]] - cg) scale, {k, 1, n}];
subn = 1000;
inds = RandomSample[Range[n], subn];
subtab = Table[tabcg[[inds[[k]]]], {k, 1, subn}];
perm = subtab;
PARMS = {};
tol = 0.03;

While[Length[perm] > 0,
 {ret, perm, parms} = Circlefit[perm, tol];
 tol = 0.02;
 Print[Length[perm], parms];
 AppendTo[PARMS, parms]
]

grtab = Table[Graphics[{Dashed,Circle[{PARMS[[k, 1]], PARMS[[k, 2]]}, PARMS[[k, 3]]]}], {k, 1, Length[PARMS]}];
gr0 = ListPlot[tabcg, AspectRatio -> 1, PlotStyle -> {Blue, Thick}];
Show[gr0, grtab, PlotRange -> All, AspectRatio -> 1, Axes -> False]

NOTE

The minimization process is executed over a randomly selected sample of $1000$ points (subtab) to avoid de burden of too long processing times. The selection of tolerances to extract features is critical. See in the script tol = 0.03 and tol = 0.02. The main procedure Circlefit[data,tol]returns at each step, {ret,perm,parms} which are respectively extracted feature, remaining points, fitting circle parameters found.