pro try, pedimg
on_error, 2

fil = filepath(root=getenv('SSW_SMEI_DAT'),subdir=['ccd','pleiades'],'pleiades.*')
fil = file_search(fil, count=n)
if n eq 0 then message, 'no files found'

firstfile = 2
n = 1;10

nview = 8 < (n-1)

nx = 384
stat = bin_read (fil[firstfile],skyimg,/aint,off=160,nx=nx)
if n gt 1 then begin
	skyimg = SuperArray(skyimg,n,/trail)
	for i=1,n-1 do begin
		stat = bin_read (fil[firstfile+i],tmp,/aint,off=160,nx=nx)
		skyimg[*,*,i] = tmp
	endfor
endif

skyimg = float(skyimg)

; The first three rows are ditched. The first two columns and the column separating
; the open and closed part of the CCD are also ditched.

pedimg  = skyimg[4:*,290:*,*]			; Covered part of the CCD
badpedx = [141]							; Hot pixel
badpedy = [196]
pedimg[badpedx,badpedy,*] = !VALUES.F_NAN

skyimg = skyimg[4:*,2:288,*]			; Open part of the CCD
badskyx = [224]							; Cold pixel
badskyy = [ 63]
skyimg[badskyx,badskyy,*] = !VALUES.F_NAN

;img = pedimg
img = skyimg

sz = size(img)
whatis, img

ped = fltarr(n,/nozero)
for i=0,n-1 do begin
	cmpimg = img[*,*,i]
	whatis, cmpimg
	ped[i] = median( cmpimg[where(finite(cmpimg))] )
	img[*,*,i] = img[*,*,i]-ped[i]
endfor

print, 'Pedestals:',strcompress(ped)

cmpimg = img

twin,/del,/all

sz = size(img)
i = max(img,loc,/nan)
print, 'Max at start: ',i,' at',[Loc mod sz[1],Loc/sz[1] mod sz[2],Loc/(sz[1]*sz[2])]

twin,/bl, /reversebw, xs=400,ys=400			; Input image, histogram
plot, histogram(img[*,*,nview],/nan) < 100

twin,/tl, /reversebw, xs=400,ys=400			; Input image, surface map
surface, img[*,*,nview],ax=5,az=20

view,in=img[*,*,nview],/stretch,/cl 		; Input image, greyscale

MaxStars = 60
Stars = replicate(-1L,n,MaxStars)

bTest = 0B
t1 = systime(1)
for i=0,n-1 do begin
	img0 = img[*,*,i]
	j = FindStars(img0,maxstars=MaxStars,threshold=10,	$
			starpix=StarPix,peakpix=PeakPix,starpos=StarPos,peakpos=PeakPos,btest=bTest)
	img[*,*,i] = img0
	Stars[i,0:n_elements(StarPos)-1] = StarPos
endfor

;openw, /get_lun, iu, 'd:\temp\stars.txt'
;printf, iu, format='('+strcompress(n,/rem)+'I6)',Stars
;free_lun, iu

t2 = systime(2)
print, t2-t1, (t2-t1)/n

;twin,/tr,/reversebw,xs=sz[1],ys=sz[2]
;px = StarPos mod sz[1]
;py = StarPos/sz[1]
;oplot, px[0:19], py[0:19], psym=2
;oplot, px[20:*], py[20:*], psym=1
;px = PeakPos mod sz[1]
;py = PeakPos/sz[1]
;oplot, px, py, psym=5
;return

i = max(img,Loc,/nan)
print, 'Max after star removal: ',i,' at',ArrayLocation(Loc,size=sz)

twin,/bc, /reversebw, xs=400,ys=400			; Input image, histogram
plot, histogram(img[*,*,nview],/nan) < 100

twin,/tc, /reversebw, xs=400,ys=400			; Input image, surface map
surface, img[*,*,nview],ax=5,az=20

view,in=img[*,*,nview],/stretch,/center			; Input image, greyscale

glitch = Find2DGlitch(img,sum=1,/exclude,loc=loc,sigma=5.0,spot=3,/remove)
;glitch = Find2DGlitch(img,sum=3,/exclude,loc=loc,sigma=5.0,spot=3,/remove,pixels=PeakPos)
t3 = systime(1)
print, t3-t2,t3-t1

if glitch gt 0 then	$
	print, '# glitches: ',n_elements(Loc), '  min: ',min(cmpimg[Loc]),'  max: ',max(cmpimg[Loc])

twin,/br, /reversebw, xs=400,ys=400
plot, histogram(img[*,*,nview],/nan) < 100

twin,/tr, /reversebw, xs=400,ys=400
surface, img[*,*,nview],ax=5

view,in=img[*,*,nview],/stretch,/cr

j = max(img,i,/nan)
print, 'New max:',j,' at', ArrayLocation(i,size=sz)

return  &  end