 |
| |
The RIEMANN procedure computes the "Riemann sum" (or its inverse) which helps implement the backprojection operator used to reconstruct the cross-section of an object, given projections through the object from multiple directions. This technique is widely used in medical imaging in the fields of computed x-ray tomography, MRI imaging, Positron Emission Tomography (PET), and also has applications in other areas such as seismology and astronomy. The inverse Riemann sum, which evaluates the projections given a slice through an object, is also a discrete approximation to the Radon transform.
Given a matrix A ( m,n ), which will contain the reconstructed slice; a vector P , containing the ray sums for a given view; and an angle Theta measured in radians from the vertical: the Riemann sum "backprojects" the vector P into A . For each element of A , the value of the closest element of P is summed, leaving the result in A . Bilinear interpolation is an option. All operations are performed in single-precision floating point.
In the reverse operation, the ray sums contained in the view vector, P , are computed given the original slice, A , and Theta . This is sometimes called "front projection".
The Riemann sum can be written:
which is the sum of the data along lines through an image with an angle of theta from the vertical.
A k -element floating-point projection vector (or matrix if the ROW keyword is specified). For backprojection (when the BACKPROJECT keyword is set), P contains the ray sums for a single view. For the inverse operation, P should contain zeros on input and will contain the ray sums for the view on output.
Set this keyword to perform backprojection in which P is summed into A . If this keyword is not set, the inverse operation occurs and the ray sums are accumulated into P .
Set this keyword to use bilinear interpolation rather than the default nearest neighbor sampling. Results are more accurate but slower when bilinear interpolation is used.
Set this keyword equal to a floating-point number specifying the center of the projection. The default value for CENTER is one-half the number of elements of P .
Set this keyword equal to a two-element floating-point vector specifying the center of rotation in the array A . The default value is [ m /2., n /2.], where A is an m by n array.
For symmetric results, given symmetric operands, COR should be set to the origin of symmetry [( m -1)/2, ( n -1)/2], and CENTER should be set to ( n -1)/2., where n is the number of elements in the projection vector, P .
Set this keyword to a value between -1 and 0 to use the cubic convolution interpolation method with the specified value as the interpolation parameter. Setting this keyword equal to a value greater than zero specifies a value of -1 for the interpolation parameter. Park and Schowengerdt (see reference below) suggest that a value of -0.5 significantly improves the reconstruction properties of this algorithm.
Cubic convolution is an interpolation method that closely approximates the theoretically optimum sinc interpolation function using cubic polynomials. According to sampling theory, details of which are beyond the scope of this document, if the original signal, f , is a band-limited signal, with no frequency component larger than w 0 , and f is sampled with spacing less than or equal to 1/2 w 0 , then f can be reconstructed by convolving with a sinc function: sinc ( x ) = sin ( p x ) / ( p x ).
In the one-dimensional case, four neighboring points are used, while in the two-dimensional case 16 points are used. Note that cubic convolution interpolation is significantly slower than bilinear interpolation.
Rifman, S.S. and McKinnon, D.M., "Evaluation of Digital Correction Techniques for ERTS Images; Final Report", Report 20634-6003-TU-00, TRW Systems, Redondo Beach, CA, July 1974.
S. Park and R. Schowengerdt, 1983 "Image Reconstruction by Parametric Cubic Convolution", Computer Vision, Graphics & Image Processing 23, 256.
Use this keyword to specify the spacing between elements of P , expressed in the same units as the spacing between elements of A . The default is 1.0.
Set this keyword to specify the P vector as a given row within a matrix, so that the sinogram array can be used directly without having to extract or insert each row. In this case, P must be an array with a first dimension equal to k , and the value of ROW must be in the range of 0 to the number of vectors of length k in P , minus one.
This example forms a synthetic image in A , computes Nviews equally spaced views, and stores the stacked projections (commonly called the "sinogram") in a matrix PP. It then backprojects the projections into the matrix B, forming the reconstructed slice. In practical use, the projections are convolved with a filter before being backprojected.
N = 100L ; Define number of columns in A.
M = 100L ; Define number of rows in A.
nviews = 100 ; Number of views.
K = CEIL(SQRT(N^2 + M^2)) ; The length of the longest projection. If filtered backprojection is used, 1/2 the length of the convolution kernel must also be added.
A = FLTARR(N, M) ; Form original slice.
A(N/2:N/2+5, M/2:M/2+5) = 1.0 ; Simulate a square object.
pp = FLTARR(K, nviews) ; Make array for sinogram.
FOR I=0, NVIEWS-1 DO RIEMANN, pp, A, I * !PI/nviews, ROW=i
;
Compute each view.
B = FLTARR(N,M) ; Initial reconstructed image.
FOR I=0, nviews-1 DO $ ; Do the backprojection for each view.
