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