# Rotational Averaging of Image Periodograms (FFTS)

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:

    img1 = Import["http://i.imgur.com/0QZs0AN.png"];

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]


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;

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


To get something like this

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.

-

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]


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


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

-
@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


Extracting a line looks similar to yours

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


-
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