function MASKMEDFILTLEV, Data, Halfwidth, Lowlev, Highlev, $
	MASK=MASK, DIF=DIF

;+
; NAME:
;	ROUTINE_NAME
;
; PURPOSE:
;	This function generates/adds to a mask of pixels that satisfy
;	a median filter criterion.  A mask is a byte array the same
;	size as Data, with value 1=good pixel and 0=bad pixel.   
;
; CATEGORY:
;	Array processing.
;
; CALLING SEQUENCE:
;
;	Result = MASKMEDFILTLEV(Data, Halfwidth, Lowlev, Highlev)
;
; INPUTS:
;	Data:		The 2D array of values to test.
;	Halfwidth:	Size parameter for median-filter box.  Box is
;				2*Halfwidth+1 elements wide and high.
;	Lowlev:		Lower boundary of good data below median in box.
;
; OPTIONAL INPUTS:
;	Highlev:	Upper boundary of good data above median in
;			box (default: Lowlev).
;
; KEYWORD PARAMETERS:
;	MASK:	The mask to which to add new bad pixels.  All old bad
;		pixels in the mask are preserved.
;
; OUTPUTS:
;	This function returns the (possibly updated) mask.
;
; OPTIONAL OUTPUTS:
;	DIF:	Difference between Data and median-filtered Data.
;
; SIDE EFFECTS:
;	Updates mask.
;
; RESTRICTIONS:
;	The edges are (necessarily) not filtered.
;
; PROCEDURE:
;	For each element, the median is calculated in a box centered
;	on it.  The median filter box is 2*Halfwidth+1 elements wide
;	and high.  Elements whose values are between median-Lowlev and
;	median+Highlev are good.
;
; EXAMPLE:
;	data = fltarr(10, 10)+10.
;	data[4,5] = 100.
;	mask=maskmedfiltlev(Data, 2, 10, 10)
;	print,mask
;	   1   1   1   1   1   1   1   1   1   1
;	   1   1   1   1   1   1   1   1   1   1
;	   1   1   1   1   1   1   1   1   1   1
;	   1   1   1   1   1   1   1   1   1   1
;	   1   1   1   1   1   1   1   1   1   1
;	   1   1   1   1   0   1   1   1   1   1
;	   1   1   1   1   1   1   1   1   1   1
;	   1   1   1   1   1   1   1   1   1   1
;	   1   1   1   1   1   1   1   1   1   1
;	   1   1   1   1   1   1   1   1   1   1
;
; MODIFICATION HISTORY:
; 	Written by:	Joseph Harrington, Cornell. 2001 Apr 25
;			jh@oobleck.astro.cornell.edu
;-

if not keyword_defined(highlev) then highlev = lowlev

if not keyword_defined(mask) then begin
  mask    = byte(data)
  mask[*] = 1
endif

med = median(data, halfwidth*2+1, /even)
dif = data - med

badlist = where(dif lt -lowlev or dif gt highlev)
if isarray(badlist) then begin
  mask[badlist] = 0
endif

return, mask
end