# Is there a fast method of generating a Standard Deviation image from a frame stack?

Imagine I have a stack of image frames, $(f_1, ..., f_N) \in F$, where each $f_i$ has the same set of $(x,y)$ dimensions. Using these frames, I'd like to create an image where each pixel in the image consists of the standard deviation of the values at the corresponding pixel position in each of the $f_i$.

Unfortunately, the following procedure is painfully slow (when applied to a $256 \times 256$ pixel image):

SDImage = ImageData@FrameStack[[1]];
For[a = 1, a <= ImageDimensions[FrameStack[[1]]][[1]], a++,
For[b = 1, b <= ImageDimensions[FrameStack[[1]]][[2]], b++,
SDImage[[a, b]] = StandardDeviation[Table[ImageData[FrameStack[[k]]][[a, b]], {k, 1, Length[FrameStack]}]];
Print[SDImage[[a, b]]];
];
];

SDImage = Image[SDImage]


Is there a faster method of doing this, or perhaps a built in tool (like in ImageJ)?

-
By the way: The code you posted in your question might be painful slow because the procedural nature of the code. Trying to avoid For loops and using Mathematica build in functions is always a good way in making the code more efficient. – g3kk0 Jul 27 '13 at 8:15
Related: Averagind 2D datasets. Or different approach for work with multiple images. Offtop: try to avoid you variable names starting with capital letters. – Kuba Jul 27 '13 at 8:16

Assuming you have loaded a stack of images and stored in a variable called imageStack you can do the following to compute the standard deviation image over the entire stack:

imageStack = Import["ExampleData/CTengine.tiff"];
xyStd = Image@Thread[StandardDeviation[ImageData /@ imageStack], {1}];


This gives you the standard deviation image in z direction of the stack. For the upper example xyStd gives the following image:

Edit

As Kuba pointed out, a better solution would be:

StandardDeviation[ImageData /@ imageStack]

-
StandardDeviation[ImageData /@ imageStack] is enough. – Kuba Jul 27 '13 at 8:02
@Kuba: right, thanks ;) – g3kk0 Jul 27 '13 at 8:11
No problem :) I've made that mistake here – Kuba Jul 27 '13 at 8:15

This is about 100 times faster than g3kk0's answer on my machine---the results are identical:

imageStack = Image[ColorCombine@Import["ExampleData/CTengine.tiff", "ImageList"], "Real"];
ImageApply[StandardDeviation, imageStack]

-
I can confirm the speedup. Thank you! – Sparse Pine Jul 31 '13 at 5:23
Nice usage of ColorCombine and ImageApply. I can confirm a speedup of about 25 fold compared to my version. – g3kk0 Jul 31 '13 at 9:51

I wondered whether the built-in StandardDeviationFilter would work here (on Alexander's imageStack).

 StandardDeviationFilter[
Image3D[
imageStack,
ColorFunction -> "LowRange"],
{Length[imageStack], 1, 1}]


It kind of works - but it's probably not the solution you're looking for.

-
That's a good suggestion - but I previously tried this, and in my particular case, the results are radically different. – Sparse Pine Jul 27 '13 at 12:30
@SparsePine Oh. – cormullion Jul 27 '13 at 12:32
Perhaps this is better StandardDeviationFilter[Image3D[imageStack], {Ceiling@Length[imageStack]/2, 0, 0}] -- the last arguments are radii. But it's several times slower than Alexander Schmitz's. – Michael E2 Jul 27 '13 at 12:42
@cormullion Sorry, I think I came off as a bit rude there. That really wasn't my intention. I just meant to say that I tried the filter, and it just didn't work very well in my particular case. Thank you for your answer! – Sparse Pine Jul 27 '13 at 12:53
@SparsePine No, I didn't think that... I was just surprised that you'd tried it. It's always worth mentioning what you've tried (and providing a small sample of your data). Image3d isn't very quick or dependable on my machine anyway... – cormullion Jul 27 '13 at 13:03