pro removeblemish,im,imout,lowlim=lowlim,uplim=uplim,width=width,sigma_mult=sigma_mult,blemishin=blemishin,noplot=noplot,$
                  blemishout=blemishout

;+
;
; NAME:
;   removeblemish
; PURPOSE:
;   Remove blemishes from DISR images
; CATEGORY:
;    data reduction
; CALLING SEQUENCE:
;    removeblemish,im,imout,lowlim,uplim,width=width,sigma_mult=sigma_mult,blemishin=blemishin,noplot=noplot,$
;                  blemishout=blemishout
; INPUTS:
;    im:  Input image
; OPTIONAL INPUT PARAMETERS:
;   NONE
; OUTPUTS:
;    imout:  Output image with blemishes removed.
; OPTIONAL OUTPUT PARAMETERS:
;   NONE
; KEYWORDS:
;   lowlim, uplim:  If the result is plotted, the image can be stretched by specifying the lowest and highest value
;                   displayed. Default is no stretching.
;           width:  Width of the region for the median filter
;      sigma_mult:  Limit for the comparison between a point and its environment. IF ABS(pixel value-median) >
;                   standard deviation*sigma_mult, the point will be replaced.
;       blemishin:  Array of size (2,number of blemishes) which may specify the points where the blemishes are expected.
;                   If such a list is provided, the points listed are replaced without doing any additional filtering.
;          noplot:  If set, the plotting of the results is disabled.
;      blemishout:  If set, the pixels which were replaced are listed in the same format as for blemishin.
; COMMON BLOCKS:
;   NONE
; SIDE EFFECTS:
;   NONE
; RESTRICTIONS:
;   NONE
; PROCEDURE:
;   Either replace the points listed in blemishin by their local environment (the size of which is psecified by width)
;   or apply a median filter to the images using point_remover.
;
; MODIFICATION HISTORY:
;   17 01 05, MK
;-


  ; set defaults if keywords were not provided
if n_elements(width) eq 0 then width = 5
if n_elements(sigma_mult) eq 0 then sigma_mult = 4.


  ; If blemishlist is given, just replace the pixels listed and return.
if n_elements(blemishin) gt 0 then begin
  imout = im
  sz = size(blemishin)
  for i = 0,sz[2]-1 do begin
     imout[blemishin[0,i],blemishin[1,i]] = $
        median(im[blemishin[0,i]-width/2:blemishin[0,i]+width/2,blemishin[1,i]-width/2:blemishin[1,i]+width/2])
  endfor
  return
endif

  ; Otherwise, run median filter
imout = point_remover(im,width,sigma_mult,/both)


  ; Find out which pixels were replaced
ind = twodindex(imout,where(imout ne im))

  ; Plot both images and mark the replaced pixels unless keyword noplot is set.
sz = size(ind)
if not keyword_Set(noplot) then begin
  loadct,0
  wset,0
  erase
  if n_elements(lowlim) eq 0 then lowlim = min(im)
  if n_elements(uplim) eq 0 then uplim = max(im)
  tvscl,im>lowlim<uplim
  wset,1
  erase

  tvscl,imout>lowlim<uplim
  loadct,3
  repeat begin
     a = get_kbrd(1)
  endrep until a eq 'n'
  for i = 0,sz[2]-1 do plots,ind[0,i],ind[1,i],psym=3,color=192,/device
endif
  ; return list of replaced pixels through keyword.
blemishout=ind
end