;+
; NAME:
;    FIXIMAGES
; PURPOSE:
;    Remove cosmic rays from a set of images, using the median and
;    then summing them in combining
; CALLING SEQUENCE:
;    fixImages, input[, outfile=outfile][,sigma=sigma][,keep=keep][,use=use]
; EXAMPLES:
;    fixImages, 'science1comb'
;        combines the images listed in the textfile 'science1comb',
;        which contains the filenames without the '.fits'
;    fixImages, ['13_3Da0008_1', '13_3Da0007_1', '13_3Da0006_1', '13_3Da0005_1']
;        combines the images named in the string array
; KEYWORDS:
;    input = string containing the filename of a textfile with the
;    filenames of the images that have to be combined without the
;    .fits, or a string-array of filenames without the .fits of images
;    to combine
;    outfile = the name of the output file without the .fits. If not
;    given, then the output file is combCorr.fits
;    sigma = the sigma level at which the cosmic rays are clipped. If
;    not set, then 10 sigma is applied
;    keep = if set, then the images without the cosmic rays are
;    kept. They will be named with corr before the .fits
;    use
;    use = an array of the size of the number of input images,
;    with 1's for inclusion of the image in the summed image and 0's
;    for the images only to be used in construction of the median
;    image.
; HISTORY:
;    Before 08/02/2006 written and created by Eduard Westra. If you
;    used this script, please contact me and acknowledge me in your
;    paper.
;    08/02/2006 EW added the use keyword
;    10/07/2007 EW In the case that no EXPTIME is found, assume that
;    it has the ESO header starting with HIERARCH and ending with
;    EXPTIME and parse it to obtain the values.
;    10/07/2007 EW Made sure that the header written to the output
;    file is that of the first file read in.
;-

FUNCTION myStdDev, x, mean=mean
xx = x
s = size(xx)
;mean = TOTAL(x, 1)/s[1]
max = REFORM(xx[0,*,*])
FOR i=1, s[1]-1 DO BEGIN
    max = max > xx[i,*,*]
ENDFOR
sum = TOTAL(xx, 1)
clippedMean = (sum-max)/(s[1]-1)

;xdiff = MAKE_ARRAY(DIMENSION=s[1:s[0]], TYPE=s[s[0]+1], /nozero)
FOR i=0, s[1]-1 DO BEGIN
    xx[i,*,*] = xx[i,*,*] - max
END
mean = clippedMean
clippedNegMean = TOTAL(xx, 1)/(s[1]-1)
FOR i=0, s[1]-1 DO BEGIN
    xx[i,*,*] = xx[i,*,*] - clippedNegMean
END
xx = xx*xx
xx = TOTAL(xx,1)
xx = xx - clippedNegMean^2
xx = xx / (s[1] - 2)
xx = SQRT(xx)
;print, "Subtract clippedNegMea: ", xx[*,695,9]
;print, "Square it             : ", xx[*,695,9]
;print, "Sum it                : ", xx[695,9]
;print, "clippedNegMean^2      : ", (clippedNegMean^2)[695,9]
;print, "Subtract clipNedMean^2: ", xx[695,9]
;print, "Divide by Ndim - 2    : ", xx[695,9]
;print, "Take the sqrt         : ", xx[695,9]
;print, ""
;print, "Max                   : ", max[695,9]
;print, "Clipped mean          : ", clippedMean[695,9]
;print, "Unclipped stddev      : ", stddev(sub)
;print, "Clipped stddev (IDL)  : ", stddev((sub[reverse(sort(sub))])[1:*])
RETURN, xx
END

FUNCTION sigClip, x, nhigh=nhigh
;IF NOT KEYWORD_SET(nhigh) THEN nhigh = 0
RETURN, STDDEV((x[reverse(sort(x))])[nhigh:*])
END

FUNCTION makeStdDev, x
;Calculates the standard deviation through the first axis of an array
;and returns an array with the standard deviation of the dimensions
;remaining.
s = SIZE(x)
nDim = s[0]
IF ((nDim LT 1) AND (nDim GE 3)) THEN MESSAGE, 'Not enough dimensions!'
result = MAKE_ARRAY(s[2:nDim],TYPE=s[nDim+1])
FOR i=0,s[2]-1 DO BEGIN
    FOR j=0,s[3]-1 DO BEGIN
        result[i,j] = STDDEV(x[*,i,j])
    ENDFOR
ENDFOR
RETURN, result
END

;PRO fixImages, inFiles, medFile, sigma=sigma, outfile=outfile, keep=keep
PRO fiximages, inFiles, sigma=sigma, outfile=outfile, keep=keep, help=help, use=use
IF KEYWORD_SET(help) THEN BEGIN & doc_library,'fiximages' & RETURN & ENDIF
IF NOT KEYWORD_SET(inFiles) THEN BEGIN & doc_library,'fiximages' & RETURN & ENDIF
IF NOT KEYWORD_SET(sigma) THEN sigma = 10.
IF NOT KEYWORD_SET(outfile) THEN outfile = 'combCorr'
IF NOT KEYWORD_SET(keep) THEN keep = 0

