pro SMEI_cr_removal, files, nbuffer=nbuffer,	$
	sigmathreshold=SigmaThreshold, minthreshold=MinThreshold, remove=Remove,	$
	silent=silent

@compile_opt.pro		; On error, return to caller

if n_elements(files) ne 0 then begin
	files = file_search(files)
	if files[0] eq '' then destroyvar, files
endif

if n_elements(files) eq 0 then begin
	files = dialog_pickfile(filter='*.nic')
	if files[0] eq '' then return
endif

if n_elements(nbuffer) eq 0 then nbuffer = 10
nfiles = n_elements(files)

if nfiles le nbuffer then begin
	message, /info, 'not enough files'
	return
endif

nstack = 2*nbuffer+1

fov = SMEI_camera(/get_structure)
stat = bin_read(files[0], tmp, /dim)			; Update image size
fov.nsize = tmp

cubedim = [fov.nsize, nfiles]

Bad = BadValue(0.0)

; The img array has room for nstack = 2*nbuffer images.
; Fill in 1st half of the buffer (0,..,nbuffer) with files 0,..,,nbuffer

img = fltarr(fov.nsize[0], fov.nsize[1], nstack, /nozero)
for i=0,nbuffer do begin
	stat = bin_read(files[i], tmp, silent=silent)
	img[*,*,i] = tmp
endfor
img[*,*,nbuffer+1:nstack-1] = Bad


rr = gridgen(fov.nsize, origin=fov.center)			; x-, y-coordinates for all pixel centers
rr = cv_coord(from_rect=rr, /to_polar)				; Polar coordinates

in_fov = Inside_Wedge(fov.limits, rr)
in_fov = where(in_fov, nfov)						; List of pixels with centers inside fov

rr = rr[*,in_fov]									; Polar coordinates of pixel centers inside fov

; rr[0:1,*,nbuffer] contains the 'nfov' unperturbed pixel locations
; rr is in polar coordinates.
; Except for rr[0:1,*,nbuffer] positions will not exactly match pixel centers.

rr = TMO_skymotion(rr, fov.axis, gridgen(nstack, range=nbuffer*[-1,1]))
rr = reform(rr, 2, nfov*nstack, /overwrite)			; Array[2,nfov*nstack]

; Note that in_fov_stack[*,nbuffer] = 1

in_fov_stack = Inside_Wedge(fov.limits, rr)			; Array[nfov*nstack], pixels inside fov; 0:no 1:yes

rr = cv_coord(from_polar=rr, /to_rect)				; Convert to rectangular

; Rectangular CCD coordinates for all 'nfov' pixels in stack
; Positions rr[0:1,i,*] follow pixel rr[0:1,i,nbuffer] across the sky.

rr = rr+SuperArray(fov.center, nfov*nstack, /trail)
rr = round(rr)										; Round values to nearest pixel centers


not_done = bytarr(nfov, nstack, /nozero)
not_done[*,0		:nbuffer ] = 1B					; Status for 'nstack' images
not_done[*,nbuffer+1:nstack-1] = 0B

; lstack runs from -nbuffer to +nbuffer for i=0,nstack-1
; i.e. lstack[nbuffer] = 0

lstack = round(gridgen(nstack, range=nbuffer*[-1,1]))

; Let i be the image currently begin processed
; Image i is stored in the image stack 'img' in position n_i = ls[nbuffer] = i mod nstack
; nbuffer preceding and following images are stored before and after position n_i:
;	preceding: n_i-nbuffer, .. , n_i-1 = ls[0:nbuffer-1]
;	following: n_i+1, .. , n_i+nbuffer = ls[nbuffer+1:nstack-1]
; (wrapping at the end of the image stack if necessary).

lstack = (nstack+lstack) mod nstack
lstack = shift(lstack, 1)

nfiles = 20

for i=0,nfiles-1 do begin

; Locations in buffer of preceding and following 'nbuffer' images in stack
; i.e. this is an index into the last dimension of 'img'.
; istack[nbuffer] gives the location of image 'i'

	lstack = shift(lstack, -1)
	n_i = lstack[nbuffer]						; Location of image i in stack (= i mod nstack)

; Check which pixels in image are not yet processed

	k_fov = where(not_done[*,n_i], m_fov)	; Array[m_fov], Index into nfov locations inside fov

	if m_fov ne 0 then begin				; Still unprocessed pixels present

; We collect m_fov unprocessed pixels from each of the 'nstack' images, i.e. m_fov*nstack elements
; Some of these will lie outside the field of view.
;
; k_fov_stack: m_fov*nstack elements; linear index into in_fov_stack, rr[0,*], rr[1,*]
; z_fov_stack: m_fov*nstack elements: index into the 3rd dimension of img (the frame number in the stack).

		k_fov_stack = k_fov#replicate(1L,nstack)+replicate(nfov,m_fov)#lindgen(nstack)
		z_fov_stack = replicate(1L,m_fov)#lstack

