# Counting radial ridges on an image

After happily using the v9 image assistant to crop elliptically an image, and then the drawing tools to put a white disk in the middle, I turned this image:

into the following one:

that can be imported with

im=Import["http://i.stack.imgur.com/NNzNM.png"]


The overall objective is to programatically count the number of those radial lines. Due to the lighting, there are parts of the image in which those lines are darker than its surroundings and others where it is lighter.

So far I haven't found a good way worth posting, so any pointer to a good solution would be appreciated. I have the feeling the image processing people will see better ways of dealing with this and I will be grateful to learn something. Thanks a lot

EDIT

A first approach with @RahulNarain's suggestion would be

int = ListInterpolation[
ImageData[ColorConvert[im, "GrayLevel"]], {{-1, 1}, {-1, 1}}];
polInt = Function[t, int[0.95 Cos[2 \[Pi] t], 0.95 Sin[2 \[Pi] t]]];

Plot[polInt[t], {t, 0, 1}, AspectRatio -> 0.2, ImageSize -> Large,
PlotRange -> Full]


Now,

ListLinePlot[Abs@Fourier[polInt@Range[0, 1 - 0.001, 0.001]],
PlotRange -> {{0, 500}, {0, 50}}]


A better zoom shows the maximum at 98

However, manual counting (could be wrong) gave me 96, and nikie's approach is suggesting 97. 1 or 2 off count could be due to the light changes making the real ridge be the local minimum at some places and local maximum at others

-
I'm by no means proficient in image manipulation, but Sharpen may be what you're looking for - nest it 10 times. –  VF1 Dec 8 '12 at 6:10
Maybe: build an interpolating function from the image, sample it over a circular path something like $(0.95\cos\theta, 0.95\sin\theta)$, and then analyze the variation of intensity as a function of $\theta$. That reduces it to a one-dimensional signal processing problem. –  Rahul Dec 8 '12 at 7:40
@RahulNarain nice idea, I'll give it a shot now :) –  Rojo Dec 8 '12 at 7:41
The big version image is rather tiny... –  Yves Klett Dec 8 '12 at 7:58
@YvesKlett haha, I thought it wasn't appropriate to upload the big version, but I can do it if you want –  Rojo Dec 8 '12 at 7:59

A simple first try would be to rotate the image, then measure the distance to the original image:

img = ColorConvert[Import["http://i.stack.imgur.com/NNzNM.png"],
"Grayscale"];

Monitor[t =
Table[{i,
ImageDistance[img, ImageRotate[img, 360°/i, Full, Background -> White]]}, {i, 10,
200}], i];

ListLinePlot[t]


Obviously, if the angle is exactly 360° / [number of radial lines], the match should be lowest, so the estimated count would be:

count = Extract[t, Position[t[[All, 2]], Min[t[[All, 2]]]]][[1, 1]]


which is 97

If I overlay 97 radial lines over your image, it seems as if the count wasn't too far of. I can't tell if it's exact, though:

center = 0.5 ImageDimensions[img];
Show[img,
Graphics[
{Red, Table[
Line[{center + 80 {Cos[\[Phi]], Sin[\[Phi]]},
center + 250 {Cos[\[Phi]], Sin[\[Phi]]}}], {\[Phi], 0, 2 \[Pi],
2 \[Pi]/count}]}]]


EDIT: I've been playing with this some more, especially with the FFT idea. First, I polar-transform the image:

polar = ImageTransformation[
img, #[[2]]*{Cos[#[[1]]], Sin[#[[1]]]} &, {500, 20},
PlotRange -> {{0, 2 \[Pi]}, {0.9, 1.0}},
DataRange -> {{-1, 1}, {-1, 1}}]


Then I've applied a windowed Fourier transform to the mean of that signal:

mean = Mean[ImageData[polar]];
window = Array[HannWindow, Length[mean], {-1.5, 1.5}];
stft = Table[(Abs[Fourier[mean*RotateLeft[window, i]]][[
80 ;; 120]]^2) // #/Max[#] &, {i, 0, 500}];
ArrayPlot[stft\[Transpose], ColorFunction -> GrayLevel,
DataRange -> {{1, 500}, {80, 120}}, FrameTicks -> True]


The windowed Fourier transform looks as if the frequency isn't constant over the whole area. Which would make sense, if the center isn't perfect or if there's an affine/perspective transformation. Sadly, I'm not sure what to do with this, but I thought I'd post it, in case it gives somebody else an idea.

-
Very smart approach! Thanks a lot! What I'm trying using Rahul's suggestion so far is giving me 98 or 99, but counting them manually gives me 96. Haven't triple checked manually, you can imagine why –  Rojo Dec 8 '12 at 8:48
If you want 96, you could use DistanceFunction -> "MutualInformationVariation". That's not very scientific, though. ;-) –  nikie Dec 8 '12 at 8:49
After your polar transform, do a column sum to make it 1d, and then MaxDetect would seem suited for the problem? –  JxB Dec 8 '12 at 14:23
+1 very original –  Vitaliy Kaurov Dec 8 '12 at 15:22
Nice edit. Perhaps you could have used SpectrogramArray for the windowed fourier transform –  Rojo Dec 8 '12 at 18:13

I think the big bright spots in the lower half of the image are confusing the Fourier transform, because they lie right in between where two bright lines should be, and so they are exactly out of phase with the rest of the signal. How about I just throw them away?*

polInt = Function[t,
If[(t > 0.11 && t < 0.16) || (t > 0.89), 0.2,
int[0.95 Cos[2 \[Pi] t], 0.95 Sin[2 \[Pi] t]]]];
Plot[polInt[t], {t, 0, 1}, AspectRatio -> 0.2, ImageSize -> Large,
PlotRange -> Full]


Then the Fourier transform has a maximum at 96, as you want.

* This procedure is entirely unscientific. Chopping data at hard boundaries is not recommended by signal processing experts. Do not try this at home. Void where prohibited.

-
Actually... Does the output of Fourier have frequency 0 (the DC component) at position 1? Then the peak here is at 95 cycles, which might be a bad thing. –  Rahul Dec 8 '12 at 11:16