;+
; NAME:
;       BPIXFIX
;
; PURPOSE:
;       Iterative bad pixel replacement funtion
;
; CALLING SEQUENCE:
;       Result = BPIXFIX(input_image, bpix_map, neighbors=neighbors)
;
; INPUTS:
;       input_image = Input image.
;
;       bpix_map = Map of bad pixels (bad pixels = 1, good pixels = 0),
;               with same dimentions as input_image.
;
; OUTPUTS:
;       Result = Input image, with bad pixels replaced by median of
;               surrounding good pixels.
;
; OPTIONAL INPUT KEYWORDS:
;       NEIGHBORS -  minumum number of good neighbors required to replace
;               a pixel in a given iteration. (default = 4)
;
; EXAMPLE:
;       Replace all pixels with value > 1e6 with the median of surroundings.
;               IDL> bpix = BYTE(image)
;               IDL> bpix(*) = 0
;               IDL> bpix(WHERE(image gt 1e6)) = 1
;               IDL> new_image = BPIXFIX(image, bpix)
;
; MODIFICATION HISTORY:
;       Written circa 1999 by A. Bouchez (Caltech).
;       Optimized, commented, changed to function.  March 2004. A. Bouchez.
;       Set default neighbors to 3 for corner bad pix. M. van Dam, Feb 2005.
;-
FUNCTION        BPIXFIX,im0,bp0,neighbors=neighbors

  if not KEYWORD_SET(neighbors) then $  ; number of good neighbors required
    neighbors = 3                       ;  to fix pixel
  sz = (SIZE(im0))[1:2]                 ; get image dimentions
  bp = bp0                              ; temporary bad pixel map
  im = im0                              ; corrected image
  nb = TOTAL(bp0)                       ; number of bad pixels remaining

;In each iteration, correct only those bad pix with sufficient good 'neighbors'
  while nb gt 0 do begin
    wb = WHERE(bp,nb)                   ; find bad pixels, make column vector
    xb = wb mod sz[0]                   ; x coord of bad pixels
    yb = wb / sz[0]                     ; y coord of bad pixels
    gp = 1 - bp                         ; temporary good pixel map
    for n=0,nb-1 do begin
      sgp = gp[(xb[n]-1)>0:(xb[n]+1)<(sz[0]-1), $
               (yb[n]-1)>0:(yb[n]+1)<(sz[1]-1)]
      if TOTAL(sgp) ge neighbors then begin
        sim = im[(xb[n]-1)>0:(xb[n]+1)<(sz[0]-1), $
                 (yb[n]-1)>0:(yb[n]+1)<(sz[1]-1)]
        im[xb[n],yb[n]] = MEDIAN(sim[WHERE(sgp)])
        bp[xb[n],yb[n]] = 0
      endif
    endfor
  endwhile

RETURN,im
END