function remove_tomo_artifacts, image, width=width, threshold=threshold, $
zingers=zingers, rings=rings, diffraction=diffraction

;+
; NAME:
;    REMOVE_TOMO_ARTIFACTS
;
; PURPOSE:
;    Removes artifacts from tomography images and sinograms.
;
; CATEGORY:
;    Tomography data processing
;
; CALLING SEQUENCE:
;    Result = REMOVE_TOMO_ARTIFACTS(Image, Width=Width, Threshold=threshold, /Zingers,
;                                    /Rings, /Diffraction)
;
; INPUTS:
;    Image:    The input array from which artifacts are to be removed.
;
; KEYWORD PARAMETERS:
;    /ZINGERS:        Set this keyword to remove zingers from the input image
;    /RINGS:            Set this keyword to remove ring artifacts from a sinogram
;    /DIFFRACTION:     Set this keyword to removed diffraction artifacts from a sinogram
;    WIDTH:            Set this keyword to adjust the size of the filter kernal used in the
;                    artifact removal. The default is 9 pixels.
;    THRESHOLD:        Set this keyword to adjust the threshold used in detecting zingers and 
;                    diffraction artifacts. The defaults are 1.2 for zingers and 0.8 for 
;                    diffraction artifacts.
;
; OUTPUTS:
;    This function returns an image with the specified artifacts removed.
;
; PROCEDURE:
;    THIS STILL NEEDS TO BE DOCUMENTED. For now, see the source code.
;
; EXAMPLE:
;    r = remove_tomo_artifacts(sinogram, /Rings)
;
; MODIFICATION HISTORY:
;    Written by:    Mark Rivers, May 13, 1998
;-

; Copy input image to output
output = image
size = size(image)
ncols = size[1]
nrows = size[2]

if (keyword_set(zingers)) then begin
    ; Default kernal width for smoothing
    if (n_elements(width) eq 0) then width=9
    ; Default threshold for detecting zingers
    if (n_elements(threshold) eq 0) then threshold=1.5
    ; Compute ratio of unsmoothed image to smoothed image
    ratio = float(image)/smooth(image, width)
    ; Find indices of pixels which are above threshold
    zinger = where(ratio gt threshold, count)
    ; Replace pixels with zingers by the average of the pixels 2 to the left and
    ; 2 to the right. We don't use nearest neighbor, since zingers sometimes hit 2
    ; adjacent pixels.
    if (count gt 0) then output[zinger] = $
        (output[zinger-2] + output[zinger+2])/2.
    endif

if (keyword_set(rings)) then begin
    ; Default kernal width for smoothing
    if (n_elements(width) eq 0) then width=9
    ; Sum all of the rows in the image, get average
    ave = total(image, 2)/nrows
    ; Get the difference between the average row and smoothed version
    ; of average row. These are the column deviations.
    diff = ave - smooth(ave, width)
    ; Subtract this difference from each row
    for i=0, nrows-1 do output[0,i] = output[*,i] - diff
endif

if (keyword_set(diffraction)) then begin
    ; Default kernal width for smoothing
    if (n_elements(width) eq 0) then width=9
    ; Default threshold for detecting diffraction peaks
    if (n_elements(threshold) eq 0) then threshold=.8
    ; Loop over each column in the image
    for i=0, ncols-1 do begin
        col = reform(output[i,*])
        ratio = col/smooth(col, width)
        bad = where(ratio lt threshold)
        col[bad] = (col[(bad-2)>0] + col[(bad+2)<(ncols-1)]) / 2.
        output[i,*] = col
    endfor
endif

return, output
end