FUNCTION EDF_COMPUTE, IMAGE

compile_opt idl2

;- Compute the empirical distribution function for a scaled integer image

loc = where((image ge 0) and (image le 32767), count)
if (count eq 0) then message, 'No valid data'
his = histogram(image[loc], min=0, max=32767)
edf = total(his, /cumulative, /double) / count
return, edf

END

;-------------------------------------------------------------------------------

FUNCTION MODIS_EDF_DESTRIPE, IMAGE, REF, MIR, VERSION=VERSION

compile_opt idl2

;- Destripe a MODIS image using a reference detector (20 detector version)
;- IMAGE is an [ncol, nrow] array of scaled integer MODIS image data
;- REF is the index of the reference detector (zero-based)
;- MIR is an [nrow] array of mirror side flags (0 or 1)

;- Set version string
rcs_id = '$Id: modis_edf_destripe.pro,v 1.6 2004/03/19 17:33:29 gumley Exp $'
if (arg_present(version) eq 1) then version = rcs_id

;- Check arguments
if (n_elements(image) eq 0) then message, 'Argument IMAGE is undefined'
if (n_elements(ref) eq 0) then message, 'Argument REF is undefined'
if (n_elements(mir) eq 0) then message, 'Argument MIR is undefined'

;- Get image dimensions and number of earth scans
dims = size(image, /dimensions)
nscan = dims[1] / 10

;- Get index for reference detector on mirror side zero
ref_ind = ref
if (mir[ref] eq 1) then ref_ind = ref + 10

;- Create the EDF for each detector
edf = fltarr(32768, 20)
for det = 0, 19 do begin

  ;- Compute row indices for this detector
  row = lindgen(nscan / 2) * 20 + det
  
  ;- Handle an odd number of scans
  if ((nscan mod 2) eq 1) then begin
    if (det le 9) then begin
      row = [row, max(row) + 20]
    endif else begin
      row = [row, max(row) - 20]
    endelse
  endif

  ;- Compute the EDF for this detector
  edf[*, det] = edf_compute(image[*, row])
  
endfor

;- Create the lookup table that maps each detector to the reference detector
lut = lonarr(32768, 20)
v = findgen(32768)
x = edf[*, ref_ind]
for det = 0, 19 do begin
  if (det ne ref_ind) then begin
    u = edf[*, det]
    linterp, x, v, u, result
    lut[*, det] = round(result)
  endif
endfor

;- Create output image
destripe = image

;- Apply the lookup table to all detectors except the reference detector
for det = 0, 19 do begin
  if (det ne ref_ind) then begin
    row = lindgen(nscan / 2) * 20 + det
    if ((nscan mod 2) eq 1) and (det le 9) then row = [row, max(row) + 20]
    destripe[*, row] = lut[destripe[*, row], det]
  endif
endfor

;- Replace bad destriped values
loc = where((destripe le 0) or (destripe ge 32767), count)
if (count gt 0) then destripe[loc] = image[loc]

;- Replace original invalid values (i.e. missing data)
loc = where(image gt 32767, count)
if (count gt 0) then destripe[loc] = image[loc]

;- Set median value of destriped image to median value of original image
loc_old = where(image le 32767, count_old)
if (count_old gt 0) then begin
  loc_new = where(destripe le 32767, count_new)
  if (count_new gt 0) then begin
    median_old = long(median(image[loc_old]))
    median_new = long(median(destripe[loc_new]))
    median_del = median_new - median_old
    if (median_del ne 0) then $
      destripe[loc_new] = destripe[loc_new] - median_del
  endif
endif

;- Return the destriped image to the caller
return, destripe

END