ImageJ

Calculate local stdDev, using a sliding box

_‹ Previous Topic Next Topic _›

Classic

List

Threaded

7 messages Options

Steve Wolf

Calculate local stdDev, using a sliding box

Attempting to calculate local stdDev, using a sliding box, on a 4096x4096 image. I’m using the following Python script, but it runs extremely slow, any suggestions on speeding it up?

from ij import IJ
from ij import ImagePlus
from ij.process import FloatProcessor
from array import zeros
from ij.process import ImageStatistics as IS

imp = IJ.getImage()
ip = imp.getProcessor()

options = IS.MEAN | IS.MEDIAN | IS.MIN_MAX | IS.STD_DEV
stats = IS.getStatistics(ip, options, imp.getCalibration())
boxW = 409
boxH = 409
imageW = 4096
imageH = 4096
i = 0

nwidth = 4096
nheight = 4096
pixels = zeros('f', nwidth * nheight)

for x in xrange(0,imageW,1):
for y in xrange (0,imageH,1):
imp.setRoi( y, x, boxW, boxH)
# rimp = imp.getRoi
stats = IS.getStatistics(ip, options, imp.getCalibration())
# print x, ",", y, ",", stats.stdDev
pixels[i] = stats.stdDev
i += 1

#generate image of local stdDev
fp = FloatProcessor(nwidth, nheight, pixels, None)
imp = ImagePlus("Local StdDev", fp)
imp.show()

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

Michael Schmid

Re: Calculate local stdDev, using a sliding box

Hi Steve,

one quick thing: avoid the median in the options, it is not needed and is the slowest operation by far. Also the MIN_MAX is not needed.

For better performance:
Do the calculation by summing up the pixel values and their square for the first for the first ROI (use a double-precision number).
variance = (sumOfSuares - sumOfValues*sumOfValues/n)/n [OR divide by (n-1) instead of n, as you prefer]
stddev = sqrt(variance)
Then, when moving the ROI, just add the values and their squares for the pixels that enter the ROI, and subtract for the pixels that leave the ROI. Then recalculate the stddev from the new sums.

Even faster: Do the sums over the boxWidth for each row (as above), and put them into arrays. Then apply the same strategy in y.

If it needn't be boxes but it can be circles, you can also use the ImageJ-builtin 'variance' filter. It needs about 55 seconds for a 4096x4096 image with radius=204 (diameter 409) on my old computer.

Michael
________________________________________________________________
On Feb 9, 2016, at 17:43, Steve Wolf wrote:

> Attempting to calculate local stdDev, using a sliding box, on a 4096x4096 image. I’m using the following Python script, but it runs extremely slow, any suggestions on speeding it up?
>
>
> from ij import IJ
> from ij import ImagePlus
> from ij.process import FloatProcessor
> from array import zeros
> from ij.process import ImageStatistics as IS
>
> imp = IJ.getImage()
> ip = imp.getProcessor()
>
> options = IS.MEAN | IS.MEDIAN | IS.MIN_MAX | IS.STD_DEV
> stats = IS.getStatistics(ip, options, imp.getCalibration())
> boxW = 409
> boxH = 409
> imageW = 4096
> imageH = 4096
> i = 0
>
> nwidth = 4096
> nheight = 4096
> pixels = zeros('f', nwidth * nheight)
>
> for x in xrange(0,imageW,1):
> for y in xrange (0,imageH,1):
> imp.setRoi( y, x, boxW, boxH)
> # rimp = imp.getRoi
> stats = IS.getStatistics(ip, options, imp.getCalibration())
> # print x, ",", y, ",", stats.stdDev
> pixels[i] = stats.stdDev
> i += 1
>
> #generate image of local stdDev
> fp = FloatProcessor(nwidth, nheight, pixels, None)
> imp = ImagePlus("Local StdDev", fp)
> imp.show()
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

Kenneth Sloan-2

Re: Calculate local stdDev, using a sliding box

In reply to this post by Steve Wolf

There are ways to stably update mean and stdDev - but you would have to roll your own.

Briefly, you keep track of Sum(x) and Sum(x*x). When you move the window, you can subtract out the column (row) you are leaving behind and add in the column(row) you are adding. to be extra clever, move in a serpentine pattern so you always move by 1 row or 1 col.

Be careful though - make sure all your accumulators stay within bounds.

At the level you are using now, I don’t think a speedup is possible, but I don’t know that for sure. You may be making as many as 5 (unlikely, but…) passes over each window (plus a sort, if Median is a true median). Perhaps I’m wrong about that - but I’ll bet there is at least 1, and perhaps 2 (I would probably use 2), complete passes over each window as you go. The package can’t know your pattern of ROIs, and can’t optimize. You can - but you have to write your own low level code to do that.

One more thing - unless I mis-read the code, you have edge effects issues, and the stdDev image is translated wrt the original image. If I’m wrong, I’m sure someone will correct me.

--
Kenneth Sloan
[hidden email]
Vision is the art of seeing what is invisible to others.

