;+
; NAME:
;	smei_frm_base
; PURPOSE:
;	Calculate pedestal and dark current for frame
; CATEGORY:
;	smei/camera/idl/misc
; CALLING SEQUENCE:
	PRO smei_frm_base, hdr, img, silent=silent,						$
		cam_full=cam_full, cam_ok=cam_ok, init=init, reset=reset,	$
		clear=clear, frm=frm, step=step, status=status, 			$
		nfill_control=nfill_control, nbad_control=nbad_control
; INPUTS:
;	hdr 		array[1]; type: structure
;					frame header
;	img 		array; type: integer
;					frame data
; OPTIONAL INPUT PARAMETERS:
;	nfill_control=nfill_control
;				scalar; type: integer; default: 10
;	/init		initialize all camera/mode combinations
;	/reset		re-initialize the current camera/mode combination only
;	/clear		clear all variables stored internally in common block
;	silent		=0: display all messages
;				=1: suppress non-essential messages
;				=2: suppress all messages
; OUTPUTS:
;	hdr 		array[1]; type: structure
;					frame header with a number of values filled, in particular
;					hdr.pedestal, hdr.dark current, hdr.squares, hdr.centerpix
;	img 		array; type: float
;					image with pedestal and dark current subtracted
;					Any negative elements after subtraction are set to zero
;					(these will mostly be pixels in the semi-covered columns
;					(5 and 315 for mode 2 data)
; OPTIONAL OUTPUT PARAMETERS:
;	cam_ok=cam_ok
;				scalar; type: byte
;				0: failed to determine pedestal and/or dark current
;				1: pedestal and dark current determined successfully
;	nfill_control=nfill_control
;				scalar; type: integer
;					# good frames used to build reference base values
;	cam_full=cam_full
;				0: still busy accumulating reference base data (based on 
;					nfill_control preceeding good frames) for current
;					camera/mode combination
;				1: reference data for preceeding nfill_control good frames
;					available
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
	@smei_roi_mask.pro		; SMEI masks
; COMMON BLOCKS:
	common smei_frm_base_save, ped_threshold, drk_threshold,			$
		ped_mean_init, drk_mean_init, ped_sigma_init, drk_sigma_init,  	$
		nfull, nbads, ifull, bfull, nbad, badv,	$
		peds_mean, peds_sigma, drks_mean, drks_sigma,$
		ped_mean , ped_sigma , drk_mean , drk_sigma , $
		drk_nout_init, drk_dout_init, drks_nout, drk_nout, drk_dout,	$
		drk_lr_ratio
; CALLS:
;	InitVar, IsType, TimeGet, TimeUnit, BadValue, smei_property
;	destroyvar
; PROCEDURE:
;	
; MODIFICATION HISTORY:
;	MAY-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, silent, 0 

InitVar, init , /key
InitVar, reset, /key
InitVar, clear, /key
InitVar, step , 1

IF clear THEN BEGIN 				; Clear all the common block variables, and return

	destroyvar, ped_threshold, drk_threshold
	destroyvar, ped_mean_init, drk_mean_init, ped_sigma_init, drk_sigma_init
	destroyvar, nfull, ifull, bfull  , nbad, badv
	destroyvar, peds_mean, peds_sigma, drks_mean, drks_sigma
 	destroyvar, ped_mean , ped_sigma , drk_mean, drk_sigma

	destroyvar, drk_nout_init, drk_dout_init, drks_nout, drk_nout, drk_dout
	destroyvar, drk_lr_ratio

	RETURN

ENDIF ELSE IF init THEN BEGIN		; Initialize common block variables, and return

	; Starting values for 'nominal' pedestal and dark current for each camera
	; These are used only until better values have been calculated from
	; the frames themselves.

	; Total # ped and dark pixels are 2048,512,64 in mode 0,1,2
	; 1st dim is camera, 2nd dim is mode

	;ped_strt = [1000.5d0,1004.0d0,1063.0d0]
	;drk_strt = [1.5d0,1.5d0,13.0d0]

	; Pedestal sigma for individual frames:
	;			Cam1		Cam2		Cam3
	; Mode 0	2.3			2.2			2.4
	; Mode 1	1.2			1.3			1.4
	; Mode 2	0.6-0.9		0.6-1.0		1.6


	; The following initializations don't need to be very accurate.
	; They are needed to get started, and are updated be values derived
	; from the data being processed.

	; Mean pedestal (close to 1000) calculated for all pixels in single frame

					; Cam 1   Cam 2  Cam 3

	ped_mean_init = 1000.0d0+	$
					[[0.5d0, 4.3d0, 63.4d0],	$	; Mode 0
					 [0.5d0, 4.3d0, 63.4d0],	$	; Mode 1
					 [0.5d0, 4.3d0, 63.4d0] 	]	; Mode 2

	; Mean dark current calculated for all pixels in single frame

					; Cam 1   Cam 2  Cam 3

	drk_mean_init = [[1.4d0, 1.4d0, 12.4d0],	$	; Mode 0
					 [1.4d0, 1.4d0, 12.4d0],	$	; Mode 1
					 [1.4d0, 1.4d0, 12.4d0]		] 	; Mode 2

	; Sigma for pedestal for all pixels in a single frame

					; Cam 1   Cam 2  Cam 3

	ped_sigma_init = [[2.30d0,2.2d0,2.4d0],		$	; Mode 0
					  [1.20d0,1.3d0,1.4d0],		$	; Mode 1
					  [0.55d0,1.0d0,1.6d0]		]	; Mode 2

	; Sigma for pedestal+dark current for all pixels in a single frame
	; (i.e. NOT just for the dark current alone)

					; Cam 1   Cam 2  Cam 3

	drk_sigma_init = [[2.6d0,2.5d0,10.5d0],		$	; Mode 0
					  [1.2d0,1.3d0,1.4d0],		$	; Mode 1
					  [1.4d0,1.0d0,1.6d0]		]	; Mode 2

	; Sigma for pedestal+dark current for all pixels in a single frame

					; Cam 1   Cam 2  Cam 3

	drk_nout_init  = [[50.0d0,50.0d0,120.0d0],	$	; Mode 0
					  [ 1.2d0, 1.3d0, 1.4d0],	$	; Mode 1
					  [ 0.9d0, 3.0d0, 1.6d0]	]	; Mode 2

	drk_dout_init  = [[50.0d0,50.0d0,120.0d0],	$	; Mode 0
					  [ 1.2d0, 1.3d0, 1.4d0],	$	; Mode 1
					  [ 0.9d0, 0.5d0, 1.6d0]	]	; Mode 2

	drk_lr_ratio	= [3.0d0,1.5d0,0.75d0]			; Mode 0,1,2

	; The next four get updated with info from the previous nfull frames

	ped_mean  = ped_mean_init
	drk_mean  = drk_mean_init
	ped_sigma = ped_sigma_init
	drk_sigma = drk_sigma_init

	drk_nout  = drk_nout_init
	drk_dout  = drk_dout_init

	; The next four thresholds never change

					; Cam 1   Cam 2  Cam 3

	ped_threshold = [[1850,1850,1850],	$			; Mode 0
					 [ 462, 462, 462],	$			; Mode 1
					 [  58,  58, 58]	]			; Mode 2

					; Cam 1   Cam 2  Cam 3

	drk_threshold = [[1000,1000,1000],	$			; Mode 0
					 [ 250, 250, 250],	$			; Mode 1
					 [  31,  31,  31]	]			; Mode 2

	InitVar, nfill_control, 10
	nfull = nfill_control		; # previous frames used for update statistical info

	InitVar, nbad_control, 2*nfull
	nbads = nbad_control

	ifull = replicate(-1,3,3)	; Counter for # previous frames processed 
	bfull = bytarr(3,3) 		; Flag set to true when nfull frames previous frames have been processed
	nbad  = lonarr(3,3)

	badv  = BadValue(ped_mean)

	; The next four contain the information for up to nfull previous frames
	; The averaging is over all pedestal or dark current pixels in a frame

	peds_mean  = replicate(badv,nfull,3,3)	; Pedestal mean
	drks_mean  = peds_mean 		; Dark current mean
	peds_sigma = peds_mean		; Pedestal standard deviation
	drks_sigma = drks_mean		; Dark current standard deviation

	drks_nout  = drks_sigma

	RETURN

ENDIF

cam  = hdr.camera-1
mode = hdr.mode

IF reset OR nbad[cam,mode] GT nbads THEN BEGIN

	message, /info, 'resetting camera'+strcompress(cam+1)+' mode'+strcompress(mode)

	; Reinitialize the current camera/mode combination.
	; i.e. erase the information for the previous nfull frames, and
	; start collecting new information from scratch.

	ped_mean [cam,mode] = ped_mean_init [cam,mode]
	drk_mean [cam,mode] = drk_mean_init [cam,mode]
	ped_sigma[cam,mode] = ped_sigma_init[cam,mode]
	drk_sigma[cam,mode] = ped_sigma_init[cam,mode] 

	drk_nout [cam,mode] = drk_nout_init [cam,mode]
	drk_dout [cam,mode] = drk_dout_init [cam,mode]

	ifull[cam,mode] = -1
	bfull[cam,mode] =  0
	nbad [cam,mode] =  0

	peds_mean [*,cam,mode] = badv
	drks_mean [*,cam,mode] = badv
	peds_sigma[*,cam,mode] = badv
	drks_sigma[*,cam,mode] = badv

	drks_nout [*,cam,mode] = badv

ENDIF


status = TimeGet(smei_property(hdr,/time),/_ydoy,upto=TimeUnit(/sec))
IF IsType(frm, /defined) THEN status = status+' ('+strcompress(frm, /rem)+')'
status = status+' c'+strcompress(cam+1,/rem)+'m'+strcompress(mode,/rem)

;print, 'SS flatfield: ', smei_property(hdr,/onboard_ff_enabled), min(img[*roi_bad[mode]]), max(img[*roi_bad[mode]])

; Calculate pedestal
; - Ignore values of zero or 512 in the pedestal columns
;	The last one to three pixels in a frame are sometimes contain
;	0 or 512. We explicitly exclude those (this is probably redundant
;	since these are also caught by the floor value)
; - Ignore values larger than 'ceiling'.
;   ('ceiling' is based on the pedestal for the previous images for this
;	camera, plus a constant positive offset)
; - Ignore value less than 'floor'.

; This test mainly catches 'crazy frames'

head_room = ([1.0d0,50000d0/65535d0])[smei_property(hdr,/onboard_ff_enabled)]
exclude = round([0,512]/head_room)

; Effectively keep ped_sigma above unity. Since the values in each frame are integers
; problems could arise if ped_sigma drops so far below unity that the range
; [floor, ceiling] contains only one possible frame value.

floor		= ped_mean[cam,mode]-3*(ped_sigma[cam,mode] > 1)
ceiling 	= ped_mean[cam,mode]+5*(ped_sigma[cam,mode] > 1)
threshold	= ped_threshold[cam,mode]

frm_ped_mean = smei_property(img, /average, /pedestal, exclude=exclude,	$
	floor=floor, ceiling=ceiling, threshold=threshold, rejected=rejected,	$
	silent=silent, sigma=frm_ped_sigma, stdev=frm_ped_stdev, count=frm_ped_nout)

CASE finite(frm_ped_mean) OF
0: BEGIN
	status = status+', count'+strcompress(frm_ped_nout[2])+'<'+strcompress(threshold,/rem)+	$
		', bad ped'+strcompress(rejected)+						$
		', nominal='+strcompress(ped_mean[cam,mode],/rem) +		$
		' ['+strjoin(strcompress([floor,ceiling],/rem),',')+']'
	message, /info, status

	frm_drk_mean  = badv
	frm_drk_sigma = badv
	frm_drk_stdev = badv
	frm_drk_nout  = [0L,0L]

	nbad[cam,mode] = nbad[cam,mode]+1

	IF NOT bfull[cam,mode] THEN BEGIN
		ped_mean [cam,mode] = smei_property(img, /average, /pedestal,	$
			sigma=frm_ped_sigma, stdev=frm_ped_stdev, silent=silent)
		ped_sigma[cam,mode] = frm_ped_sigma
	ENDIF
END
1: BEGIN

	; Calculate dark current, and subtract from image
	; This does not (yet?) take the slope of the dark current from row to row
	; into account
	; - Ignore values of zero in the dark current. The zeros in the pedestal
	;   pixels at the end of the frame (top row, last three columns) do not
	;   appear to extend into the dark current columns. We check for them
	;   anyway. NOTE: don't subtract the pedestal from the frame and calculate
	;   the dark current from the subtracted image. This screws up the zero
	;   exclusion (the subtraction might introduce some valid zeros).
	; - Ignore values larger than 'ceiling'
	;   ('ceiling' is based on the dark current for the previous frame for this
	;   camera, plus a constant positive offset)
	; - Ignore values less than 'floor'. For mode=2 an occasional funny values
	;   appears in the very last pixel. These can not all be explained by
	;   binning pixels with one or more pixels set to zero. 
	;	We never ever expect anything very far below the pedestal.

	floor     = (frm_ped_mean+drk_mean[cam,mode]-5*drk_sigma[cam,mode]) > (frm_ped_mean-3*ped_sigma[cam,mode]) 
	ceiling   =  frm_ped_mean+drk_mean[cam,mode]+5*drk_sigma[cam,mode]
	threshold = drk_threshold[cam,mode]

	frm_drk_mean = smei_property(img, /average, /dark_current,			$
		exclude=exclude, floor=floor, ceiling=ceiling, threshold=threshold,	$
		rejected=rejected, silent=silent, sigma=frm_drk_sigma,				$
		stdev=frm_drk_stdev, count=frm_drk_nout)


	CASE finite(frm_drk_mean) OF

	0: BEGIN

		status = status+', count'+strcompress(frm_drk_nout[2])+'<'+			$
			strcompress(threshold,/rem)+', ped'+strcompress(frm_ped_mean)+	$
			', bad drk'+string(rejected-frm_ped_mean,format='(F7.3)')

		nbad[cam,mode] = nbad[cam,mode]+1		; Count # sequential bad frames

	END

	1: BEGIN

		frm_drk_mean = frm_drk_mean-frm_ped_mean
		status = status+', ped'+strcompress(frm_ped_mean)+', drk'+strcompress(frm_drk_mean)

		CASE frm_drk_mean GT 0 OF
		0: BEGIN
		
			status = status+', negative, bad drk'
			frm_drk_mean = badv

		END
		1: BEGIN

			lr_ratio =	$
				smei_property(img, /average, /dark_current, /right,	$
					/sandwich, /high, exclude=exclude, floor=floor,		$
					ceiling=ceiling, silent=2)-							$
				smei_property(img, /average, /dark_current, /left ,	$
					/sandwich, /high, exclude=exclude, floor=floor,		$
					ceiling=ceiling, silent=2)

			CASE finite(lr_ratio) OF
			0: BEGIN
				status = status+', no L/R, bad drk'
				frm_drk_mean = badv
			END
			1: BEGIN
				lr_ratio = lr_ratio/frm_drk_mean
				IF abs(lr_ratio) GT drk_lr_ratio[mode] THEN BEGIN
					status = status+', L/R'+string(lr_ratio,format='(F4.1)')+	$
						'>'+string(drk_lr_ratio[mode],format='(F3.1)')+', bad drk'
					frm_drk_mean = badv
				ENDIF
			END
			ENDCASE

		END
		ENDCASE

		IF bfull[cam,mode] THEN BEGIN

			IF finite(frm_drk_mean) THEN BEGIN
				tmp = drk_nout[cam,mode]+20*(drk_dout[cam,mode] > 0.5)
				;print, 'checking for ',20*drk_dout[cam,mode], ' outliers'
				IF tmp GT 0 AND frm_drk_nout[1] GT tmp THEN BEGIN
					status = status+', outliers'+strcompress(frm_drk_nout[1])+'>'+strcompress(tmp)+', bad drk'
					frm_drk_mean = badv
				ENDIF
			ENDIF

			IF finite(frm_drk_mean) THEN BEGIN
				tmp = drk_mean[cam,mode]-10*stddev(drks_mean[*,cam,mode],/nan)
				IF frm_drk_mean LT tmp THEN BEGIN
					status = status+', low drk <'+strcompress(tmp)+', bad drk'
					frm_drk_mean = badv
				ENDIF
			ENDIF

		ENDIF

	END
	ENDCASE

	CASE finite(frm_drk_mean) OF
	0: BEGIN
		status = status+', nominal='+string(drk_mean[cam,mode],format='(F6.3)')+	$
						' ['+strjoin(string([floor,ceiling]-frm_ped_mean,format='(F6.3)'),',')+']'
		message, /info, status

		nbad[cam,mode] = nbad[cam,mode]+1		; Count # sequential bad frames

		IF NOT bfull[cam,mode] THEN	$
			drk_mean[cam,mode] = smei_property(img, /average,	$
				/dark_current, silent=2)-frm_ped_mean
	END
	1: BEGIN

		; At this point we know that numbers for pedestal and dark current have
		; been obtained. Use them to build statistics from last ncontrol good frames.

		ifrm = (ifull[cam,mode]+step) mod nfull

		peds_mean [ifrm,cam,mode] = frm_ped_mean
		drks_mean [ifrm,cam,mode] = frm_drk_mean

		peds_sigma[ifrm,cam,mode] = frm_ped_sigma
		drks_sigma[ifrm,cam,mode] = frm_drk_sigma

		drks_nout [ifrm,cam,mode] = frm_drk_nout[1]

		ped_mean [cam,mode] = mean(peds_mean [*,cam,mode], /nan)
		drk_mean [cam,mode] = mean(drks_mean [*,cam,mode], /nan)

		ped_sigma[cam,mode] = mean(peds_sigma[*,cam,mode], /nan)
		drk_sigma[cam,mode] = mean(drks_sigma[*,cam,mode], /nan)

		tmp = drks_nout[*,cam,mode]
		drk_nout [cam,mode] = mean(tmp, /nan)
		IF total(finite(tmp)) GT 1 THEN drk_dout[cam,mode] = stddev(tmp, /nan)

		frm_drk_stdev = sqrt( frm_drk_stdev*frm_drk_stdev-frm_ped_stdev*frm_ped_stdev)

		status = status+	$
			','+string(frm_ped_sigma,format='(F6.3)')+strjoin(strcompress(frm_ped_nout))+	$
			','+string(frm_drk_sigma,format='(F6.3)')+strjoin(strcompress(frm_drk_nout))+	$
			string(lr_ratio,format='(F4.1)')

		; We are going to accept pedestal and dark current only if bfull[cam,mode] = 1
		; i.e. if they were submitted to all test and passed.

		bWasFull = bfull[cam,mode]
		ifull[cam,mode] = ifrm
		bfull[cam,mode] = bfull[cam,mode] OR ifrm EQ nfull-1

		CASE bWasFull OF
		0: BEGIN

			status = status+' (rejected)'

			IF silent LE 1 THEN message, /info, status

			frm_ped_mean = badv
			frm_drk_mean = badv

		END
		1: BEGIN

			nbad[cam,mode] = 0

			IF silent LE 1 THEN message, /info, status

			; Subtract pedestal and dark current from frame

			img -= frm_ped_mean+frm_drk_mean

			IF silent EQ 0 THEN BEGIN
				pos = *roi_fov[mode]				; List of pixels defining fov
				tmp = where(img[pos] LT 0, n)

				IF n GT 0 THEN BEGIN
					pos = pos[tmp]
					n = where_common(pos,*roi_bad[mode],absent=tmp)
					n = n_elements(tmp)-(tmp[0] EQ -1)
					IF n GT 0 THEN BEGIN
						pos = pos[tmp]
						message, /info, strcompress(n)+	$
							' pixs zeroed in fov after ped+drk removal'
						pos = ArrayLocation(pos, dim=*roi_siz[mode])
						pos = [pos, img[pos[0,*],pos[1,*]]]
						print, pos[*]
					ENDIF
				ENDIF
			ENDIF

			img = img > 0
		END
		ENDCASE

	END
	ENDCASE

END
ENDCASE

hdr.pedestal     = frm_ped_mean
hdr.dark_current = frm_drk_mean

cam_ok   = finite(frm_ped_mean) AND finite(frm_drk_mean)
cam_full = bfull[cam,mode]

CASE cam_ok OF
0: BEGIN
	hdr.square_residual = BadValue(hdr.square_residual)
	hdr.center_residual = BadValue(hdr.center_residual)
END
1: BEGIN
	hdr.square_residual = smei_property(img, /average, /squares  , silent=2)
	hdr.center_residual = smei_property(img, /average, /centerpix, silent=2)
END
ENDCASE

hdr.ped_sigma = frm_ped_stdev
hdr.drk_sigma = frm_drk_stdev

RETURN  &  END