; NAME:
;       pro interp.pro
;
; PURPOSE:
;       This program identifies statistical outliers and replaces them
;       with interpolated values. It is not very sophisticated.
;
; CALLING SEQUENCE:
;       interp, inFile, outFile [, sigrej= [, /help]]
;
; INPUTS:
;       inFile  - A 2-dimensional FITS image
;       
; OUTPUTS
;       outFile - A 2-dimensional FITS image
;
; KEYWORD PARAMETERS:
;       help             - Print a help message
;       sigrej           - Sigma clipping threshold in sigma
;
; EXAMPLE
;       Fill in 2.5-sigma outliers with interpolated values
;       IDL> interp, 'theFile.fits', 'theFixedFile.fits', sigrej=2.5
;
;
; REFERENCE:
;
;
; MODIFICATION HISTORY:
;       Written by:
;       B.J. Rauscher, 28 April 2003
;
;-

pro interp, inFile, outFile, sigrej=sigrej, help=help

       ;; Print help message if requested
       if (keyword_set(help)) then begin
           stop, 'Usage: interp, inFile, outFile [, sigrej= [, /help]]'
       endif

       ;; Set defaults
       if (not keyword_set(sigrej)) then sigrej = 3.0

       pix_counter = 0L

       ;; Read in the data
       fits_read,inFile,data,header

       for loopCounter = 0L, 4 do begin

           print, 'Starting iteration ', loopCounter

           ;; Median filter
           medData = data - median(data, 5, dimension=2)

           ;; Calculate noise in the input image
           djs_iterstat, medData, sigrej=sigrej, sigma=sigma, mask=mask
           print, 'Sigma = ', sigma

           xsize = n_elements(data[*,0])
           ysize = n_elements(data[0,*])
           for y = 0L, ysize - 1 do begin
               for x = 0L, xsize - 1 do begin

                   if (mask[x,y] eq 0) then begin

                       pix_counter = pix_counter + 1L

                       sum = 0.0
                       count = 0L
                       for _y = long(y - 1), long(y + 1) do begin
                           for _x = long(x - 1), long(x + 1) do begin
                               if (_x eq x and _y eq y) then continue
                               if (_x eq xsize or _y eq ysize) then continue
                               if (_y eq -1 or _x eq -1) then continue
                               sum = sum + data[_x,_y]
                               count = count + 1L
                           endfor
                       endfor
                       data[x,y] = sum / count

                   endif
               endfor
           endfor

       endfor
       
       print, 100. * float(pix_counter) / xsize / ysize, ' Percent of pixels replaced'
       fits_write, outFile, data

end