FUNCTION DIcal_decompress, in, fitsHdr, METHOD=method, FN=fn, FLAGS=flags, NOISE=noise

;;------------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Restores and image that has been compressed. This is done through a 
;;	LUT identified in fitsHdr. To determine which LUT to use, a SQL
;;	database is querried to find the file that contains the LUT. The
;;	file is a FITS file that contains a 256x2 array. The first column
;;	is the lowest 14-bit value for the bin while the second column is
;;	the highest 14-bit value for the bin.
;;
;; CALLING SEQUENCE:
;;	out = DIcal_decompress(in, fitsHdr, FN=fn)
;;
;; REQUIRED INPUTS:
;;	in - The image to which this pipeline element is going to be applied
;;	fitsHdr - The FITS header for the image
;;
;; OUTPUTS:
;;	RETURN - the image after going through this calibration step
;;
;; OPTIONAL INPUT KEYWORDS:
;;	METHOD - The method used to distribute the data inside an 
;;		uncompressed bin. Options are:
;;		1 - Average (default)
;;		2 - Normal distribution
;;		3 - Uniform distribution
;;	FN - Filename of the decompression LUT. If not specified, the filename
;;		is found using a database
;;	FLAGS - Flag map that will be modified
;;
;; OPTIONAL OUTPUT KEYWORDS:
;;	NOISE - The quantization noise for each pixel squared
;;
;; EXAMPLE:
;;      IDL> imgOut = DIcal_decompress(imgIn, fitsHdr)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	getSQLServer - Info on how to get to the database
;;	DI_sqlQuery - Sends queries to the SQL server
;;	getLocFn - Gets a file and stores it locally
;;	SXPAR - Read the FITS header
;;	FXADDPAR - Change values of header items
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;   2005-06-07  M. Desnoyer    Added the ultra compression flag and fixed LUT numbering
;;
;;------------------------------------------------------------------------------

;; Do error handling 
;CATCH, error
error=0

IF error NE 0 THEN BEGIN
	CATCH, /CANCEL
	message, 'Decompress - ' + !ERROR_STATE.MSG, /NONAME
ENDIF

@dical_flags

;; Set the distribution method
IF NOT keyword_set(method) THEN meth = 1 ELSE meth = method

;; Get the date stamp on the image
date = SXPAR(fitsHdr,'OBSDATE', COUNT=c1)

;; Get the LUT used
lutNum = SXPAR(fitsHdr,'IMGH001', COUNT=c2)

;; Get the instrument
inst = SXPAR(fitsHdr,'INSTRUME', COUNT=c3)

;; Make sure the header contained all the parameters
IF (c1 EQ 0) OR (c2 EQ 0) OR (c3 EQ 0) THEN $
	message, 'Invalid FITS header', /NONAME

;; Check the flag maps
dim = size(in,/dimensions)
IF NOT keyword_set(flags) THEN flags = bytarr(dim[0],dim[1])

;; Get the filename of where the LUT is stored
IF NOT keyword_set(fn) THEN BEGIN

	;; The LUT filenames will be stored in the SQL database so they can be 
	;; searched easily by time
	serv = getSQLServer()
	db = 'di_calib' 

	;; Generate the SQL query to get the filename containing the LUT
	tblNm = 'DECOMPRESS' 
	sel = ['Filepath'] 
	cond = "(DATE <= '"+date+"') and (INSTRUMENT='"+strtrim(inst,2) $
		+"') and (Idx="+strtrim(lutNum,2)+ $
		") order by DATE DESC, Version desc limit 1"

	;; Query the database and make sure it returns useful data
	webfn = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)
	IF min(dim) EQ 0 THEN $;; The database didn't have anything
		message, 'Could not find valid decompression LUT', /NONAME

	;; Get the local filename of the LUT to use
	fn1 = getLocFn(webfn, subdir='DECOMPRESS')
ENDIF ELSE fn1 = fn

;; Get the LUT. 
lut = readFITS(fn1[0])

;; Find any pixels that are ultra compressed. As of now, this means IR data in LUT 2 and 3
;; between values 250-255
IF strtrim(inst,2) EQ 'HRIIR' AND lutNum GE 2 THEN BEGIN
	w = where(in GE 250 AND in LE 255,cnt)
	IF cnt GT 0 THEN flags[w] = flags[w] OR FLAG_ULTRACOMPRESS
ENDIF

;; Perform the decompression
out = 0>in<255
lowout = lut[out,0]
hiout = lut[out,1]

;; Calculate the quantization noise squared, which is the bin size squared divided by 12
sz = size(out,/dimensions)
noise = (hiout-lowout+1.0)^2/12
noise = reform(noise,sz[0],sz[1])

;; Choose the method of distribution
CASE meth OF
	;; Average of the bin
	1: dist = 0.5
	
	;; Normal distribution centered at the middle of the bin
	2: dist = 0>(randomn(systime(1),n_elements(out))/4+0.5)<1

	;; Uniform distribution
	3: dist = randomu(systime(1),n_elements(out))

	ELSE: message, 'Invalid distribution method', /noname
ENDCASE

;; Generate the output with the chosen distribution
out = round(lowout + dist*(hiout-lowout+1)-0.5)
out = reform(out,sz[0],sz[1])

;; If it's IR, also have to invert the signal
IF strtrim(inst,2) EQ 'HRIIR' THEN BEGIN
	out = 16383-out
	w_wrap = where(out GT 16383,wrapCnt)
	IF wrapCnt GT 0 THEN out[w_wrap] = out[w_wrap]-16384
ENDIF

;; Update header to say image is no longer compressed
fxaddpar, fitsHdr, 'COMPRESS', 'DECOMPRESSED'
FXADDPAR, fitsHdr, 'CMPRESSN', 'T'
fxaddpar, fitsHdr, 'LUTTABLE', $
	strupcase(strmid(fn1[0], strpos(fn1[0], path_sep(), /reverse_search)+1))
fxaddpar, fitsHdr, 'CMPRMETH', fix(meth)

RETURN, out

END