;+
; NAME:
;	MEDIM
;
; PURPOSE:
;	Compute the pixel-by-pixel median of a set of FITS images.
;
; EXPLANATION:
;	Reads in a set of FITS images, and computes their pixel-by-pixel
;	median.  Subarrays may be used to reduce system memory requirements
;	(e.g. use 2x2 subarrays for large sets of 1024^2 pixel images.
;
; CALLING SEQUENCE:
;	Result = MEDIM(Filenames, Headers, subim=, reject=, exten_no=, /pharo)
;
; INPUTS:
;       Filenames = String vector containing the names of the FITS files
;               to be read.   If the filename has a ".gz" extension, it
;               will be treated as a gzip compressed file.   If it has
;               a ".Z" extension, it will be treated as a Unix compressed file.
;
; OUTPUTS:
;       Result = Pixel-by-pixel median image.
;
; OPTIONAL OUTPUT:
;       Headers = String array containing the headers from the FITS file.
;
; OPTIONAL INPUT KEYWORDS:
;       SUBIM = 2-element vector defining the subframe factor to be used
;		For example, setting subfr=[2,2] will read files 4 times,
;		keeping only one fourth of each image in memory each pass.
;
;	REJECT = Fraction of highest pixels to reject from median.  Default=0.
;
;       EXTEN_NO - non-negative scalar integer specifying the FITS extension to
;               read.  For example, specify EXTEN = 1 or /EXTEN to read the
;               first FITS extension.
;
;	/PHARO - Keyword passed to RFITS, which calls READPHARO instead of
;		READFITS to correctly interpret 512x512x4 PHARO image cubes.
;
;       /SILENT - Normally, will display subframe progress to the terminal.
;                 The SILENT keyword will suppress this.
;
; PROCEDURES USED:
;	Functions: READFITS(), READPHARO()
;
; MODIFICATION HISTORY:
;	Create 10/03 by A. Bouchez, W.M. Keck Observatory
;	Modified to use RFITS.pro; 4/04 A. Bouchez, W.M. Keck Observatory
;	Added PHARO keyword; 5/07 A. Bouchez, Caltech Optical Obs.
;	Corrected "reject" keyword behaviour; 5/24/07 A. Bouchez.
;-
FUNCTION	MEDIM, filenames, headers, subim=subim, reject=reject, $
  exten_no=exten_no, pharo=pharo, silent=silent

  hd0 = HEADFITS(filenames[0], /silent)
  sz = [SXPAR(hd0,'NAXIS1'), SXPAR(hd0,'NAXIS2'), N_ELEMENTS(filenames)]
  if KEYWORD_SET(pharo) then $
    sz[0:1] = sz[0:1] * 2

  if not KEYWORD_SET(silent) then $
    MESSAGE, /info, 'Computing median of ' + $
      STRING(sz[[2,0,1]],f='(I0.0,X,I0.0,"x",I0.0)') + ' frames.'

  if not KEYWORD_SET(subim) then $			; set default subarray
    if sz[2] gt 20 then subim = [2,2] else subim = [1,1]
  ssz = sz / subim

  itype = SXPAR(hd0,'BITPIX')
  case itype of
      8: med = BYTARR(sz[0],sz[1],/nozero)
     16: med = INTARR(sz[0],sz[1],/nozero)
     32: med = LONARR(sz[0],sz[1],/nozero)
     64: med = LON64ARR(sz[0],sz[1],/nozero)
    -32: med = FLTARR(sz[0],sz[1],/nozero)
    -64: med = DBLARR(sz[0],sz[1],/nozero)
  endcase

  if KEYWORD_set(reject) then $
    MESSAGE, /info, 'Excluding ' + STRTRIM(ROUND(sz[2]*reject),2) + $
      ' frames from median.'
  nq = subim[0] * subim[1]
  for q=0,nq-1 do begin
    if (MAX(subim) gt 1) and (not KEYWORD_SET(silent)) then $
      MESSAGE, /info, '  Computing subframe ' + $
        STRING([q+1,nq],f='(I0.0," of ",I0.0)')

    c = [q mod subim[0], q / subim[1]] * ssz		; lower-left pixel
    subfr = [c[0]+[0,ssz[0]-1],c[1]+[0,ssz[1]-1]]
    im = RFITS(filenames, headers, subfr=subfr, $
           exten_no=exten_no, /silent, pharo=pharo)
    for i=0,ssz[0]-1 do begin				; take pix-by-pix
      for j=0,ssz[1]-1 do begin				;  median of subarray
        if KEYWORD_SET(reject) then begin
          pix = im[i,j,*]
          spix = SORT(pix)
          med[c[0]+i,c[1]+j] = MEDIAN(pix[spix[0:ROUND(sz[2]*(1-reject))-1]])
        endif else $
          med[c[0]+i,c[1]+j] = MEDIAN(im[i,j,*])
      endfor
    endfor
  endfor
 
RETURN,med
END