pro FrameMoments, Frames, Average, Sigma, sumwidth=SumWidth, exclude=Exclude, pixels=Pixels

@compile_opt.pro			; On error, return to caller

;+
; NAME:
;	FrameMoments
; PURPOSE:
;	Calculate average and standard deviations sigma
; CATEGORY:
;	Statistics
; CALLING SEQUENCE:
;	FrameMoments, Frames, sumwidth=SumWidth, Average, Sigma
; INPUTS:
;	Frames		3D-array, any type, but should be float to use !VALUES.F_NAN
;					to indicate missing pixels.
;					stack of 2D frames combined in 3D array; the last dimension counts
;					the number of frames
;					Frame elements set to the value !VALUES.F_NAN are ignored
; OPTIONAL INPUT PARAMETERS:
;	sumWidth	scalar, type integer, default: 3. (single frame) or 1. (multiple frames)
;					defines a box of sumwidth by sumwidth pixels
;					(should be ODD; argument is passed to the IDL 'smooth' function).
;	/exclude	if set, the statistical moment for each pixel in the Frames
;					array is calculating, excluding the pixel itself
; OUTPUTS:
;	Average, Sigma
;				arrays of same structure as Frames, type float
;					mean, standard deviation
;					If /exclude is not set all frames in the Average and Sigma
;					arrays will be identical
; CALLS:
;	ArrayLocation, SuperArray
; PROCEDURE:
;	For each pixel the statistical moments are calculated by summing over
;	sumwidth^2 x (# frames) pixels.
;	The summing needed to calculate the statistical moments is done using
;	'total' (to sum over frames) and 'smooth' (to sum over sumwidth^2 pixels in each frame)
; MODIFICATION HISTORY:
;	OCT-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

sz = size(Frames)
if sz[0] ne 2 and sz[0] ne 3 then message, 'Argument Frames must be 2-dim or 3-dim array'

if n_elements(SumWidth) eq 0 then SumWidth = 3 else SumWidth = fix(SumWidth[0])
width = (SumWidth-1)/2

Exclude = keyword_set(Exclude)

; Calculate sums over SumWidth^2*sz[3] pixels.
; Pixels set to !VALUES.F_NAN are ignored in the summing
; The array type is float to avoid truncation problems with the 'smooth' function

N = finite(Frames)
if (where(N))[0] eq -1 then message, 'no valid input data in argument Frames

; I think the following line is unnecessary, but I'm not quite sure.
; Since it quadruples the amount of memory needed for N it might be worth checking out.

N = float(N)							; 0=Nan, 1=Valid number

F = float(Frames)						; Contains NaN of Inf where N = 0

nP = n_elements(Pixels)

if nP ne 0 then begin					; Check individual pixels
	SN = fltarr(nP)
	S1 = fltarr(nP)
	S2 = fltarr(nP)

	pp = ArrayLocation(Pixels,size=sz)	; X,Y,Z-coordinates

	for iP=0,nP-1 do begin				; Loop over all pixels: ..
		i1 = (pp[0,iP]-width) > 0			; .. sum over frames and SumWidth*SumWidth pixels
		i2 = (pp[0,iP]+width) < (sz[1]-1)
		j1 = (pp[1,iP]-width) > 0
		j2 = (pp[1,iP]+width) < (sz[2]-1)
		SN[iP] = total(N[i1:i2,j1:j2,*])
		i  = F[i1:i2,j1:j2,*]
		S1[iP] = total(i  , /nan)
		S2[iP] = total(i*i, /nan)
	endfor

	if sz[0] eq 3 then begin			; If more than one frame: replicate sz[3] times
		SN = SuperArray(SN, sz[3], /trail)
		S1 = SuperArray(S1, sz[3], /trail)
		S2 = SuperArray(S2, sz[3], /trail)
	endif

	iP = Pixels
	if sz[0] eq 3 then begin		; Position of pixels in all frames
		iP = Pixels# replicate(1L,sz[3])+replicate(1L,nP)#(sz[1]*sz[2]*lindgen(sz[3]))
		iP = reform(iP, nP*sz[3], /overwrite)
	endif

	if Exclude then begin				; Exclude the pixel itself
		N = N[iP]
		F = F[iP]

		i = where(N eq 0)
		if i[0] ne -1 then F[i] = 0.	; Replace NaN and Inf by 0.

		SN = SN-N						; Subtract 1 where Frames value finite
		S1 = S1-F
		S2 = S2-F*F
	endif

	Average = replicate(!values.f_nan, sz[sz[0]+2])
	Average = reform(Average, sz[1:sz[0]], /overwrite)
	Average[iP] = S1/SN

	Sigma = replicate(!values.f_nan, sz[sz[0]+2])
	Sigma = reform(Sigma, sz[1:sz[0]], /overwrite)
	Sigma[iP] = sqrt((SN*S2-S1*S1)/(SN*(SN-1)))

endif else begin
	SN = N								; Accumulate # points; 0 or 1 (no NaN's)
	S1 = F								; Accumulate sum X; NaN's are kept for 'smooth'
	S2 = F*F							; Accumulate sum X^2

	if sz[0] eq 3 and sz[3] gt 1 then begin	; Sum over all frames
		SN = total(SN, 3)
		S1 = total(S1, 3, /nan)
		S2 = total(S2, 3, /nan)
	endif

	if width gt 0 then begin			; Sum over SumWidth^2 pixels

; 'Smooth' behaves unexpected if both keywords /nan and /edge are used at the same time.
; If there is an NaN value away from the edge, the missing point is apparently filled in
; by some interpolation scheme before the average over SumWidth^2 pixels is calculated,
; i.e. the NaN is not actually ignored. This means that this keyword combination can't be
; used to calculate the sums needed here. To be able to get the correct result with only
; the /nan keyword it is necessary to create a bigger array by 2*width.

		VN = replicate(0.			,sz[1]+2*width,sz[2]+2*width)
		V1 = replicate(!values.f_nan,sz[1]+2*width,sz[2]+2*width)
		V2 = V1

		w1 = sz[1]+width-1
		w2 = sz[2]+width-1
		VN[width:w1,width:w2] = SN						; Copy S* into center of the center of V*
		V1[width:w1,width:w2] = S1
		V2[width:w1,width:w2] = S2

		j = SumWidth*SumWidth							; # pixels in box
		i = smooth(float(VN ne 0), SumWidth)*j			; # finite values in box

		VN = smooth(VN, SumWidth      )*j
		V1 = smooth(V1, SumWidth, /nan)*i
		V2 = smooth(V2, SumWidth, /nan)*i

		SN = VN[width:w1,width:w2]						; Copy the center from V* into S*
		S1 = V1[width:w1,width:w2]
		S2 = V2[width:w1,width:w2]
	endif

	if sz[0] eq 3 then begin							; Replicate sz[3] times to match structure of Frames
		SN = SuperArray(SN, sz[3], /trail)
		S1 = SuperArray(S1, sz[3], /trail)
		S2 = SuperArray(S2, sz[3], /trail)
	endif

	if Exclude then begin					; Remove one pixel from sums
		SN = SN-N							; Subtract 1 where Frames value finite

		i = where(N eq 0)
		if i[0] ne -1 then F[i] = 0.		; Frames value where finite; 0 where NaN

		S1 = S1-F
		S2 = S2-F*F
	endif

	Average = S1/SN
	Sigma = sqrt((SN*S2-S1*S1)/(SN*(SN-1)))

endelse

;x0 = 285  &  y0 = 258
x0 = 258  &  y0 = 285
dx = 2  &  dy = 2

;print, '-----'
;print, S1 [x0-dx:x0+dx,y0-dy:(y0+dy) < (sz[2]-1)]
;print
;print, SN [x0-dx:x0+dx,y0-dy:(y0+dy) < (sz[2]-1)]
;print
;print, frames [x0-dx:x0+dx,y0-dy:(y0+dy) < (sz[2]-1)]
;print
;print, average[x0-dx:x0+dx,y0-dy:(y0+dy) < (sz[2]-1)]
;print
;print, sigma  [x0-dx:x0+dx,y0-dy:(y0+dy) < (sz[2]-1)]

return  &  end