function gentle_filter, image, second_max = second_max

;+
;
; NAME: gentle_filter
;
; PURPOSE:
;      This function filters images to remove outliers somehwat less
;      aggressively than more traditional filtering (like, median or
;      local average smoothing). The filter restricts the value of
;      every pixel to lie between the min and max values of its eight neighbors.
;      This function is useful for filtering images with few, strong outlier
;      pixels while preserving the overall sharpness of the image.
;
;      This function is especially useful for images produced by the
;      PROBA2/SWAP routines p2sw_long_movie and swap_mosaic.
;
; INPUTS:
;      IMAGE:   The two-dimensional image array to be filtered
;
; RETURNS:
;      OUTPUT_IMAGE: The resulting filtered image
;
; KEYWORD PARAMETERS:
;      SECOND_MAX:   Causes the program to filter on the second
;                    brightest and second faintest pixel in the neighborhood.
;
; MODIFICATION HISTORY:
;      Written by:   D. B. Seaton, ROB, 2-Jun-2015.
;       Jul 2015     Improved efficiency of second_max option. (dbs)
;
;-
  
  ;; Use second max if keyword set
  second_max = keyword_set(second_max)

  ;; Get properties of image to be filtered
  im_size = size(image)

  ;; Generate an array that can quickly return the neighborhood max/min for
  ;; every pixel in the image. Note that the original pixel is
  ;; excluded from the neighborhood.
  test_im = [[[shift(image, 1, 0)]], [[shift(image, 1, 1)]], [[shift(image, 0, 1)]], $
             [[shift(image, -1, 1)]], [[shift(image, -1, 0)]], [[shift(image, -1, -1)]], $
             [[shift(image, 0, -1)]], [[shift(image, 1, -1)]]]

  ;; Find the neighborhood max for each pixel
  if ~(second_max) then $
     max_im = max(test_im, min = min_im, dim = 3) $
  else begin $
     ;; If second_max is set, find primary max/min and eliminate these
     ;; values by setting them to NaNs.  Search again to find
     ;; secondary max and min values in the neighborhood.
     max_im = max(test_im, min = min_im, dim = 3, sub_max, subscript_min = sub_min)
     test_im(sub_max) = !values.f_nan
     test_im(sub_min) = !values.f_nan
     max_im = max(test_im, min = min_im, dim = 3, /nan)
  endelse

  ;; Generate a new image in which no pixel is brighter or dimmer than
  ;; the brightest or dimmest surrounding pixel
  out_image = (image > min_im) < max_im
  
  return, out_image
	
end