Solution methods for the resulting equations (original) (raw)

Next: Reconstruction from noisy projections Up: Algebraic Reconstruction Methods Previous: Attenuation Correction

Depending on the number of projections collected, the system of equations to be solved may be either over- or underdetermined. In our case, due to limitations on the amount of memory space, the system was underdetermined (W had fewer rows than columns). The problem was therefore reformulated as an optimization problem,

Here, denotes an estimate of the intensity vector , is a LaGrange multiplier, and the matrix Q is chosen in order to control the behavior of the solution. Setting the derivative of the above to zero, we obtain the equation

For the purposes of this project, it was found that choosing Q as the identity matrix produced satisfactory results, although other choices were also investigated (see the section on Wiener filtering). By choosing Q as the identity matrix, we obtain the solution to (16) that has the smallest norm, i.e., a reasonably smooth solution.

Given the size of the system matrix (e.g., in the case of a image and as used in this study), it is not feasible to use methods such as Gauss elimination or LU decomposition to solve this system. This is due both to the large number of operations required to invert or factor such a matrix, and to the fact that the inverted or factored matrix would be more dense, i.e., utilize more memory than the original matrix. In fact, it is desirable to use a method which does not even require that the matrix product be computed explicitly. For this project, we chose to use the conjugate gradient method [5,6], which is an iterative method frequently used in solving large systems of equations arising, e.g., in the context of numerical solution of partial differential equations. The conjugate gradient method only requires that one is able to compute the result of multiplying by a vector. Thus, may be computed as , which avoids computing (and storing) the denser . Computation times on the Sparc Ultras on Owlnet were generally less than 10 minutes (CPU time) for a reconstruction from 120 projections.

Next: Reconstruction from noisy projections Up: Algebraic Reconstruction Methods Previous: Attenuation Correction

Anders Johan Nygren Thu May 8 12:28:25 CDT 1997