Algebraic Reconstruction Methods

Algebraic reconstruction methods The problem of solving for the density (actually, for the linear attenuation coefficient) of each location in the image can also be viewed as a set of simultaneous equations. Each ray integral (or summation, in the finite case we are dealing with here) provides one equation. The sum of the attenuation coefficients for the pixels (or voxels) along the ray, each multiplied by a weighting factor that takes into account the actual path length of that ray through the pixel, is equal to the measured absorption. Figure 12.13 illustrates the relationship between the pixels and the ray-integral equations. The number of unknowns in this set of equations is the number of pixels in the image of the slice through the specimen. The number of equations is the number of ray integrals, which is generally the number of detectors used along each projection profile times the number of view angles. This is a very large number of equations, but many of the weights are zero (most pixels are not involved in any one particular ray-integral equation). Furthermore, the number of equations rarely equals the number of unknowns. But fortunately there are a number of practical and well-tested computer methods for solving such sets of sparse equations when they are under- or over determined. It is not our purpose here to compare the various solution methods. A suitable understanding of the method can be attained using the simplest of the methods, known as the algebraic reconstruction technique or ART (Gordon 1974). In this approach, the equations are solved iteratively. The set of equations can be written as Amn xn = bm where n is the number of voxels, m is the number of projections, and A is the matrix of weights that correspond to the contribution of each voxel to each ray path (which can be precalculated for any particular instrument and geometry). The voxel values are the x values and the projection measurements are the b values. The classic ART method calculates each iterative set of x values from the preceding ones as xk+1 = xk +Ai bi Ai x A 2k ( ± λ ) || i || The value of λ, the relaxation coefficient, generally lies between 0 and 2, and controls the speed of convergence. When λ is very small, this becomes equivalent to a conventional least squares solution. Practical considerations, including the order in which the various equations are applied, are dealt with in detail in the literature (Censor 1983, 1984). Figure 12.14 shows a simple example of this approach. The 16 × 16 array of voxels has been given density values from 0 to 20 as shown in Figure 12.14b, and three projection sets at view angles of 0, 90, and 180° were calculated for the fan-beam geometry shown in Figure 12.14a. For an array of 25 detectors, this gives a total of 75 equations in 256 unknowns. Starting with an initial guess of uniform voxels (with density 10), the results after 1, 5, and 50 iterations are shown. The void areas and internal square appear rather quickly, and the definition of boundaries gradually improves. The errors — particularly in the corners of the image, where fewer ray equations contain any information, and at the corners of the internal dense square, where the attenuation value changes abruptly — are evident. Still, considering the extent to which the system is underdetermined, the results are rather good. Kacmarz’s method for this type of

solution is illustrated in Figure 12.15 for the very modest case of three equations and two unknowns, with λ = 1. Beginning at some initial guess, for instance that all of the pixels have the same attenuation value, one of the equations is applied. This is equivalent to moving perpendicular to the line representing the equation. This new point is then used as a starting point to apply the next equation, and so on. In the real case, the equations do not all meet in a perfect point because of finite precision in the various measurements, counting statistics, machine variation, etc.; thus there is no single point that represents a stable answer. Instead, the solution converges toward a region that is mostly within the region between the various lines and then oscillates there. However, in a high-dimensionality space with some noisy equations, it is possible for the solution to leave this region and wander away after much iteration. In real cases with many dimensions, the convergence may not be very fast. The greatest difficulty in using the iterative algebraic technique is deciding when to stop. Logically, we would like to continue until the answer is as good as it can get, but without knowing the “truth,” it is not possible to determine this stopping point exactly. Some methods examine the change in the calculated image after each iteration and attempt to judge from that when to stop (for instance, when the normalized total variation in pixel values falls below some arbitrary limit, or when it begins to increase from the previous iteration). This method is prone to serious errors in a few cases, but is used nonetheless. It should be noted that the penalty for continuing the iteration is not simply the computational cost, but also the possibility that, for some sets of data, the answer may start to diverge (leave the bounded region near the crossover point). This condition is, of course, highly undesirable. Given the drawbacks to the algebraic approach and the relative simplicity and straight forward approach of the filtered back-projection method, why would we use this method? There are several potential advantages of algebraic methods such as ART. First, the filtered back-projection method, and the Fourier transform method that it embodies, require that the number of views be rather large and that they be equally spaced so that the frequency space is well filled with data. Missing angles, or entire sets of angles that may be unattainable due to physical limitations, present problems for filtered back-projection and introduce significant artifacts. ART methods can still produce an acceptable reconstruction. There may be a lack of detail in portions of the reconstructed image that are under sampled by the projections, but the artifacts do not spread throughout the entire image. In fact, acceptable reconstructions are often obtained with only a very few views. Another advantage to ART is the ability to apply constraints. For instance, it is possible, in a filtered back-projection or Fourier transform method, to calculate negative values of density (attenuation) for some voxels because of the finite measurement precision. However, such values have no physical meaning. In the iterative algebraic method, any such values can be restricted to zero. In the schematic diagram of Figure 12.15, this amounts to restricting the solution to the quadrant of the graph with positive values. In fact, any other prior knowledge can also be applied. If it is known that the only possible values of density and attenuation in the specimen correspond to specific materials, then the values can be easily constrained to correspond. Any geometric information, such as the outside dimensions of the object, can also be included (in this case, by forcing the voxels outside the object boundaries to zero density).

It is even possible to set up a grid of voxels that are not all of the same size and spacing. This Set up might allow, for instance, the use of a fine voxel spacing in the interior of an object, where great detail is desired, but a much coarser grid outside (or vice versa). This would still allow the calculation of the contribution of the outside material to the ray integrals, but it would reduce the number of unknowns to produce a better solution for any given number of views and projections. Sets of non square pixels or non cubic voxels can also be used when these are needed to conform to specific object shapes and symmetries. The flexibility of the algebraic method and its particular abilities to use a priori information, often available in an industrial tomography setting, compensates for its slowness and requirements for large amounts of computation. The calculation of voxel weights (the A matrix) can be tedious, especially for fan-beam or other complex geometries, but no more so than back projection in such cases, and it is a one-time calculation whose results can be stored and used for many reconstructions using the same geometry. The use of solution methods other than the simple iterative approach described here can provide improved stability and convergence.

Shiven Gandhi SKG Creative Designing | Organizing | Overseas Educational Counseling (Biomedical)