;+
; NAME: stddevarr
; PURPOSE: compute the standard deviation for each pixel in a cube
;
;NOTES:
; sd = stddevarr(array,/badval,alsobad=<>)
;
; Just as medarr returns the median at each pixel of an image cube,
; this routine returns the standard deviation for each pixel in a
; cube. This is useful when combining a stack of images if you want
; an estimate of the error at a given pixel.
;
; You can specify a badval which is used to mask out all pixels
; in the cube that have that value and only take the standard deviation
; of the remaining ones. This is useful if the images do not overlap
; completely and there are blank regions on some of the layers - just set
; the blank regions to BADVAL and set the badval keyword likewise.
;
; INPUTS:
; cube an image cube, M*N*K
; KEYWORDS:
; BADVAL set this equal to the pixel value of pixels to ignore.
; ALSOBAD in case you have a second bad pixel value.
; DIM= which dimension to take the stddev over. Default=3
; /DIVSQRT Divide by the square root of the number of good valuse at
; each pixel. In other words, return the standard deviation of
; the mean, instead of just the standard deviation of the data
; set.
; NOTE: If you want to compute the population standard deviation
; you should instead divide by the sqrt(n-1), but that's usually
; not what you want in image analysis contexts.
; /MEDIAN Compute Std. Dev. around the median, not the mean.
; OUTPUTS:
; returns an M*N array containing the stddev at each pixel
;
; HISTORY:
; Began 2002-06-21 14:51:12 by Marshall Perrin, based on
; errmap_neg.pro by Erik Rosolowsky
; 2003-07-16 Modifed keyword checking to allow '0' to be specified
; as a bad value.
; 2006-04-17 added dim= keyword; can now work on 2D images.
; 2006-06-01 Added recursive chunking of large cubes
;-
function stddevarr,cube,badval=badval,alsobad=alsobad,dim=dim,divsqrt=divsqrt, $
median=median
sz = size(cube)
if ~(keyword_set(dim)) then dim=3
if dim eq 3 then nz = sz[3] else nz=1
if dim eq 3 and n_elements(cube) gt 3e6 then begin
; break the problem down recursively for large arrays.
; This prevents it from using absolutely humongous amounts of
; memory. Cutoff at 3 million elements is chosen arbitrarily.
; and will certainly vary depending on the machine in question
; and properties of the array!
nx1 = sz[1]/3
; print,"Breaking down into smaller chunks: nx = "+strc(nx1)
return, [stddevarr( cube[0:nx1,*,*], badval=badval,alsobad=alsobad,dim=dim,divsqrt=divsqrt ), $
stddevarr( cube[nx1+1:2*nx1,*,*], badval=badval,alsobad=alsobad,dim=dim,divsqrt=divsqrt ), $
stddevarr( cube[2*nx1+1:*,*,*], badval=badval,alsobad=alsobad,dim=dim,divsqrt=divsqrt )]
endif
if (n_elements(badval) gt 0) then begin
if (n_elements(alsobad) gt 0)then begin
wgood=where((cube ne badval) and (cube ne alsobad) and finite(cube),goodcount)
endif else begin
wgood=where(cube ne badval and finite(cube),goodcount)
endelse
if goodcount le 0 then message,"No good pixels!"
endif else wgood = where(finite(cube),goodcount)
; are all pixels bad? if so, stddev = 0
if goodcount eq 0 then return, fltarr(sz[1],sz[2])
goodmask = cube*0
goodmask[wgood]=1
if (size(goodmask))[0] lt dim then return,fltarr(sz[1],sz[2])
tg= total(goodmask,DIM,/NaN)
if keyword_set(median) then means = median(cube*goodmask,dim=3)$
else means = total(cube*goodmask, DIM,/NaN)/tg
meancube = rebin(means,sz[1],sz[2],nz)
result = sqrt( total(((cube-meancube)^2)*goodmask, DIM,/NaN)/(tg-1) )
if keyword_set(divsqrt) then result /= sqrt(tg)
return, result
end