Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I need your help to detect the radius of the blue circle in images like this:

Mathematica graphics.

Reading these questions: 1, 2, 3, I figured that this code taken from @halirutan answer to question n°3 could help:

imga=Import[http://i.stack.imgur.com/moV5x.png];
circle = SelectComponents[
MorphologicalComponents[
LaplacianGaussianFilter[ColorNegate@imga, 2], 0.003`], 
"Count", # > 400 &];
Colorize[circle]

And I have this output:

Mathematica graphics

And now I'm stuck, if the circle were closed it would be simple, how can I find the best circle completing it? Obviously if you have a completely different solution it is fine (for example I thought it could be possible to take only the left semicircle, but I don't know how to do it).

Please note that I have to process about 100 images like that so if possible I would need the program to run in a reasonable amount of time and with no manual adjustments.

As a bonus question: since I have to determine also a possible error on the radius, would it be possible to find the outer and the inner circles? and the eccentricity of the closed figure (if it is possible to close it)?

P.S.: since I just started working with mathematica for image processing I would also like to know a few good references or books about it. I have read a few mathematica books but none of them treated this argument in detail.

share|improve this question
    
possible duplicate of Derive a smooth circle with cusp from an image –  Murta Feb 20 at 11:47
1  
@Murta: It's not a duplicate (I think): The answers in that question all use the center of mass of the marked pixels (the centroid). When the circle isn't closed or the "thickness" varies, the center of mass is not the center of the best-fit circle. –  nikie Feb 20 at 11:56
    
@nikie, I agree. Close retracted. –  Murta Feb 20 at 12:15
    
That was exactly my point, thank you! –  John Feb 20 at 14:28
    
If you can grab a few points on the partial circle (at least 3 but preferably several more) then you can set up a linear least squares problem to find the center and radius (yes, I mean linear). –  Daniel Lichtblau Feb 20 at 19:39
add comment

2 Answers

up vote 8 down vote accepted

You're close: You can use ComponentMeasurements to calculate measurements for each component. (The same measurements that you can use in SelectComponents, too). If you can get away with the center of mass of the marked pixels, it'really quite simple:

imga = Import["http://i.stack.imgur.com/moV5x.png"];
binary = MorphologicalComponents[
  LaplacianGaussianFilter[ColorNegate@imga, 2], 0.003`]; 
components = ComponentMeasurements[
  binary, {"Count", "AreaRadiusCoverage", "Centroid", 
   "MeanCentroidDistance", "MaxCentroidDistance", 
   "MinCentroidDistance"}, #1 > 400 && #2 < 0.01 &];

The "AreaRadiusCoverage" measurement helps to select components where most pixels are near the border. This returns a list of rules, one for each component that matches the filter function:

{193 -> {4516, 0., {222.283, 211.7}, 181.482, 200.342, 167.181}}

Where 193 is the component index, 4516 is the pixel count, 0. is AreaRadiusCoverage and so on.

The found circle looks like this:

Show[HighlightImage[imga, Image[1 - Abs[binary - components[[1, 1]]]]],
 Graphics[
  {
   Yellow, Dashed,
   components /. 
     {(_ -> {count_, areaRadiusCoverage_, center_, radii__}) 
     :> (Circle[center, #] & /@ {radii})}}]]

enter image description here


If the center of mass is not good enough, you can use the "Mask" component measurement to get an image of the marked pixels in each component:

imga = Import["http://i.stack.imgur.com/moV5x.png"];
binary = MorphologicalComponents[
  LaplacianGaussianFilter[ColorNegate@imga, 2], 0.003`]; components = 
 ComponentMeasurements[
  binary, {"Count", "AreaRadiusCoverage", 
   "Mask"}, #1 > 400 && #2 < 0.01 &];

The "mask" value returned looks like this:

Image[components[[1, 2, 3]]]

enter image description here

Getting the coordinates of the white pixels in this image is easy:

whitePixels = PixelValuePositions[Image[components[[1, 2, 3]]], 1];

Now, what we're looking for is a circle with center (cx,cy) and radius r that minimized:

err = Total[
  Map[
      Function[whitePixel, (Norm[{cx, cy} - whitePixel] - r)^2], 
      whitePixels]];

i.e. the squared distance of the white pixels from the circle.

We can simply pass this error function to FindMinimum:

optimalCircle = FindMinimum[err, {cx, cy, r}]

To get a better estimate of the circle:

Show[HighlightImage[imga, Image[components[[1, 2, 3]]]],
 Graphics[
  {
   Yellow, Dashed,
   Circle[{cx, cy}, r] /. optimalCircle[[2]]}]]

enter image description here

(Finding the min/max radius should be straightfoward from here - you have the center and the white pixel coordinates.)

share|improve this answer
    
Thank you for you answer, works great! I managed to write the code of the first part, but as you said I didn't want to use the centroid because it is imprecise in this situation. Can you explain a little more why you set the "AreaRadiusCoverage" to be <0.01 and what exactly it does? –  John Feb 20 at 14:31
    
@John: I think AreaRadiusCoverage (for each connected component) takes a disc centered at the centroid, with the same area as the component, and measures the overlap between that disc and the component. For a ring, where all pixels are far from the centroid, this should always be small. Other suitable measurements might be the ConvexCoverage, or maybe something like (MaxCentroidDistance-MinCentroidDistance)/MeanCentroidDistance –  nikie Feb 20 at 15:15
add comment

rotate the image by some degree and calculate the correlation of both images. it should be small for uncentered circles and large when many points overlap. repeat for as many degrees as you need to be sure.

share|improve this answer
1  
You can add a simple example to illustrate the answer. –  Sektor Feb 20 at 21:08
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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