FUNCTION Find2DGlitch, Frames, sumwidth=SumWidth, spotwidth=SpotWidth,	$
		sigmathreshold=SigmaThreshold, minthreshold=MinThreshold,		$
		exclude=Exclude, remove=Remove, loc=Loc,				$
		group=group, ngroup=nGroup, pgroup=pGroup, lgroup=lGroup, fgroup=fGroup,	$
		quiet=quiet	;,	pixels=Pixels

;+
; NAME:
;	Find2DGlitch
; PURPOSE:
;	Find 'cosmic rays' in a set of 2D frames
; CATEGORY:
; CALLING SEQUENCE:
;	Result = Find2DGlitch(Frames [,sumwidth=SumWidth, /exclude,	$
;				sigmathreshold=SigmaThreshold, loc=Loc])
; INPUTS:
;	Frames			array[n,l,m]; type: any (float or double needed to use NaN option)
;						stack of 2D frames combined in 3D array; the last dimension counts
;						the number of frames (if it is absent then a dummy dimension of 1
;						is added).
;						Frame elements set to the value !VALUES.F_NAN or D_NAN are ignored
; OPTIONAL INPUT PARAMETERS:
;	SumWidth		scalar, type integer, no default
;						defines a box of SumWidth by SumWidth pixels used by procedure FrameMoments
;	SpotWidth		scalar, type integer, no default
;						defines a box of SpotWidth by SpotWidth pixels used by procedure CleanGlitchBox
;	sigmathreshold	scalar, any type, default: 4.
;						see PROCEDURE
;	minthreshold	scalar, any type, default: 0
;						see PROCEDURE
;	/exclude		if set, then each pixel in the Frames array is tested based on statistical
;						moments calculated omitting the pixel itself
;	/remove			if set the glitches are replaced by average values calculated
;						after removal of the glitches
;	/group			activates grouping of pixels that are close together.
; OUTPUTS:
;	Result			scalar, type long integer
;						number of glitches found
; OPTIONAL OUTPUT PARAMETERS:
;	loc=Loc			1-dim array [*], type long integer
;						location of the glitches
;						glitch positions are identified by a one-dimensional array index
;	/quiet			if set, the list of glitches found is not printed to the screen
;
;	ngroup=ngroup
;	pgroup=pgroup
;	lgroup=lgroup	unmodified return arguments from a call to the GroupPixels procedure
;	fgroup=fgroup	contains total # adus in each group (after subtraction of underlying average)
; INCLUDE:
	@compile_opt.pro	; On error, return to caller
; CALLS:
;	FrameMoments, ArrayLocation, GroupPixels, CleanGlitchBox, BadValue
; PROCEDURE:
;	If there are nT frames of dimension nX by nY then Frames has dimensions nX,nY,nT
;
;	For each of the nX*nY pixels an average and standard deviation is calculated by averaging
;	all nT boxes of SumWidth^2 pixels (total of SumWidth^2*nT pixels).
;
;	A 'cosmic ray' is defined as a pixel where the frame value is more than SigmaThreshold
;	standard deviations above the average.
;	If MinThreshold is set than the frame value must be at least MinThreshold above the
;	average to count as a glitch. This second test can be used to avoid that lots of points
;	are flagged when the standard deviation becomes real small.
; MODIFICATION HISTORY:
;	OCT-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, quiet, /key
InitVar, group, /key

szin = size(Frames)
IF szin[0] NE 2 AND szin[0] NE 3 THEN RETURN, 0	; Only accept 2-dim and 3-dim arrays

Mask = where(finite(Frames))					; Make sure there are valid pixels
IF Mask[0] EQ -1 THEN RETURN, 0					; If not, return zero ('no glitches found')

CASE szin[0] EQ 2 OF							; Add dummy dimension of 1 to Frames
0: sz = szin
1: BEGIN
	Frames = reform(Frames, [szin[1:2],1], /overwrite)
	sz = size(Frames)
END
ENDCASE

InitVar, SigmaThreshold, 4.0
InitVar, MinThreshold  , 0.0
InitVar, Remove, /key
												; Calculate statistics
FrameMoments, Frames, sumwidth=SumWidth, exclude=Exclude, Average, Sigma;, pixels=Pixels

; Reduce potential glitches to actual glitches.
; 'Mask' lists all pixels with finite values as indices into Frames
; 'Loc' will list the glitches as an index into Mask

nLoc =  0
Loc  = -1

nGroup =  0
pGroup = -1
lGroup =  0

Mask = where(finite(Sigma))
IF Mask[0] NE -1 THEN BEGIN
	Loc = where(Frames[Mask]-Average[Mask] GT (SigmaThreshold*Sigma[Mask] > MinThreshold),nLoc)
	;Loc = where(Sigma[Mask] GT 0. and Frames[Mask]-Average[Mask] GT SigmaThreshold*Sigma[Mask],nLoc)

	IF nLoc GT 0 THEN BEGIN						; Glitches found
		Loc = Mask[Loc]							; Convert from index into Mask to index into Frames
												; Check neighbouring pixels in same frame where glitch is found
		nLoc = CleanGlitchBox(Frames,Loc,spotwidth=SpotWidth,frac=frac)

		IF group THEN BEGIN
			Loc = GroupPixels(Loc, sizeframe=sz, range=[2,2,0],	$
					ngroup=nGroup, pGroup=pGroup, lgroup=lGroup)
			;iP = ArrayLocation(Loc[pGroup],size=sz)
			;fP = Frames[Loc[pGroup]]

			fGroup = fltarr(nGroup, /nozero)
			fmin   = fltarr(nGroup, /nozero)
			fmax   = fltarr(nGroup, /nozero)

			FOR i=0,nGroup-1 DO BEGIN
				fval = Loc[pGroup[i]:pGroup[i]+lGroup[i]-1]
				fval = Frames[fval]-Average[fval]
				fmin[i] = min(fval, max=tmp)
				fmax[i] = tmp
				fGroup[i] = total(fval)
			ENDFOR

			iP = ArrayLocation(Loc[pGroup],size=sz)

			FOR i=0,(1-quiet)*nGroup-1 DO print, 'CR  '  +strcompress(i+1		)+				$
				' at ['		+strcompress(iP[0,i])+','	 +strcompress(iP[1,i]	)+']'+			$
				' of frame'	+strcompress(iP[2,i])+' ('	 +strcompress(lGroup[i]	)+' pixels)'+	$
				', min'		+strcompress(fmin[i])+', max'+strcompress(fmax[i]	)+', tot'+strcompress(fGroup[i])

		ENDIF

		IF Remove THEN BEGIN					; Remove glitches
			Frames[Loc] = BadValue(Frames) 		; Recalculate average
			FrameMoments, Frames, sumwidth=SumWidth, Average;,pixels=Loc
			Frames[Loc] = Average[Loc]			; Set glitches to new average
		ENDIF

	ENDIF

ENDIF

IF szin[0] EQ 2 THEN Frames = reform(Frames, /overwrite)

RETURN, nLoc  &  END