FUNCTION MEDSMOOTH,ARRAY,WINDOW
;+
; NAME:
;       MEDSMOOTH
;
; PURPOSE:
;       Median smoothing of a vector, including points near its ends.
;
; CALLING SEQUENCE:
;       SMOOTHED = MEDSMOOTH( VECTOR, WINDOW_WIDTH )
;
; INPUTS:
;       VECTOR  = The (1-d numeric) vector to be smoothed
;       WINDOW = Odd integer giving the full width of the window over which 
;               the median is determined for each point.     (If WINDOW is
;               specified as an even number, then the effect is the same as
;               using WINDOW+1)   
;
; OUTPUT:
;       Function returns the smoothed vector
;
; PROCEDURE:
;       Each point is replaced by the median of the nearest WINDOW of points.
;       The width of the window shrinks towards the ends of the vector, so that
;       only the first and last points are not filtered. These points are 
;       replaced by forecasting from smoothed interior points.
;
; EXAMPLE:
;       Create a vector with isolated high points near its ends
;       IDL> a = randomn(seed,40) & a[1] = 10  & a[38] = 10
;       Now do median smoothing with a 7 point window 
;       IDL> b = medsmooth(a,7)
;       Note that, unlike MEDIAN(), that MEDSMOOTH will remove the isolated
;       high points near the ends.
; REVISION HISTORY:
;       Written, H. Freudenreich, STX, 12/89
;       H.Freudenreich, 8/90: took care of end-points by shrinking window.
;       Speed up using vector median when possible  W. Landsman February 2002
;-

 LEND = N_ELEMENTS(ARRAY)-1
 IF (LEND+1) LT WINDOW THEN BEGIN
   message,/CON, $
         'ERROR - Size of smoothing window must be smaller than array size'
   RETURN,ARRAY
 ENDIF

 OFFSET = FIX(WINDOW/2)

 smoothed = median(array, window )

; Fix the ends:
 NUMLOOP = (WINDOW-1)/2 - 1
 IF NUMLOOP GT 0 THEN BEGIN
   FOR J=1,NUMLOOP DO BEGIN 

     LEN = 2*J+1
     SMOOTHED[J] = MEDIAN(ARRAY[0:LEN-1])
     SMOOTHED[LEND-J] =  MEDIAN(ARRAY[LEND-LEN+1:LEND]) 

   ENDFOR
ENDIF

; Now replace the very last and first points:
 Y0 = 3.*ARRAY[0]-2.*ARRAY[1]         ; Predicted value of point -1
 SMOOTHED[0] = MEDIAN([Y0,ARRAY[0],ARRAY[1]])
 Y0 = 3.*ARRAY[LEND]-2.*ARRAY[LEND-1] ; Predicted value of point LEND+1
 SMOOTHED[LEND] = MEDIAN([Y0,ARRAY[LEND],ARRAY[LEND-1]])
               
 RETURN,SMOOTHED
 END