FUNCTION POINT_REMOVER,MAP,WIDTH,SIGMA_MULT,ITMAX, BOTH=WHATEVER, FLOOR=MINV
;+
; NAME:
;	POINT_REMOVER
;
; PURPOSE:
;	Remove isolated HIGH SIGNAL/NOISE elements from a 2D array, replacing
;	them with a local median. This routine iterates until convergence or
;	a user-supplied maximum number. This is NOT a low-pass filter. Only
;	points passing the S/N cut are affected. 
;
; CALLING SEQUENCE:
;	NewMap = POINT_REMOVER( Map, Width, Sigma_Mult, ItMax,[ /BOTH, FLOOR= ])
;
; INPUT ARGUMENT:
;	MAP = the array to de-point
;	WIDTH = the width of the square moving neighborhood centered on the 
;		pixel in question. If half the pixels within this neighborhood
;		 are vacant, nothing is done.
;
; OPTIONAL INPUT ARGUMENT:
;	SIGMA_MULT = the number of standard deviations by which a pixel must 
;		exceed the local background for it to be replaced. DEFAULT = 2.0
;	ITMAX = the maximum number of iterations. DEFAULT = 20
;
; OPTIONAL INPUT KEYWORDS:
;	BOTH - if set, the absolute value of the difference between a pixel and 
;		its neighborhood is used, so that low pixels may also be 
;		replaced. Else, only + high points will be replaced.
;	FLOOR = the value representing an empty pixel. Default = 0.0
;
; OUTPUT:
; 	POINT_REMOVER returns the array, with point sources removed.
;
; REVISION HISTORY:
;	Written, H.T. Freudenreich, HSTX, ?/90 or 91
;	HF, 3/94, to include FLOOR and USE_ABS keywords
;-

; Fill in default parameters:
IF N_PARAMS() LT 4 THEN ITMAX=20
IF N_PARAMS() LT 3 THEN SIGMA_MULT=2.0
IF N_PARAMS() LT 2 THEN WIDTH=7

IF KEYWORD_SET(WHATEVER) THEN USE_ABS=1 ELSE USE_ABS=0
IF KEYWORD_SET(MINV)     THEN EMPTY=MINV ELSE EMPTY=0.0

; Set some needed constants:
SYZ=SIZE(MAP)
NX=SYZ(1)
NY=SYZ(2)
OFF=FIX(WIDTH)/2
MINPIX = WIDTH*WIDTH/2 ; the minimum size of the neighborhood

AMAP    = MAP
ITNUM   = 0

again:
ITNUM = ITNUM + 1
NUM_REJ = long(0)
TEMP = AMAP
FOR I=OFF,NX-OFF-1 DO BEGIN
  FOR J=OFF,NY-OFF-1 DO BEGIN
;   Extract the neighborhood:
    A=TEMP(I-OFF:I+OFF,J-OFF:J+OFF)
;   Are there enough non-zero points?
    Q=WHERE(A GT EMPTY,N) 
    IF N GT MINPIX THEN BEGIN
       A=A(Q)
;      Sort the data, get the median and interquartile range:
       A=A(SORT(A))
       N1=N/4
       N3=3*N1
       N2=N/2
;      The median:
       IF N MOD 2 NE 0 THEN BCKGD=A(N2) ELSE BCKGD=.5*(A(N2-1)+A(N2))
;      The interquartile range:
       SIGMA=(A(N3)-A(N1))/1.35
       IF SIGMA LT EMPTY THEN SIGMA=TOTAL(ABS(A-BCKGD))/.8/N
       HITE= TEMP(I,J)-BCKGD
       IF USE_ABS EQ 1 THEN HITE=ABS(HITE)
       IF HITE GT (SIGMA_MULT*SIGMA) THEN BEGIN
          AMAP(I,J) = BCKGD
          NUM_REJ = NUM_REJ + long(1)
       ENDIF
    ENDIF
  ENDFOR
ENDFOR
PRINT,'Iteration number ',ITNUM,'. Rejected ',NUM_REJ,' pixels.'
IF (NUM_REJ GT 1) AND (ITNUM LT ITMAX) THEN GOTO,AGAIN
       
RETURN,AMAP
END