I tried to detect the rim of the following cup as an ellipse.

some cup

I've tried the solutions given in How to find circular objects in an image? to detect the ellipse, but the detection result was not quite what I wanted. The reason might be that ChanVeseBinarize function can't separate the rim from the body. Therefore SelectComponents will view the whole cup as one component and consequently detect nothing.

Then I tried to use edge detection. After manipulating the parameters, I can detect the edge of the rim like this:

edge detection attempt

But I can't find an algorithm that detect the line of an ellipse. How might I do that?

  • $\begingroup$ What I can think of right now: try to isolate the ellipse e.g. through Image[ ImageData@EdgeDetect[pic] - ImageData@EdgeDetect[pic, "StraightEdges" -> 0.2]] and then apply MorphologicalComponents and then on each component (one will be the ellipse, you may be able to delete the other somehow) use ComponentMeasurements. 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$ – C. E. May 21 '13 at 13:34
  • $\begingroup$ @RodLm It would be easier If he did that himself when editting his question. Now it's been 5 times editted. Wouldn't like to make it a community-wiki question. $\endgroup$ – Artes May 21 '13 at 13:46
  • $\begingroup$ Well you will need to use a fuzzy ellipse. The basic idea is to generate a number of ellipses until one of them is such that most of its points are within a certain distance of a white point of your edge detection output. Now the question is reduced to how to effectively generate ellipses. $\endgroup$ – SMeznaric May 21 '13 at 14:05
  • 1
    $\begingroup$ Substantial progress towards a solution can be found by reviewing these related posts. A Hough transform looks attractive. $\endgroup$ – whuber May 21 '13 at 15:03
  • 3
    $\begingroup$ I'm not sure about the hough transform: an ellipse has 5 degrees of freedom, that's a huge parameter space to update for each point. $\endgroup$ – Niki Estner May 21 '13 at 17:01

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:

  1. 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)
  2. Fit the shape to these points
  3. 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.)
  4. 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["http://i.stack.imgur.com/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[
      (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 = 
    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:

  Evaluate[conic[{x, y}] == 0 /. result[[2]]], {x, 0, w}, {y, 0, h}, 
  ContourStyle -> {Red, Thick}],

enter image description here

Disclaimer: The result will not always look like this, as the algorithm is randomized and the algorithm can get stuck in local minima.

  • 6
    $\begingroup$ +1 That's a real tour de force, beautifully described. Impressive. $\endgroup$ – whuber May 21 '13 at 22:27
  • $\begingroup$ Code is great, and descriptions we can learn from even greater! $\endgroup$ – cormullion May 22 '13 at 7:25
  • $\begingroup$ A great detailed tutorial. $\endgroup$ – Ali May 27 '13 at 17:42

Too late and not so general as @Nikie but I will still like to answer in order to highlight a method involving generalized Eigen values! Lets take the image.

(* Detecting ellipsoidal edges and removing the straight line component *)
i = ImageSubtract @@ {EdgeDetect[img, 10],EdgeDetect[img, 10, "StraightEdges" -> .2]};
{wid, hig} = ImageDimensions[i];
(* Taking control of the two straight line segments *)
lines = ImageLines[EdgeDetect[img, 10, 0.025], .15,"Segmented" -> True];
(* Find the mid point of longest line and its hight *)
clip = Ceiling@Abs@Mean[Last@Flatten@Differences@
          First@Sort[#,(EuclideanDistance @@@ #1 >EuclideanDistance @@@ #2)&]&/@lines];
(*Trim the bottom content using above clip and extract the points*)
padim = ImagePad[ImageTake[i, clip], {{0, 0}, {hig - clip, 0}}, Black];
list = Position[Transpose@ImageData[padim, "Bit", DataReversed -> True], 1];

Lets see the above steps of image processing results clockwise starting from top left.

image = Show[img, Graphics[{Thick, Blue, Line /@ lines}],ImageSize -> 300];
GraphicsGrid[{{img, ColorNegate@i}, {image, ColorNegate@padim}},ImageSize -> 800,
Frame -> All, Spacings -> 0]

enter image description here

At this point we can use EigenSystem efficiently as very simplistically discussed in here!

{X, Y} = Transpose[list];
(* design matrix *)
d = Transpose@{X^2, X Y, Y^2, X , Y, ConstantArray[1, Length@X]};
s = Transpose@d . d;
(* constraints *)
c = ConstantArray[0, {6, 6}];
c[[1, 3]] = -2;
c[[2, 2]] = 1;
c[[3, 1]] = -2;
{val, vec} = Eigensystem[{N[s, 200], c}];
(* parameter for a general conic*)
{a, b, c, d, e, f} = vec // Last;

Now we are ready to check out the fitting.

(* Visualize fit *)
ptEig = Show[img, 
a x^2 + b x y + c y^2 + d x + e y + f == 0, {x, 0, 1385}, {y, 0,1023},
ContourStyle -> {Thickness[0.004], Blue}, 
PlotPoints -> 200, MaxRecursion -> 4, PlotRange -> All, 
Frame -> None, Background -> None], 
Graphics[{PointSize[0.001], Red, Point[#] & /@ list}], 
ImageSize -> 400]

enter image description here


This didn't work as well as I'd hoped, but it's quite interesting to see what went wrong.

i = Import["http://dailycoffeenews.com/wp-content/uploads/2011/09/disposable-coffee-cup.jpeg"];

j = ColorConvert[i, "GrayScale"];
k = Thinning[
             EdgeDetect[j], 1], 1]]]

which gives us a trace of the cup:


(Aargh, I've lost the rim. Pressing on anyway:)

pts = PixelValuePositions[k, "Max"];

left = First@Sort[pts, #1[[1]] < #2[[1]] &];
right = Last@Sort[pts, #1[[1]] < #2[[1]] &];
top = First@Sort[pts, #1[[2]] > #2[[2]] &];
centre = Mean[{left, right}];
rim = Disk[centre, {First[centre - left], Last[top - centre]}];

 Graphics[{Opacity[0.5], Brown, rim}]]


It's close, but not close enough. And now I can see where I went wrong - I've assumed the cup is aligned with the camera and geometrically accurate, but I don't think it is. (And what exactly is the rim anyway?) I wonder whether this method would work if the cup was aligned.


Here's another approach that uses some of the built-in morphological (image processing) functions. Begin by reading in the figure, detecting the edges, and then making them thick:

fat = Dilation[EdgeDetect[i], 20]

enter image description here

Get just the outline of the cup:

thin = Dilation[Thinning[fat], 1]

enter image description here

(it's actually connected when looking at it full size). Now use MorphologicalComponents to get the different pieces of the image

components = MorphologicalComponents[ColorNegate[thin]]; 
components // Colorize

enter image description here

And finally ask for some relevant properties of the two regions. The centroids of the two regions are:

ComponentMeasurements[components, "Centroid", BorderComponents -> False]

{2 -> {701.478, 783.066}, 3 -> {696.272, 454.265}}

So region number 2 is the purplish cup rim and region 3 is the greenish side of the cup. A variety of other measurements can be easily extracted:

ComponentMeasurements[components, {"Length", "Width", "SemiAxes", 
      "Orientation", "Elongation", "Eccentricity"}, BorderComponents -> False]

{2 -> {619.841, 179.126, {309.921, 89.5632}, -3.14045, 0.711013, 0.957333}, 
 3 -> {619.414, 556.459, {309.707, 278.229}, -1.39043, 0.101637, 0.439253}}

The best ellipsoidal fit for the rim has the given length, width, semiaxes, orientation, elongation, and eccentricity. The other is the best ellipsoidal fit for the cup sides (which is not a good fit at all).


My approach to to similar, but somewhat fuzzier, problem, finding an ellipse as a border of galaxy galaxy

binarized = MorphologicalBinarize[galaxy, {0.165, 0.93}]

Finding the parameters for binarize can be tricky, these values do well for this specific case.


components = MorphologicalComponents[binarized];
params = ComponentMeasurements[
    components, {"Centroid", "SemiAxes", "Orientation"}] // Values;
 Graphics[{Red, {
      Rotate[Circle[#[[1]], #[[2]]], #[[3]]]
   } & /@ params}]

The result:

enter image description here


Your Answer