> On Feb 9, 2016, at 10:43 , Steve Wolf <[hidden email]> wrote:
>
> Attempting to calculate local stdDev, using a sliding box, on a 4096x4096 image. I’m using the following Python script, but it runs extremely slow, any suggestions on speeding it up?
>
>
> from ij import IJ
> from ij import ImagePlus
> from ij.process import FloatProcessor
> from array import zeros
> from ij.process import ImageStatistics as IS
>
> imp = IJ.getImage()
> ip = imp.getProcessor()
>
> options = IS.MEAN | IS.MEDIAN | IS.MIN_MAX | IS.STD_DEV
> stats = IS.getStatistics(ip, options, imp.getCalibration())
> boxW = 409
> boxH = 409
> imageW = 4096
> imageH = 4096
> i = 0
>
> nwidth = 4096
> nheight = 4096
> pixels = zeros('f', nwidth * nheight)
>
> for x in xrange(0,imageW,1):
> for y in xrange (0,imageH,1):
> imp.setRoi( y, x, boxW, boxH)
> # rimp = imp.getRoi
> stats = IS.getStatistics(ip, options, imp.getCalibration())
> # print x, ",", y, ",", stats.stdDev
> pixels[i] = stats.stdDev
> i += 1
>
> #generate image of local stdDev
> fp = FloatProcessor(nwidth, nheight, pixels, None)
> imp = ImagePlus("Local StdDev", fp)
> imp.show()
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

Saalfeld, Stephan

Re: Calculate local stdDev, using a sliding box

In reply to this post by Steve Wolf

Try

Plugins > Integral Image Filters > Standard Deviation

in Jython

from ij import IJ

IJ.run("Standard Deviation", "block_radius_x=204 block_radius_y=204")

or using the API directly and creating a new image as you did in your
code

from mpicbg.ij.integral import BlockStatistics
from ij import IJ, ImagePlus

boxW = 409
boxH = 409

imp = IJ.getImage()
ip = imp.getProcessor().duplicate().convertToFloatProcessor()
b = BlockStatistics(ip)
b.std(boxW/2, boxH/2)
imp = ImagePlus("Local StdDev", ip)
imp.show()

More info

http://fiji.sc/Integral_Image_Filters

In general, do not use an interpreted language to loop over pixel but
only to stitch together higher level pipelines.

Cheers,
Stephan

On Tue, 2016-02-09 at 11:43 -0500, Steve Wolf wrote:

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

Stefan Helfrich-2

Re: Calculate local stdDev, using a sliding box

In reply to this post by Steve Wolf

Dear Steve,

you could also use ImageJ Ops [1] to achieve your goal. Just copy the following to a new Jython script in the Script Editor and run it with an open image. It will ask for a windowSpan, which is the span of your sliding window (span==3 --> 7x7 window).

# @Dataset input
# @Integer windowSpan
# @OpService ops
# @UIService ui

from net.imglib2.algorithm.neighborhood import RectangleShape
from net.imglib2.outofbounds import OutOfBoundsMirrorFactory
from net.imglib2.outofbounds.OutOfBoundsMirrorFactory import Boundary

# Generate output image that has same dimensions as input
output = input.getImgPlus().factory().create(input, input.getImgPlus().firstElement())

# Apply mean filter using a RectangleShape with a span of windowSpan pixels
ops.filter().mean(output, input, RectangleShape(windowSpan, False), OutOfBoundsMirrorFactory(Boundary.SINGLE))

# Display the output
ui.show("Output", output)

All the best,
Stefan

[1] http://imagej.net/ImageJ_Ops

On Tue, 9 Feb 2016 11:43:57 -0500, Steve Wolf <[hidden email]> wrote:

>Attempting to calculate local stdDev, using a sliding box, on a 4096x4096 image. I’m using the following Python script, but it runs extremely slow, any suggestions on speeding it up?
>
>
>from ij import IJ
>from ij import ImagePlus
>from ij.process import FloatProcessor
>from array import zeros
>from ij.process import ImageStatistics as IS
>
>imp = IJ.getImage()
>ip = imp.getProcessor()
>
>options = IS.MEAN | IS.MEDIAN | IS.MIN_MAX | IS.STD_DEV
>stats = IS.getStatistics(ip, options, imp.getCalibration())
>boxW = 409
>boxH = 409
>imageW = 4096
>imageH = 4096
>i = 0
>
>nwidth = 4096
>nheight = 4096
>pixels = zeros('f', nwidth * nheight)
>
>for x in xrange(0,imageW,1):
> for y in xrange (0,imageH,1):
> imp.setRoi( y, x, boxW, boxH)
># rimp = imp.getRoi
> stats = IS.getStatistics(ip, options, imp.getCalibration())
># print x, ",", y, ",", stats.stdDev
> pixels[i] = stats.stdDev
> i += 1
>
>#generate image of local stdDev
>fp = FloatProcessor(nwidth, nheight, pixels, None)
>imp = ImagePlus("Local StdDev", fp)
>imp.show()
>
>--
>ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

Steve Wolf

Re: Calculate local stdDev, using a sliding box

In reply to this post by Steve Wolf

Thank you. I hadn't used Fiji before, just the standard ImageJ product. I didn't realize the Integral Image Filters existed either, I was simply able to use that. However, I did find that it doesn't seem to work on Floats; I converted the image to 16bit and processed the image.

Thanks again for the help!

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

Saalfeld, Stephan

Re: Calculate local stdDev, using a sliding box

It works on floats, you will have to adjust the visible min and max
range (brightness and contrast) to see the result.

Cheers,
Stephan

On Wed, 2016-02-10 at 13:17 -0500, Steve Wolf wrote:
> Thank you. I hadn't used Fiji before, just the standard ImageJ product. I didn't realize the Integral Image Filters existed either, I was simply able to use that. However, I did find that it doesn't seem to work on Floats; I converted the image to 16bit and processed the image.
>
> Thanks again for the help!
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html