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 have an image of concentric circles, I would like to find the radii of the circles (only the innermost few are important).

I've had a go using what I could find in previous posts, but am a bit confused about which method I should be using - if I need MorphologicalComponents, or whether to use SelectComponents "count" or "equivalentradius", or colornegate etc. Sometimes the circles are broken (especially when I binarize), so I need to look for incomplete circles too..

So far I have:

i = Import["http://i.imgur.com/oTTM9MG.jpg"];
b = Binarize[i, {0.3, 1}];
m = MorphologicalComponents[b];
c = SelectComponents[m, {"Count", "Holes"}, 
   1000 < #1 < 20000 && #2 > 0 &] // Colorize
ComponentMeasurements[c, {"Centroid", "EquivalentDiskRadius"}]

or, using :

i = Import["http://i.imgur.com/oTTM9MG.jpg"];
disk = ColorNegate[Binarize[i, {0.3, 1}]];
rings = ComponentMeasurements[
   disk, {"Centroid", "EquivalentDiskRadius"}, 
   450 <= #1[[1]] <= 550 && 325 <= #1[[2]] <= 375 && 
     10 <= #2 <= 500 &];
Show[{disk, 
  Graphics[{{Red, Circle[rings[[1, 2]][[1]], rings[[1, 2]][[2]]]}}]}]

Am I making it harder than it is? Could someone bump in the right direction?

both this, How to find circular objects in an image? and this, Finding the centroid of a disk in an image were helpful (but I need to expand it to multiple circles and fit partial circles).

Any help would be appreciated. Thanks!

enter image description here

share|improve this question
    
This reminds me of electron diffraction rings. What is the picture of? –  rcollyer Jun 14 '13 at 15:47
    
Its from a Fabry-Perot etalon en.wikipedia.org/wiki/Fabry%E2%80%93P%C3%A9rot_interferometer (being used to look at the anomalous Zeeman effect in mercury :) ) –  Jeff Jun 14 '13 at 16:01
    
Looked a bit like electron or x-ray diffraction to me too but the central spot from the incident beam looked a bit too weak. –  s0rce Jun 14 '13 at 18:34
    
(I guess I should add, the F.P. etalon is being imaged with a B&W digital camera..) –  Jeff Jun 17 '13 at 15:14

2 Answers 2

up vote 18 down vote accepted

What you could do is apply an edge filter and find the threshold which binarizes your image best:

i = Import["http://i.imgur.com/oTTM9MG.jpg"];
edges = LaplacianGaussianFilter[ColorNegate[i], 2];
Manipulate[Binarize[edges, t], {t, 0, .1}]

After that you could select all objects with a certain radius or you throw out all small objects with a specific "Count". You have to decide then what you prefer as radius. I thought maybe the "MeanCentroidDistance" gives a quite stable measure

circles = 
  SelectComponents[
   MorphologicalComponents[LaplacianGaussianFilter[ColorNegate@i, 2], 0.0056`],
   "Count", # > 300 &];
Colorize[circles]
ComponentMeasurements[circles, "MeanCentroidDistance"]
(*
{26 -> 271.952, 33 -> 262.778, 129 -> 221.202, 157 -> 209.482, 
 329 -> 154.293, 398 -> 136.493}
*)

enter image description here

share|improve this answer
    
Thanks very much! This is along the same lines I was going but much more sophisticated (and your approach works!). I'll try this with all my data (I'm afraid the data may get a little noisy, but I'll get to that later). Thanks! –  Jeff Jun 14 '13 at 16:11
    
Hi, I've played around with this quite abit and like the solution. However, MorphologicalComponents requires the circles to be connected (I think?) Is there a way to use only partial circles? e.g. if I use only the middle section of the previous image, i.imgur.com/JkRlEvW.png –  Jeff Jun 17 '13 at 15:16
    
@Jeff If you have circles which are not complete, than you'll run indeed into problems. Have you tried to work with a circular Hough transform? Mathematica provides a similar approach in ImageLines for linear structures. For arbitrary circles it is computational even more expensive but if you like I can post another answer with some code. –  halirutan Jun 18 '13 at 23:23

You can use the center detection from this answer:

Basically, the center you're looking for is a point in the image for which every gradient points towards or away from the center. That means we can minimize this error term:

squaredError = 1/2 ({cx - x, cy - y}.{-gy, gx})^2;

which leads to a linear equation system:

errDerivative = Expand[D[squaredError, {{cx, cy}}]];
linearSystem = {{D[errDerivative, cx], 
   D[errDerivative, cy]}, -errDerivative /. {cx -> 0, cy -> 0}}

Now we can simply insert the gradients from the image into this equation system:

gradientX = ImageData@GaussianFilter[img, 1, {0, 1}];
gradientY = ImageData@GaussianFilter[img, 1, {1, 0}];
xArr = Array[N[#2] &, Dimensions[gradientX]];
yArr = Array[N[#1] &, Dimensions[gradientX]];
ls = Total[
   linearSystem /. {gx -> gradientX, gy -> gradientY, x -> xArr, 
     y -> yArr}, {-2, -1}];
center = LinearSolve @@ ls;

The center location is in indices (starting with 1) and image processing functions want coordinates (starting with 0, at the bottom left corner), so I have to convert the coordinates:

center[[1]] -= 1;
center[[2]] = Length[gradientX] - center[[2]];

Then I can apply a polar transform

maxRadius = 250;
polar = ImageTransformation[img, 
  center + {Cos[#[[1]]], Sin[#[[1]]]}*#[[2]] &, {360, maxRadius}, 
  DataRange -> Full, 
  PlotRange -> {{0, 360 \[Degree]}, {1, maxRadius}}]

enter image description here

The mean brightness for each row gives a measure of the strength of each radius:

radiusStrength = Mean /@ ImageData[polar, DataReversed -> True];
peakX = Position[
    MapThread[#1 > #2 && #1 > #3 &, {radiusStrength, 
      RotateLeft[radiusStrength], RotateRight[radiusStrength]}], 
    True][[All, 1]];
peaks = SortBy[Transpose[{peakX, radiusStrength[[peakX]]}], Last];
ListLinePlot[radiusStrength, Epilog -> {Red, Point[peaks[[-4 ;;]]]}]

enter image description here

The radii with the highest strengths correspond nicely with the radii you're looking for:

Show[img, 
 Graphics[{Red, Dashed, Circle[center, #] & /@ peaks[[-4 ;;, 1]]}]]

enter image description here

Note that while this is more complex than a 'ComponentMeasurements'-based solution it is much more robust. You don't have to adjust any thresholds or parameters. And it should work fine, even if parts of the circles are occluded or hardly visible.

share|improve this answer
    
+1 Great answer! - is the slight waviness in the polar image showing the fact that the original image isn't circular or is it an optical illusion or inaccuracy? –  cormullion Jun 14 '13 at 14:04
    
@cormullion: I'm not 100% sure. Since the center is a least squares estimate using all gradients in the image, the error could be due to gradients that don't point towards the center. Or maybe my index -> coordinate conversion is off by one? –  nikie Jun 14 '13 at 14:09
    
The coordinate conversion was off by one. Fixed it. I love MMA for image processing, but the distinction between indices and coordinates is ugly. Especially as there is no built-in function to convert them. –  nikie Jun 14 '13 at 14:19
    
Stuff like this is awesome. Can you recommend any resources to get started with MMA image processing? I'd like to play around with it sometime. –  Reid Jun 14 '13 at 14:33
1  
@Reid This online event should be worth catching. –  cormullion Jun 14 '13 at 16:46

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.