# Find radii of concentric circles in image

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


or, using :

i = Import["http://i.imgur.com/oTTM9MG.jpg"];
disk = ColorNegate[Binarize[i, {0.3, 1}]];
rings = ComponentMeasurements[
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!

• 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

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}
*)


• 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}}]


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,
True][[All, 1]];
peaks = SortBy[Transpose[{peakX, radiusStrength[[peakX]]}], Last];
ListLinePlot[radiusStrength, Epilog -> {Red, Point[peaks[[-4 ;;]]]}]


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]]}]]


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.

• +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? – Niki Estner 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. – Niki Estner 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
• @Reid This online event should be worth catching. – cormullion Jun 14 '13 at 16:46