IF N_ELEMENTS(inFiles) EQ 1 THEN BEGIN
    list = inFiles
    OPENR, lun, list, /GET_LUN
    inFiles = STRARR(FILE_LINES(inFiles))
    READF, lun, inFiles
    CLOSE, lun
    FREE_LUN, lun
ENDIF
IF NOT KEYWORD_SET(use) THEN use = uintarr(N_ELEMENTS(inFiles)) + 1

header  = HEADFITS(inFiles[0]+'.fits')
baseExp = FXPAR(header, 'EXPTIME')
IF baseExp LE 0 THEN BEGIN
;Assume ESO headers
    print, "Could not find EXPTIME. Assuming ESO headers..."
    baseExpTmp = header[WHERE(STRMATCH(header, '*EXPTIME*') EQ 1)]
    baseExp    = DOUBLE((STREGEX(baseExpTmp, '=(.*)/', /EXTRACT, /SUB))[1])
ENDIF
xsize  = FXPAR(header, 'NAXIS1')
ysize  = FXPAR(header, 'NAXIS2')

image = DBLARR(N_ELEMENTS(inFiles), xsize, ysize)
scale = DBLARR(N_ELEMENTS(inFiles))

FOR i=0,N_ELEMENTS(inFiles)-1 DO BEGIN
    hdr      = HEADFITS(inFiles[i]+'.fits')
    exp      = FXPAR(hdr, 'EXPTIME')
    IF exp LE 0 THEN BEGIN
    ;Assume ESO headers
        print, "Could not find EXPTIME. Assuming ESO headers..."
        expTmp = hdr[WHERE(STRMATCH(hdr, '*EXPTIME*') EQ 1)]
        exp    = DOUBLE((STREGEX(expTmp, '=(.*)/', /EXTRACT, /SUB))[1])
    ENDIF
    ;scale[i] = medExp/exp
    scale[i] = baseExp/exp
    print, "Reading in ("+(use[i] EQ 1 ? "Y" : "N")+") "+inFiles[i]+'.fits... Exptime: '+STRTRIM(STRING(exp),2)+' Scale '+STRTRIM(STRING(scale[i]), 2)
    ;image[i,*,*] = scale[i]*MRDFITS(inFiles[i]+'.fits', /SILENT, ROWS=rows)
    image[i,*,*] = scale[i]*MRDFITS(inFiles[i]+'.fits', /SILENT)
ENDFOR
;stdDev = CMAPPLY('USER:sigClip', image, 1, FUNCTARGS={nhigh : 1})
stdDev = myStdDev(image, mean=median)

newImage = DBLARR(N_ELEMENTS(inFiles), xsize, ysize)
FOR i=0,N_ELEMENTS(inFiles)-1 DO BEGIN
    out = REFORM(image[i,*,*])
    sig = (out - median)/stdDev
    idx = WHERE(sig GT sigma, count)
    IF (count NE 0) THEN BEGIN
        out[idx] = median[idx]
    ENDIF
    IF use[i] EQ 1 THEN BEGIN
        ;Start writing the corrected files
        ;First read the header
        hdr = ''
        hdr = HEADFITS(inFiles[i]+'.fits')
        ;Write the output
        IF (keep NE 0) THEN BEGIN
            print, "Writing "+inFiles[i]+'corr.fits...'
            MWRFITS, FLOAT(out/scale[i]), inFiles[i]+'corr.fits', hdr
            ;MWRFITS, sig, inFiles[i]+'sigcorr.fits', hdr
        ENDIF
        newImage[i,*,*] = out/scale[i]
    ENDIF
ENDFOR
;totExp = TOTAL(1./scale)*medExp
idx = WHERE(use EQ 1)
totExp = TOTAL(1./scale[idx])*baseExp
FXADDPAR, header, 'EXPTIME', totExp
FXADDPAR, header, 'HISTORY', 'Created with fiximages.pro on '+systime()
FXADDPAR, header, 'HISTORY', 'Created with the following input files:'
FOR i=0, N_ELEMENTS(inFiles)-1 DO BEGIN
    IF use[i] EQ 1 THEN BEGIN
        FXADDPAR, header, 'HISTORY', inFiles[i]
    ENDIF
ENDFOR
FXADDPAR, header, 'HISTORY', 'giving a total exposure time of: '+STRTRIM(STRING(totExp),2)+' sec'
print, "Writing "+outfile+".fits with total exposure time of "+STRTRIM(STRING(totExp),2)+" sec"
;MWRFITS, FLOAT(CMAPPLY('+', newImage, 1)), outfile+'.fits', header
MWRFITS, FLOAT(TOTAL(newImage, 1)), outfile+'.fits', header
END