Iterative method for solution of radiation emission/transmission matrix equations

0.07 MB·English

by  Clinton DeW. Van Siclen

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
allowsconstraintsunrelatedtotheradiationmeasurementstobeincorporatedintothereconstructed
image. AcomparisonismadewiththewidelyusedMaximum-LikelihoodExpectation-Maximization
(MLEM) algorithm.
1
1
0 Imaging by radiation emission or transmission effec- (n) (n−1) 1 di
µ =µ M . (5)
n 2 teixvaemlyppler,odinuctehseacsaesteoofflicnoedaerdeqaupaetritounrsetiombaegisnoglv,etdh.eFsoor- j j iMijPi di(n−1) ij!
a lutionisa“reconstructed”setofradiationsources,while Note that this last equPation can be written
J in the case of x-ray interrogation, the solution is a set

4 of attenuation coefficients for the voxels comprising the
di M
] volume through which the x-ray beam passes. µ(n) =µ(n−1)(cid:28)di(n−1) ij(cid:29)i (6)
h The linear equations have the form j j M
p h ijii
J
- where the last factor is essentially a weighted average of
p d = M µ (1)
m i ij j alld /d(n−1). Thusthe set µ(n) approachesa solution
Xj=1 i i j
s.co dwihsterriebuthtieosnetre{cdoir}dceodrraetspaonddetsetcototrhe(araddeiatetciotnorinptiexneslitiys e{ffµejc}t,bbyyrreeqquuiririninggPalilddi(i(nn))=n→Pdii.odi at each iteration; in
c labeled by the index i), the set µj is the solution, and The iteration procedure alternates between use of Eq.
{ }
si thematrixelementMij connectstheknown di totheun- (3)andEq. (5)untilalldi(n) areasclosetodi asdesired.
y known µj. Typicallythe matrixM is non-squaresothat For the first (n = 1) iteration, an initial set µ(0) is
h µ cannot be obtained by standard matrix methods. j
j
p ({An}dnote that,whenthe setofequationsislarge,itcan chosen, which produces the set d(0) accordinng tooEq.
[ i
bedifficulttoascertaina priori whethertheequationset (3). These values are used innEq. o(5), so producing

1 is over- or under-determined.) the set µ(1) . And so on... That a final set µ(n)
v In any case the matrix equation d = Mµ may be j j
9 solvedby the iterative method that is derivedasfollows. is a soluntionoµj to the matrix equation d =nMµ ios
81 Clearly this method will feature a relation between µj(n) verifiedbyche{ckin}gthatalldi(n) =di towithina desired
0 and µ(n−1), where n is the iteration number. Consider tolerance.
. j Some cautions andopportunities follow fromthis sim-
01 the two equations for µ(jn) and µj(n−1), ple derivation of Eq. (5). A caution is that, in the event
1 theequationsetisunder-determined,differentinitialsets
d(n) = M µ(n) (2)
1 i j ij j µ(0) will lead to different final sets µ that satisfy
: j { j}
v P nthe moatrix equation. The corresponding opportunity is
Xi di(n−1) = jMijµj(n−1) (3) that this problem may be mitigated to some extent by
the addition, to the original set of equations, of linear
r and rewrite the latter as P
a equations that further constrainthe µ (perhaps derived
j
d = di M µ(n−1). (4) from,forexample,independentknowledgeofsomeofthe
i d(n−1) j ij j contents of a container under interrogation). In general
i
P the di appearingin aconstraintequationwillhavenoth-
Then the relationship between µ(n) and µ(n−1) is ob- ing to do with radiation intensity.
j j
tained by setting d(n) = d : The form of any added constraints, and the initial
i i i i (0) (n) (n)
choice µ , must allow all µ µ and d d
P P j j → j i → i
M µ(n) = di M µ(n−1) monotonnicalloy. In particular, care should be taken when
i j ij j i d(n−1) j ij j ! aconstrainthasoneormorecoefficientsMij <0,asthat
(cid:16) (cid:17) i
P P P P affectsthe denominator M inEq. (5)(astraightfor-
i ij
ward fix may be to reduce the magnitudes of all M
ij
µ(n) M = µ(n−1) di M coefficients and di in thaPt constraint equation by a mul-
j j i ij j( j i d(n−1) ij!) tiplicative factor). In any event, the acceptability of a
n o i
P P P P
2
constraint equation is easily ascertained by monitoring may be regardedas the sought-after reconstruction.
the behavior d(n) d for that constraint. Anothercautionfollowsfromthefactthatthe denom-
i → i
Note that all solutions µj to a set of equations inator iMij inEq. (5)isthesumofallelementsincol-
that includes additional con{stra}ints with d > 0 and all umn j of matrixM. This iterativemethod canofcourse
i P
M 0 are accessible from sets µ(0) of initial values, beusedwithoutexplicitlyconvertingasetoflinearequa-
ij ≥ j tions into amatrix equation(orseveral sets ofequations
and further that any set µ(0) nwill poroduce a solution into a single matrix equation), but in that event very
j
µj . Thissuggeststhat,nforthoisimplementationofcon- careful attention must be paid to get the factors iMij
s{tra}ints, a superposition of many solutions may give a right.
P
good“probabilistic”reconstruction. Toachievethis,con- It may be noticed that Eq. (5) is similar to the so-
sider that the innumerable solutions to the set of equa- called Maximum-Likelihood Expectation-Maximization
tionsmayberegardedaspointsinaJ-dimensionalspace (MLEM)algorithm(seerefs. [1]and[2]forderivationsof
(J is the number of elements in a solution µ ). These the latter, and see numerous papers in the recent imag-
j
pointsmustmore-or-lesscluster,producing{aclu}stercen- ing literature for applications of it). The MLEM, which
troid thatisitselfasolution. Whilethecentroidsolution isderivedfromphysicalconsiderationshavingtodowith
µ(c) hasnointrinsicspecialstatus(asallclusterpoints radiation emission and detection, purports to find the
j
set µ that maximizes the probability P ( d µ ),
rnepresoentequallylikelyreconstructions),itmaybetaken { j} { i}|{ j}
whichistheprobabilityofrealizingtheobservedset d
to represent the particular set of equations. The cluster { i}
given a set µ . This is in contrast to the derivation
size, which indicates the degree to which solutions are { j}
above,whichleadstoEq. (5)assimply amethodtofind
similar to one other, should decrease as constraints are
asolutiontoamatrixequation(asetoflinearequations).
added. A logical measure of the cluster size is
Thederivationpresentedheremakesclearhowtheitera-
tiveprocedureshouldbeimplementedforanapplication,
σ = (x x ) (x x ) 1/2 (7)
cluster h c− k · c− k ik and allows constraints to be added to the original set
of equations (those produced by the imaging exercise)
where x is the centroid vector and x is the vector
c k thereby enabling a more-accurate reconstruction when
corresponding to the kth solution. Thus the quantity
the original equation set is under-determined.
σ /√J, which represents the standard deviation of
cluster This work was supported in part by the INL Labora-
the innumerable values of an arbitrary element µ , is a
j tory Directed Research and Development Program un-
usefulmeasureofthe variationamongsolutions µ . In
j der DOE Idaho Operations Office Contract DE-AC07-
{ }
generalitis desirablethat the variationamong solutions
05ID14517.
be muchless thanthe variationwithin the centroidsolu-
tion, which is
2 1/2
σµ(c) = µ(jc)−µ(c) (8) ∗ [email protected], [email protected]
(cid:28)(cid:16) (cid:17) (cid:29)j [1] L. A. Shepp and Y. Vardi, “Maximum likelihood recon-
struction for emission tomography,” IEEE Trans. Med.
whereµ(c) = µ(c) . Inthatcase(σ J−1/2 σ(c)) Imag. MI-1, 113 (1982).
j j cluster ≪ µ [2] K. Lange and R. Carson, “EM reconstruction algorithms
the centroid sDolutioEn is little changed by additional con- for emission and transmission tomography,” J. Comput.
straints, so suggesting that the centroid solution µ(c) Assist. Tomography 8, 306 (1984).
j
n o