function sinogram, input, angles, $
                acc_values = acc_values, $
                air_values = air_values, $
                backlash = backlash, $
                center = center, $
                tweak_center = tweak_center, $
                cog = cog, $
                debug = debug

;+
; NAME:
;        SINOGRAM.PRO
; PURPOSE:
;        To convert raw tomography data into a sinogram.
; CATEGORY:
; CALLING SEQUENCE:
;        result = SINOGRAM(INPUT, ANGLES, [keywords])
; INPUTS:
;    INPUT
;        An array of raw tomography data. INPUT(I, J) is the intensity at
;        position I for view angle J. Each row is assumed to contain at least
;        one air value at each end for normalization.
;    ANGLES
;        An array of the angles of each row of the input. Units are radians.
; OPTIONAL INPUT PARAMETERS:
;    NONE
; KEYWORD PARAMETERS:
;    ACC_VALUES=acc_values
;        The number of values to be discarded at the beginning and end of each
;        row. These values are typically discarded because the stage was in its
;        acceleration/deceleration phase or because there are simply an
;        unecessarily large number of air values at the ends of each row.
;        This keyword is only useful for first-generation CT data.
;        The default value is 0.
;    AIR_VALUES=air_values
;        The number of air values to be averaged together at the beginning and
;        and of each row, after discarding the ACC_VALUES. This averaging is
;        done to decrease the statistical uncertainty in the air values.
;        The default value is 10.
;    COG=cog
;        This keyword is used to return the measured and fitted
;        center-of-gravity data for the sinogram. The center-of-gravity data are
;        very useful for diagnosing problems such as backlash, beam hardening,
;        detector saturation, etc. COG is dimensioned (n_angles, 2), where
;        n_angles is the number of angles in the input array. COG(*,0) is the
;        measured center-of-gravity data. COG(*,1) is the fitted data. The
;        following command can then be given after the SINOGRAM command
;            IDL> PLOT, COG(*,0)
;            IDL> OPLOT, COG(*,1)
;        to see is the sinogram data are reasonable.
;    /CENTER
;        Used to specify that the output is to be shifted left or right so that
;        the fitted rotation axis of the sinogram is centered exactly on the
;        center column of the output array. This operation is used to correct
;        for the fact that the rotation axis may not have been perfectly
;        centered when the data were collected.
;        The default is not to center the output.
;
;    TWEAK_CENTER=tweak_center
;        Used to specify a "tweak" to the center. This tweak value is a floating
;        point value (in pixels) which is added to the center which is determined automatically
;        by the /CENTER keyword.
;
;    /BACKLASH
;        Used to specify that even rows of the image are to be shifted left or
;        right so that the fitted rotation axis of the even rows is the same
;        as that for the odd rows. This is used to correct for backlash on
;        images done with first generation CT scans, when the scanning is
;        done bidirectionally.
;        The default is not to correct backlash.
;    /DEBUG
;        Used to turn on debugging output. The default is not to print debugging output.
;
; RETURN:
;    The output array containing the corrected sinogram. It is always of
;    type FLOAT.
;
; PROCEDURE:
;        This routine creates a sinogram from raw tomography data. It does the
;        following:
;    -    Converts to an odd number of columns (if necessary) by discarding the
;        last column
;    -    Discards unwanted or uneeded pixels from the left and right edges of the
;        image. These could be values collected during motor acceleration or
;        extra air values which will simply slow down the reconstruction.
;        "ACC_VALUES" pixels are discarded from both the left and right edges
;        of the array.
;    -    Averages the air values for "air_values" pixels on the left and right
;        hand sides of the input.
;    -    Logarithmation. output = -log(input/air). The air values are
;        interpolated between the averaged values for the left and right hand
;        edges of the image for each row.
;    -    Backlash correction (optional) If /BACKLASH is specified then motor
;        backlash is corrected for. This is done by fitting the
;        center-of-gravity separately for the even and odd rows of the image
;        and then shifting the even rows so the the rotation axes are the same.
;        The measured center-of-gravity is fitted to a sinusoidal curve
;        of the form Y = A(0) + A(1)*SIN(X + A(2)).
;            A(0) is the rotation axis
;            A(1) is the amplitude
;            A(2) is the phase
;        The fitting is done using routine CURVE_FIT in the User Library.
;        The shifting is done using routine POLY_2D which can shift by
;        fractional pixels.
;    -    Centering of the rotation axis (optional). If /CENTER is specified then
;        the image is shifted so that the rotation axis obtained by fitting the
;        center-of-gravity data for all rows in the image coincides exactly
;        with the center column. If TWEAK_CENTER is specified then this value is
;        added to the fitted center before the shifting is done.
;        The shifting is done using routine POLY_2D which can shift by
;        fractional pixels.
; MODIFICATION HISTORY:
;        Created 21-OCT-1991 by Mark Rivers.
;        MLR 25-NOV-1994:
;            Removed BSIF_COMMON
;            Made TOMO_HEADER a structure passed to the routine
;            This structure must contain the fields .ANGLE, .WIDTH, .STEP_SIZE
;        MLR 13-MAY-1998
;            Converted from a procedure to a function, added ANGLES parameter, removed
;            TOMO_HEADER parameter. Added TWEAK_CENTER and DEBUG keywords.
;-