; Check which of the m_fov*nstack elements in k_fov_stack lie inside the field of view.
; k: index into k_fov_stack (also index into f_fov_stack array of size m_fov*nstack)

		k = where(in_fov_stack[k_fov_stack])		; Index into k_fov_stack

; Pixel values in pixels taking sky motion into account for m_fov unprocessed pixels.
; f_fov_stack[*,nbuffer] contains the pixels for current image i
; f_fov_stack[*,0:nbuffer-1] are from 'nbuffer' preceding images
; f_fov_stack[*,nbuffer+1:nstack-1] are from 'nbuffer' following images

		f_fov_stack = replicate(Bad, m_fov, nstack)	; Array for pixel values

; kk_in: subset of k_fov_stack above. These are the subset (< m_fov*nstack) of elements
; inside the field of view. This again is a linear index into in_fov_stack, rr[0,*], rr[1,*].

		tmp = k_fov_stack[k] 						; Subset of k_fov_stack
		f_fov_stack[k] = img[rr[0,tmp],rr[1,tmp],z_fov_stack[k]]

		ftot = round(total(finite(f_fov_stack),2))	; Array[m_fov]

; Only process pixels with more than 'nbuffer' points in the stack.
; tmp can be used as an index to k_fov, and reduces the number of pixels begin processed.

		tmp = where(ftot gt nbuffer, m_fov)

		if m_fov gt 0 then begin
			k_fov		= k_fov[tmp]				; Array[m_fov], Index into 'nfov' locations in fov
			k_fov_stack = k_fov_stack[tmp,*] 		; Array[m_fov,nstack]
			z_fov_stack = z_fov_stack[tmp,*]
			f_fov_stack = f_fov_stack[tmp,*]

; Update not_done status

			k   = where(in_fov_stack[k_fov_stack])	; Index into k_fov_stack

			tmp = k_fov_stack[k] 					; Subset of k_fov_stack

			done = bytarr(fov.nsize[0], fov.nsize[1], nstack)
			done[rr[0,tmp],rr[1,tmp],z_fov_stack[k]] = 1B
			done = reform(done, fov.nsize[0]*fov.nsize[1], nstack, /overwrite)
			done = done[in_fov,*]
			done = where(done)

			not_done[done] = 0B

; Check f_fov_stack for glitches
; The dummy 2nd dimension in f_fov_stack forces Find2DGlitch to treat it as a stack of
; 'nstack' images of size 'm_fov' by 1 pixels. Note that sumwidth and spotwidth both
; MUST be set to 1.

			f_fov_stack = reform(f_fov_stack, m_fov, 1, nstack, /overwrite)

			nCR = Find2DGlitch(f_fov_stack, sumwidth=1, spotwidth=1, /exclude,	$
				sigmathreshold=SigmaThreshold, minthreshold=MinThreshold, remove=Remove, loc=Loc)

			f_fov_stack = reform(f_fov_stack, m_fov, nstack, /overwrite)

; The return array Loc is a linear index into f_fov_stack and k_fov_stack

			if nCR gt 0 then begin
				Loc = ArrayLocation(Loc, dim=[m_fov,nstack])
				Loc = [rr[*,k_fov_stack[Loc[0,*],Loc[1,*]]], Loc[1,*]-nbuffer+i]

; Convert index into image cube with dimensions cubedim (nfov.nsize[0] x nfov.nsize[1] x nfiles)
; from 3D index to linear index.

				Loc = ArrayLocation(Loc, dim=cubedim, /onedim)

				if n_elements(Locations) eq 0 then	$
					Locations = Loc					$
				else								$
					Locations = [Location, Loc]

			endif

		endif
	endif

; Read next image i+nbuffer+1

	n   = i+nbuffer+1
	n_i = n mod nstack

	case n ge nfiles of
	0: begin
		stat = bin_read(files[n], tmp, silent=silent)
		img[*,*,n_i] = tmp
	end
	1: img[*,*,n_i] = Bad
	endcase

	not_done[*,n_i] = n lt nfiles

endfor

; Locations is an array of pixels into the full image cube of 'nfiles' images.
; We run this through GroupPixels to order the glitches into groups.

if n_elements(Locations) ne 0 then begin
	Locations = GroupPixels(Locations, dim=dim, range=[2,2,0],	$
		ngroup=nGroup, pGroup=pGroup, lgroup=lGroup)

; Still need to collect information about total ADUs in glitches
; Maybe also a histogram of glitches

endif


return  &  end