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.
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.
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.