if n_elements(air_values) eq 0 then $
    air_values = 10
if n_elements(acc_values) eq 0 then $
    acc_values = 0
if n_elements(tweak_center) eq 0 then $
    tweak_center = 0.

ncols = n_elements(input(*,0))
nrows = n_elements(input(0,*))
; If there are an even number of columns in the image throw out the last one
if (ncols mod 2) eq 0 then ncols=ncols-1

; Throw out the acceleration values, convert data to floating point
output = float(input(acc_values:ncols-acc_values-1, *))
ncols = n_elements(output(*,0))

cog = fltarr(nrows) ; Center-of-gravity array
linear = findgen(ncols) / (ncols-1)
lin2 = findgen(ncols) + 1.
weight = fltarr(nrows) + 1.
for i=0, nrows-1 do begin
    air_left = total(output(0:air_values-1,i)) / air_values
    air_right = total(output(ncols-air_values:ncols-1,i)) / air_values
    air = air_left + linear*(air_right-air_left)
    output(0,i) = -alog(output(*,i)/air > 1.e-5)
    cog(i) = total(output(*,i) * lin2) / total(output(*,i))
endfor
odds = where((indgen(nrows) mod 2) ne 0)
evens = where((indgen(nrows) mod 2) eq 0)
x = angles
a = [(ncols-1)/2., $ ; Initial estimate of rotation axis
    (max(cog) - min(cog))/2., $ ; Initial estimate of amplitude
    0.] ; Initial estimate of phase
cog_fit = curvefit(x(odds), cog(odds), weight(odds), a, sigmaa, $
                    function_name='sine_wave')
cog_odd = a(0) - 1.
if (keyword_set(debug)) then print, format='(a, f8.2, a, f8.2)', $
'Fitted center of gravity for odd rows = ', cog_odd, ' +-', sigmaa(0)
cog_fit = curvefit(x(evens), cog(evens), weight(evens), a, sigmaa, $
                    function_name='sine_wave')
cog_even = a(0) - 1.
if (keyword_set(debug)) then print, format='(a, f8.2, a, f8.2)', $
    'Fitted center of gravity for even rows = ', cog_even, ' +-', sigmaa(0)
if (n_elements(backlash) ne 0) then begin
    back = cog_even - cog_odd
    if (keyword_set(debug)) then print, format='(a,f8.2)', $
        'Backlash (even rows shifted right relative to odd rows) = ', back
    P = [[back, 0.],[1., 0.]]
    Q = [[0., 1.],[0., 0.]]
    temp = poly_2d(output(*, evens), P, Q, 1)
    for i=0, nrows-1, 2 do begin
        output(0, i) = temp(*, i/2)
        cog(i) = total(output(*,i) * lin2) / total(output(*,i))
    endfor
endif
cog_fit = curvefit(x, cog, weight, a, sigmaa, $
                    function_name='sine_wave')
cog_mean = a(0) - 1.
error_before = cog_mean - (ncols-1)/2.
if (keyword_set(debug)) then print, format='(a, f8.2, a, f8.2)', $
    'Fitted center of gravity = ', cog_mean, ' +-', sigmaa(0)
if (keyword_set(debug)) then print, format='(a,f8.2)', $
        'Error before correction (offset from center of image) = ', error_before
if (n_elements(center) ne 0) then begin
    P = [[error_before+tweak_center, 0.],[1., 0.]]
    Q = [[0., 1.],[0., 0.]]
    output = poly_2d(output, P, Q, 1)
    for i=0, nrows-1 do begin
    cog(i) = total(output(*,i) * lin2) / total(output(*,i))
endfor
    cog_fit = curvefit(x, cog, weight, a, sigmaa, $
                        function_name='sine_wave')
    cog_mean = a(0) - 1.
    error_after = cog_mean - (ncols-1)/2.
    if (keyword_set(debug)) then print, format='(a, f8.2, a, f8.2)', $
        'Fitted center of gravity after correction= ', cog_mean, ' +-', sigmaa(0)
    if (keyword_set(debug)) then print, format='(a,f8.2)', $
        'Error after correction (offset from center of image) = ', error_after
endif

cog = [[cog], [cog_fit]]

if (keyword_set(debug)) then print, "Sinogram used average of "+string(air_values)+" pixels for air"
if (keyword_set(debug)) then print, "Skipped "+string(acc_values)+" acceleration pixels"
if (n_elements(backlash) ne 0) then begin
    if (keyword_set(debug)) then print, "Backlash corrected "+string(back)+" pixels"
endif
if (n_elements(center) ne 0) then begin
    if (keyword_set(debug)) then print, "Center corrected "+string(error_before)+" pixels"
endif

return, output
end