next up previous contents
Next: EM (expectation maximation) Up: Iterative methods Previous: Iterative methods

ART (algebraic reconstruction technique)

This is an extension of the Kaczmarz (1937) method for solving linear systems. It has been introduced in imaging by Gordon et al. (1970). We describe it in a more general context. We consider the linear system

  equation1941

where tex2html_wrap_inline4092 are bounded linear operators from the Hilbert space H into the Hilbertspace tex2html_wrap_inline4096 . With tex2html_wrap_inline4098 a positive definite operator we define an iteration step tex2html_wrap_inline4100 as follows:

  eqnarray1945

If tex2html_wrap_inline4102 (provided tex2html_wrap_inline4104 is onto) and tex2html_wrap_inline4106 , then tex2html_wrap_inline4108 is the orthogonal projection of tex2html_wrap_inline4110 onto the affine subspace tex2html_wrap_inline4112 . For dim tex2html_wrap_inline4114 this is the original Kaczmarz method. Other special cases are the Landweber method (p = 1, tex2html_wrap_inline4118 and fixed-block ART (dim tex2html_wrap_inline4120 finite, tex2html_wrap_inline4122 diagonal, see Censor and Zenios (1997)). It is clear that there are many ways to apply (4.3) to (4.1), and we will make use of this freedom to our advantage.

One can show (Censor et al. (1983), Natterer (1986)) that (4.3) converges provided that

displaymath4124

and tex2html_wrap_inline4126 . This is reminiscent of the SOR-theory of numerical analysis. In fact we have tex2html_wrap_inline4128 where tex2html_wrap_inline4130 is the k-th SOR iterative for the linear system tex2html_wrap_inline4134 with

displaymath4136

If (4.2) is consistent, ART converges to the solution of (4.2) with minimal norm in H.

Plain convergence is useful, but we can say much more about the qualitative behaviour and the speed of convergence by exploiting the special structure of the image reconstruction problems at hand. With R the Radon transform in tex2html_wrap_inline3076 we can put

displaymath4144

displaymath4146

where w is the weight function tex2html_wrap_inline4150 and tex2html_wrap_inline4152 . One can show that the subspaces

displaymath4154

displaymath4156

tex2html_wrap_inline4158 the Gegenbauer polynomials of degree m (Abramowitz and Stegun (1970)) are invariant subspaces of the iteration (4.3). This has been discovered by Hamaker and Solmon (1978). Thus it suffices to study the convergence on each subspace tex2html_wrap_inline4162 separately. The speed of convergence depends drastically on tex2html_wrap_inline4164 and - surprisingly - on the way the directions tex2html_wrap_inline4166 are ordered. We summerize the findings of Hamaker and Solmon (1978), Natterer (1986) for the 2D case where tex2html_wrap_inline4168 .

1. Let the tex2html_wrap_inline4170 be ordered consecutively, i.e. tex2html_wrap_inline4172 , tex2html_wrap_inline3236 . For tex2html_wrap_inline4164 large (e.g. tex2html_wrap_inline4106 ) convergence on tex2html_wrap_inline4162 is fast for m < p large and slow for m small. This means that the high frequency parts of f (such as noise) are recovered first, while the overall features of f become visible only in the later iterations. For tex2html_wrap_inline4164 small (e.g. tex2html_wrap_inline4192 ) the opposite is the case.

This explains why the first ART iterates for tex2html_wrap_inline4106 look terrible and why surprisingly small relaxation factors (e.g. tex2html_wrap_inline4192 ) are used in tomography.

2. Let tex2html_wrap_inline4170 be distributed at random in tex2html_wrap_inline4200 and tex2html_wrap_inline4106 . Then the convergence is fast on all tex2html_wrap_inline4162 , m < p. The same is true for more sophisticated arrangements of the tex2html_wrap_inline4170 , such as for p = 18 (Hamaker and Solmon (1978))

displaymath4212

or similarly for p = 30 in Herman and Meyer (1993).

The practical consequence is that a judicious choice of directions may well speed up the convergence tremendously. Often it suffices to do only 1-3 steps if the right ordering is chosen. This has been observed also for the EM iteration (see below).

Note that the tex2html_wrap_inline4162 with tex2html_wrap_inline4218 are irrelevant since they describe those details in f which cannot be recovered from p projections because they are below the resolution limit. This is a result of the resolution analysis in Natterer (1986).


next up previous contents
Next: EM (expectation maximation) Up: Iterative methods Previous: Iterative methods

Frank Wuebbeling
Thu Sep 10 10:51:17 MET DST 1998