An iterative method is derived for image reconstruction. Among other attributes, this method allows constraints unrelated to the radiation measurements to be incorporated into the reconstructed image. A comparison is made with the widely used Maximum-Likelihood Expectation-Maximization (MLEM) algorithm.
Deep Dive into Iterative method for solution of radiation emission/transmission matrix equations.
An iterative method is derived for image reconstruction. Among other attributes, this method allows constraints unrelated to the radiation measurements to be incorporated into the reconstructed image. A comparison is made with the widely used Maximum-Likelihood Expectation-Maximization (MLEM) algorithm.
arXiv:1101.0819v1 [physics.comp-ph] 4 Jan 2011
Iterative method for solution of radiation emission/transmission matrix equations
Clinton DeW. Van Siclen∗
Idaho National Laboratory, Idaho Falls, Idaho 83415, USA
(Dated: 4 January 2011)
An iterative method is derived for image reconstruction. Among other attributes, this method
allows constraints unrelated to the radiation measurements to be incorporated into the reconstructed
image. A comparison is made with the widely used Maximum-Likelihood Expectation-Maximization
(MLEM) algorithm.
Imaging by radiation emission or transmission effec-
tively produces a set of linear equations to be solved. For
example, in the case of coded aperture imaging, the so-
lution is a “reconstructed” set of radiation sources, while
in the case of x-ray interrogation, the solution is a set
of attenuation coefficients for the voxels comprising the
volume through which the x-ray beam passes.
The linear equations have the form
di =
J
X
j=1
Mijµj
(1)
where the set {di} corresponds to the radiation intensity
distribution recorded at a detector (a detector pixel is
labeled by the index i), the set {µj} is the solution, and
the matrix element Mij connects the known di to the un-
known µj. Typically the matrix M is non-square so that
{µj} cannot be obtained by standard matrix methods.
(And note that, when the set of equations is large, it can
be difficult to ascertain a priori whether the equation set
is over- or under-determined.)
In any case the matrix equation d = Mµ may be
solved by the iterative method that is derived as follows.
Clearly this method will feature a relation between µ(n)
j
and µ(n−1)
j
, where n is the iteration number. Consider
the two equations for µ(n)
j
and µ(n−1)
j
,
d(n)
i
= P
jMijµ(n)
j
(2)
d(n−1)
i
= P
jMijµ(n−1)
j
(3)
and rewrite the latter as
di =
di
d(n−1)
i
P
jMijµ(n−1)
j
.
(4)
Then the relationship between µ(n)
j
and µ(n−1)
j
is ob-
tained by setting P
id(n)
i
= P
idi:
P
i
P
jMijµ(n)
j
= P
i
di
d(n−1)
i
P
jMijµ(n−1)
j
!
P
j
n
µ(n)
j
P
iMij
o
= P
j
(
µ(n−1)
j
P
i
di
d(n−1)
i
Mij
!)
µ(n)
j
= µ(n−1)
j
1
P
iMij
P
i
di
d(n−1)
i
Mij
!
.
(5)
Note that this last equation can be written
µ(n)
j
= µ(n−1)
j
di
d(n−1)
i
Mij
i
⟨Mij⟩i
(6)
where the last factor is essentially a weighted average of
all di/d(n−1)
i
. Thus the set
n
µ(n)
j
o
approaches a solution
{µj} by requiring P
id(n)
i
= P
idi at each iteration; in
effect, by requiring all d(n)
i
→di.
The iteration procedure alternates between use of Eq.
(3) and Eq. (5) until all d(n)
i
are as close to di as desired.
For the first (n = 1) iteration, an initial set
n
µ(0)
j
o
is
chosen, which produces the set
n
d(0)
i
o
according to Eq.
(3).
These values are used in Eq.
(5), so producing
the set
n
µ(1)
j
o
. And so on. . . That a final set
n
µ(n)
j
o
is a solution {µj} to the matrix equation d = Mµ is
verified by checking that all d(n)
i
= di to within a desired
tolerance.
Some cautions and opportunities follow from this sim-
ple derivation of Eq. (5). A caution is that, in the event
the equation set is under-determined, different initial sets
n
µ(0)
j
o
will lead to different final sets {µj} that satisfy
the matrix equation. The corresponding opportunity is
that this problem may be mitigated to some extent by
the addition, to the original set of equations, of linear
equations that further constrain the µj (perhaps derived
from, for example, independent knowledge of some of the
contents of a container under interrogation). In general
the di appearing in a constraint equation will have noth-
ing to do with radiation intensity.
The form of any added constraints, and the initial
choice
n
µ(0)
j
o
, must allow all µ(n)
j
→µj and d(n)
i
→di
monotonically. In particular, care should be taken when
a constraint has one or more coefficients Mij < 0, as that
affects the denominator P
iMij in Eq. (5) (a straightfor-
ward fix may be to reduce the magnitudes of all Mij
coefficients and di in that constraint equation by a mul-
tiplicative factor). In any event, the acceptability of a
2
constraint equation is easily ascertained by monitoring
the behavior d(n)
i
→di for that constraint.
Note that all solutions {µj} to a set of equations
that includes additional constraints with di > 0 and all
Mij ≥0 are accessible from sets
n
µ(0)
j
o
of initial values,
and further that any set
n
µ(0)
j
o
will produce a solution
{µj}. This suggests that, for this implementation of con-
straints, a superposition of many solutions may give a
good “probabilistic” reconstruction. To achieve this, con-
sider that the innumerable solutions to the set of equa-
tions may be regarded as points in a J-dimensional space
(J is the number of elements in a solution {µj}). These
points must more-or-less cluster, producing a cluster cen-
troid that is itself a solution. While the centroid solution
n
µ(c)
j
o
has no intrinsic special status (as all cluster points
represent equally likely reconstructions), it may be taken
to represent the particular set of equations. The clu
…(Full text truncated)…
This content is AI-processed based on ArXiv data.
