;+
; NAME: ircal_fixpix
; PURPOSE:
; 	This is a front-end to ns_fixpix.pro which knows about the
; 	IRCAL bad corners and preserves them in the bad pixel mask.
;
; INPUTS:
; KEYWORDS:
; 		/maskcorners	if set, the bad corners as found in the input
;						will be set to NaN in the output as well.
;
;		/stars			Ignore areas around bright stars. This prevents
;						identifying undersampled PSFs as bad pixels
;						WARNING! This has not really been well tested yet.
; OUTPUTS:
;
; HISTORY:
; 	Began 2003-11-25 18:19:40 by Marshall Perrin 
; 	2006-05-25	Added /stars option
; 	2006-06-09	reworked /stars option to use label_regions, and
; 				other algorithm improvements. Added starmask=.
;-

PRO ircal_fixpix,imgs,badpix,maskcorners=maskcorners,fixnans=fixnans,$
	stars=stars,iter=iter,starmask=starmask


  	sz = size(imgs)
    ; keep track of the bad regions in the corners.
    i1 = imgs[*,*,0]
    wnan = where(finite(i1,/NaN),nanct) & if nanct gt 0 then i1[wnan]=-1e6 ; can't use search2d on NaNs directly
    badleft = search2d(i1,0,0, i1[0,0],i1[0,0])
    badright = search2d(i1,255,0, i1[255,0],i1[255,0])


	if keyword_set(stars) then begin
		; check for the cores of bright stars and DON'T let them be marked as
		; bad pixels
		; better algorithm: find stars, then shrink each region to only
		; have a few pixels being ignored in it?
		message,/info, "Locating stellar peaks to mask out during "
		message,/info, "bad pixel identification."

		if sz[0] eq 2 then nz = 1 else nz = sz[3]

		;s = smooth(imgs,[7,7,1],/nan,/edge)
		; try median rather than mean. This will further reduce the
		; effect of hot pixels.
		medsize=7
		s = fltarr(sz[1],sz[2],nz) & for i=0,nz-1 do s[0,0,i] = median(imgs[*,*,i],medsize)
		s[*,0:medsize,*] =0
		s[*,sz[2]-1-medsize:*] =0
		;s2 = median(imgs,[7,7,1]);,/edge)

		;stop
		if arg_present(starmask) then starmask = bytarr(sz[1],sz[2],nz)

		for i=0L,nz-1 do begin
			myim = imgs[*,*,i]
			mmm,myim,s1,s2
			thresh = s1+5*s2
			mymask = s[*,*,i] gt thresh
			if total(mymask) eq 0 then continue ; no stars in this slice
			blobs = label_region(mymask)
			; Get population and members of each blob:  
			h = HISTOGRAM(blobs, REVERSE_INDICES=r)  
			; ignore region 0, that's the background
			FOR j=1, N_ELEMENTS(h)-1 DO BEGIN  
				;Find subscripts of members of region i.  
				p = r[r[j]:r[j+1]-1] 
				region = myim[p]
				rmax = max(region,/nan)
				;if ~finite(m) then continue; ignore all NaN regions
				whigh = where(region gt rmax/2,highct)
				; find pixels which are both
				;  1. in this region
				;  2. greater than 50% of the region maximum
				;wstar = cmset_op(p,"AND", where(myim gt rmax/2),count=count)
				;wstar = p[whigh]
				; could make this a bit faster by just storing the indices
				; into the final wstar array directly?
				if highct gt 0 then stackpush,wstar,p[whigh] + sz[1]*sz[2]*i  ; starmask[ p[whigh] + sz[1]*sz[2]*i ]=1
			endfor
			;imdisp,mymask,/stop
			;starmask[*,*,i] = mymask

		endfor
		;wstar = where(starmask, starct)
		imgs0 = imgs
	endif
		
	; now fix pixels
    message,/info,"Running ns_fixpix to clean up bad pixels"
   	imgs =ns_fixpix(imgs,/NaN,iter=iter)

	if keyword_set(stars) then begin
		; put back in any stellar cores
		;imdisp,s gt thresh,/axis,title="Regions w/ stars: not fixing pixels here"
		starct = n_elements(wstar)
		if starct gt 0 then begin
			imgs[wstar] = imgs0[wstar]
			if arg_present(starmask) then starmask[wstar]=1
		endif
		
	endif
	
    ; some NaNs may survive the cleaning:
	if keyword_set(fixnans) then begin
		wnan = where(imgs ne imgs,nanct)
		if nanct gt 0 then imgs[wnan] = median(imgs)
	endif
    i1 = imgs[*,*,0]
    badpix = i1 eq i1
    ; also throw away the ircal bad corners, which should remain masked out.
    badpix[badleft] = 0
    badpix[badright]=0
    ;badpix[255,*]=0 ; last row seems bad too.

	if keyword_set(maskcorners) then begin
		; mask the corners to NaNs in the output as well
		wbad = where(badpix eq 0,badct)
		npix = sz[1]*sz[2]
		if badct gt 0 then for i=0,sz[3]-1 do imgs[wbad+npix*i]=!values.f_nan
	endif
	
	ncorners = n_elements(badleft)+n_elements(badright)
    message,/info,"    "+strc(fix(total(badpix eq 0)))+" bad pixels remain after cleaning"
	message,/info,"    including "+strc(ncorners)+" in the bad corners."

end