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'm writing a code to do rotational averages of image FFT's. I have a working version but it is pretty slow. Wondering if anyone has any thoughts on a way to do this more efficiently.

Data looks like this generally:

rawData FFT

    img1 = Import["http://i.imgur.com/0QZs0AN.png"];
    img = ImageAdjust[ImagePeriodogram[img1], 1];

    imgpts[r_, \[CapitalDelta]_, imagedimension_] := 
      If[r != 0 , Table[{r Cos[2 \[Pi] t] + imagedimension/2, 
    r Sin[2 \[Pi] t] 
    + imagedimension/2} // N, {t, 0, 1, 
    1/(\[CapitalDelta] r)}], {{imagedimension/2, imagedimension/2},
    {imagedimension/2, imagedimension/2}} // N];

    RotationalAverage[FFT_, imgdim_, rmax_, \[CapitalDelta]_, inter_] := 
      Table[{r, Mean[ImageValue[FFT, #, Resampling -> inter] & /@ 
      imgpts[r, \[CapitalDelta], imgdim]]}, {r, 0, rmax}];

    imagevals = RotationalAverage[img, 1024, 140, 16, "Bicubic"];

I decided to make each circle average over bins of the same arc-length. Im not sure if this is the "right" way to do this, but it seemed like a good idea. This gets me the data I want and I usually make a plot at this point just to check.

    ListPlot[imagevals]

plot for checking

Then i stitch the images together to get a picture.

    ImgData = Reverse@Table[If[N[Sqrt[x^2 + y^2]] < 140,         
    Interpolation[imagevals][Sqrt[x^2 + y^2] // N], 0], {x, 0, 128}, 
    {y, 0, 128}] // Quiet;

    ImageAdjust[ImageAssemble[{{Image[Reverse@Transpose@ImgData], 
    Image[ImgData]}, {Image[Transpose@ImgData],Image[Reverse@ImgData]}}], 1]

To get something like this

final product

What I'm ultimately after i guess is some optimization of accurate and fast, which is pretty vague but i guess i don't know how else to say it.

share|improve this question

2 Answers 2

up vote 10 down vote accepted

We can transform the image into polar coordinates, after which averaging across angles is trivial.

polarTransform[img_, rmax_] := 
 With[{size = Max@ImageDimensions[img]}, 
  ImageTransformation[img, 
   Function[{r, t}, {r Cos[t], r Sin[t]}] @@ # &, {rmax + 1, 2 Pi rmax}, 
   DataRange -> {{-size/2, size/2}, {-size/2, size/2}}, 
   PlotRange -> {{-1/2, rmax + 1/2}, {0, 2 Pi}}, Resampling -> "Bicubic"]]
rmax = 140;
polarImg = polarTransform[img, rmax];
imageVals = Mean /@ Transpose@ImageData@polarImg;
ListPlot[imageVals, DataRange -> {0, rmax}, PlotRange -> All]

enter image description here

On my machine, this runs about 20 times faster than the corresponding part of your code, probably because ImageTransformation is a highly optimized internal function. Maybe you'd get similar performance if you Compile'd your code?

Anyway, heck, now we can use ImageTransformation again to create the rotationally averaged image.

w = 128;
ImageTransformation[Image[{imageVals}], {Norm[#], 0} &, {2 w + 1, 2 w + 1}, 
 DataRange -> {{-1/2, rmax + 1/2}, {0, 2 Pi}},
 PlotRange -> {{-w - 1/2, w + 1/2}, {-w - 1/2, w + 1/2}}, Resampling -> "Bicubic"]

enter image description here

All the $+\frac12$'s and $-\frac12$'s in the DataRange and PlotRange are so that the pixels themselves lie at integer coordinates; see this question and compare the documentation for ImageValue and PixelValue. So, for example, by setting the range of $r$ to be $\{-1/2, r_\max + 1/2\}$ in the first ImageTransformation, we get samples at $r = 0, 1, ..., r_\max$ instead of $\frac12, 1\!\frac12, ..., r_\max - \frac12$.

share|improve this answer
    
@RahulNarian I've been looking at this for a while and comparing it to my code. One thing that is bothering me is this plot of the difference between the methods as a function of radius. There seems to be some sort of systematic effect at small radii that I don't understand. –  bshev Apr 15 at 18:04
    
@bshev: Good catch! My previous code was off by half a pixel. Please see the updated answer. There is much less of a difference with your results now. –  Rahul Narain Apr 22 at 4:08

How about using nesting ImageRotate:

t = NestList[ImageRotate[#, 5 Degree, Full] &, img, 71];
ImageData /@ t;
out = Total[%] // Image // ImageAdjust

Mathematica graphics

Extracting a line looks similar to yours

idata = ImageData[out];
ListPlot[idata[[All, 512]], PlotRange -> {{512, 652}, {0.4, 1}}, 
 Frame -> True] 

Mathematica graphics

share|improve this answer
    
Not sure this is faster than the method in the question. –  Rahul Narain Apr 13 at 3:32
    
@RahulNarain I agree I haven't tested it, but I am also performing the operation on the entire image, and from what I gather in the OP's code, only the data ~140 pixels from the center are being analyzed. –  bobthechemist Apr 13 at 11:49

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.