Basically, you want to fit a shape to a set of points with outliers. One common algorithm to do this is RANSAC (random sample consensus). The basic outline of this algorithm is:
- Select N points at random (where N is the minimum number of points required for fitting the shape, i.e. 2 for a line, 3 for a circle and so on)
- Fit the shape to these points
- Repeat 1&2, pick the "best" shape (where "best" usually means closest to a randomly sampled test point - I'll use a simpler criterion below.)
- Select the points close to that shape and fit an improved estimate to those points
Let's try this. (Or a slight modification of it, anyway.)
First, I perform standard edge extraction:
img = Import["https://i.sstatic.net/H63BK.jpg"];
edges = DeleteBorderComponents[EdgeDetect[img,5]];
edgePos = N[PixelValuePositions[edges, 1]];
Next, I define the shape I want to fit - in this case, a conic:
conic[{x_, y_}] := a[1] + a[2]*x + a[3]*y + a[4]*x*y + (1 + a[5])*x^2 + (1 - a[5])*y^2
conicNormal[{x_, y_}] = D[conic[{x, y}], {{x, y}}]
isValidEllipse[solution_] := Abs[a[4] /. solution[[2]]] < 0.1
isValidEllipse is a simple ad-hoc criterion to filter out completely deformed conics.
Normally, we'd need 5 points to fit a conic. The image contains about 50% outliers, so the chance to randomly select 5 points on the ellipse we're looking for is about $0.5^5$, which is not good. We'd need a lot of samples to be reasonably sure that at least one 5-tuple of points contains only non-outlier points.
But fortunately, the location of the edge is not the only information we have: We also know the local orientation at each pixel, which should be perpendicular to the conic's normal vector.
orientation = GradientOrientationFilter[img, 3];
Now 3 points and their orientations give us 6 equations, so we have an overdetermined equation system:
fitEllipse[samplePoints_] := Quiet[FindMinimum[
Total[
Flatten[
{
(conic /@ samplePoints)^2,
With[{\[Alpha] = PixelValue[orientation, #]},
(conicNormal[#].{Cos[\[Alpha]], Sin[\[Alpha]]})^2] & /@ samplePoints
}]], Array[a, 5]]]
This function returns not only the best-fit conic's coefficients, but also the residual error, which is a cheap way to compare between different fitted conics. The assumption is: If randomly sampled 3 points are all parts of the ellipse, the residual error will be low. If the points don't belong to one conic, the residual error will be much higher.
potentialSolutions =
SortBy[Select[Table[fitEllipse[RandomSample[edgePos, 3]], {n, 100}],
isValidEllipse], First];
result = potentialSolutions[[1]];
(There is room for improvement here: the ellipse you're looking for might not be contained in these 100 samples, or it might not be the one with the lowest residual error. A smarter algorithm would take e.g. 50 samples, take the best one or two of these and count the number of nearby points. If enough points are close to it, keep that ellipse, otherwise, take 50 new samples. But that's just a straightforward programming exercise left to the reader ;-) .)
Next step: Find points nearest to the conic. I've tried to calculate the geometric distance exactly (using FindMinimumValue) but that's very slow. A simple, fast approximation is to find N points on the ellipse and simply use Nearest to estimate the distance:
ellipsePts =
Select[Flatten[
Table[{i, y} /.
Solve[conic[{x, y}] == 0 /. result[[2]] /. x -> i, y], {i, 0, w,
10}], 1], Total[Abs[Im[#]]] == 0 &];
nf = Nearest[ellipsePts];
Now, we can calculate the distance of every edge pixel to this conic and pick the ones that are closer than some threshold:
distances = Norm[nf[#][[1]] - #] & /@ edgePos;
closestPoints = Pick[edgePos, # < 10 & /@ distances];
...and improve the ellipse estimate based on these points:
result = fitEllipse[closestPoints];
Repeat the last steps until convergence. (Possible improvements: You could try to reduce the distance threshold in each iteration. Or you could add a "weight" to the curve fitting and give points closer to the current estimate a higher weight when calculating the next estimate.)
The result (after 5 iterations) looks like this:
Show[
img,
ContourPlot[
Evaluate[conic[{x, y}] == 0 /. result[[2]]], {x, 0, w}, {y, 0, h},
ContourStyle -> {Red, Thick}],
ListPlot[closestPoints]]
Disclaimer: The result will not always look like this, as the algorithm is randomized and the algorithm can get stuck in local minima.
Image[ ImageData@EdgeDetect[pic] - ImageData@EdgeDetect[pic, "StraightEdges" -> 0.2]]and then applyMorphologicalComponentsand then on each component (one will be the ellipse, you may be able to delete the other somehow) useComponentMeasurements. Look in the documentation for the properties that it can give you, that you may want, e.g. best fit ellipse semiaxes and so on. $\endgroup$