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.

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:

Mathematica graphics

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]

Mathematica graphics

Now,

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

Mathematica graphics

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

share|improve this question
    
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 Narain 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
show 4 more comments

2 Answers 2

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]

Mathematica graphics

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

Mathematica graphics


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

Mathematica graphics

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]

Mathematica graphics

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.

share|improve this answer
    
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
add comment

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]

enter image description here

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

enter image description here

enter image description here


* 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.

share|improve this answer
2  
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 Narain Dec 8 '12 at 11:16
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.