If you have trouble viewing this, try
the pdf of this post. You can
download the code used to produce the figures in this post.
Edited May 12, 2012 to clarify that constant covariance log data and photon count data give same formula for the CRLB.
SNR with PHA vs. number of bins
In this post, I will discuss the signal to noise ratio (SNR) of a photon counting detector with pulse height analysis (PHA). I will show that it approaches the ideal fullspectrum SNR as the number of bins gets large. I will also show that we can get quite close to the ideal value even with a small number of bins, which was one of the main points in
my paper.
I will begin by revisiting a result from my
last post that the CramèrRao lower bound (CRLB) with multivariate normal log of photon count data assuming a constant covariance is nearly equal to the accurate matrix that includes the variation of the covariance. The derivation becomes problematic with a large number of bins since the mean value in each bin may be small enough so the normal approximation to the counts is not valid. I will show an alternate derivation that uses the Poisson model of the count data, which is valid for small counts. The result is formally the same as with the normal approximation so the CRLB formula can be applied for any number of bins. This result is new and was not included in
my paper.
Imaging task for SNR with Aspace data
Summarizing my previous development, to determine the SNR with energy spectrum information I analyzed the imaging task shown in Fig.
1↓. I introduced the Aspace description by using the attenuation coefficients of the background
μ_{b} and feature
μ_{f} materials as the basis function. With this basis set, the
a vectors for the feature and background materials are simply
a_{f} = ⎡⎢⎣
1
0
⎤⎥⎦, a_{b} = ⎡⎢⎣
0
1
⎤⎥⎦.
The Avector in the background region is
A_{b} = t_{b}a_{b} while in the feature region it is
A_{f} = (t_{b} − t_{f})a_{b} + t_{f}a_{f} . We can apply the statistical detection theory model described in
my previous post by assuming that the Avector data are described by a multinormal distribution with expected values under two hypotheses:
⟨A_{H0}⟩ = A_{b} and
⟨A_{H1}⟩ = A_{f} and a covariance
C_{A}. From my last post, the performance will depend only on the signal to noise ratio
where
δA = ⟨A_{H1}⟩ − ⟨A_{H0}⟩ = t_{f}(a_{f} − a_{b}).
δA
=
⟨A_{H1}⟩ − ⟨A_{H0}⟩
=
t_{f}(a_{f} − a_{b})
=
t_{f}⎡⎢⎣
1
− 1
⎤⎥⎦
In general, the covariance of the Aspace data depends on the estimator used. But we know that the ’minimum’ covariance (suitably defined for matrix data) is the CRLB, so we can determine the optimal performance by using it for the covariance
C_{A}. In a
previous post, I showed that we could use a linearized model with multivariate normal noise data. With this model, the CRLB is
where
M = (∂L)/(∂A) is the gradient of log of the xray measurements in Aspace evaluated at the operating point of the model and
C_{L} is the covariance of the log data with noise. The SNR with the imaging task in Fig.
1↓ is
CRLB with complete energy information
I would like to apply this model to show that the SNR approaches the full spectrum as the number of bins becomes large. However, as mentioned, with a large number of bins the data may not satisfy the assumptions of a large number of counts. In that case, I will use a Poisson model, which is valid for small mean counts, to show that the same formula (
2↑) can be used in all cases.
I will derive the CRLB for independent multivariate Poisson data. As I have discussed in
a previous post, this is a good model for the PHA counts if the deadtime is negligible compared to the mean interarrival time of the photons. Table
1↓ lists the notation I use in this post.
Table 1 Notation
photon counts

N = ⎡⎢⎢⎢⎢⎢⎣
n_{1}
n_{2}
⋯
n_{K}
⎤⎥⎥⎥⎥⎥⎦


mean counts

⟨N⟩ = ⎡⎢⎢⎢⎢⎢⎣
⟨n_{1}⟩
⟨n_{1}⟩
⋯
⟨n_{K}⟩
⎤⎥⎥⎥⎥⎥⎦ = ⎡⎢⎢⎢⎢⎢⎣
g_{1}( A)
g_{2}( A)
⋯
g_{K}(A)
⎤⎥⎥⎥⎥⎥⎦ = g(A)


covariance

C = ⎡⎢⎢⎢⎣
⟨n_{1}⟩
0
0
0
⋯
0
0
0
⟨n_{K}⟩
⎤⎥⎥⎥⎦


atten. coeff.

M = ⎡⎢⎢⎢⎢⎢⎣
μ_{11}
μ_{12}
μ_{21}
μ_{22}
⋯
⋯
μ_{K1}
μ_{K2}
⎤⎥⎥⎥⎥⎥⎦


It is convenient to use the Fisher information matrix
F, which is the inverse of the CRLB. As shown in Theorem 3.2 of Kay
[1],
F has components which are the negative of the expected value of the second derivative of the log of the likelihood
ℒ = log[Pr(N;A)]
F_{ij} = − ⟨(∂^{2}ℒ)/(∂A_{i}∂A_{j})⟩
The bin counts
n_{k} are independent Poisson random variables, which have a distribution function
Pr(n_{k};A) = (e^{ − ⟨nk⟩}⟨n_{k}⟩^{nk})/(n_{k}!) = (e^{ − gk(A)}g_{k}(A)^{nk})/(n_{k}!)
where the expected value ⟨n_{k}⟩ = g_{k}(A) is assumed to be a function of the Avector. Since the n_{k} are independent the joint probability is the product of the distributions and the log of the likelihood is
ℒ = ^{K}⎲⎳_{k = 1}[ − g_{k} + n_{k}log(g_{k}) − log(n_{k}!)].
To evaluate the Fisher information, we need the derivatives. The first derivative can be expressed as a matrix product of the gradient of
g, the inverse of the photon count covariance
C^{ − 1}, and the difference of the measured data and the expected values. Since the count data are independent, their covariance matrix and its inverse are
C = ⎡⎢⎢⎢⎣
g_{1}
0
⋯
0
g_{K}
⎤⎥⎥⎥⎦, C^{ − 1} = ⎡⎢⎢⎢⎣
^{1}⁄_{g1}
0
⋯
0
^{1}⁄_{gK}
⎤⎥⎥⎥⎦
(∂ℒ)/(∂A_{j})
=
⎲⎳⎡⎣(∂g_{k})/(∂A_{j})⎛⎝(n_{k} − g_{k})/(g_{k})⎞⎠⎤⎦
=
(∂g^{T})/(∂A_{j})⎡⎢⎢⎢⎣
^{1}⁄_{g1}
0
⋯
0
^{1}⁄_{gK}
⎤⎥⎥⎥⎦(N − g)
=
(∂g^{T})/(∂A_{j})C^{ − 1}(N − g)
Note that this covariance is valid even with small number of counts per bin.
The second derivative is
(∂^{2}ℒ)/(∂A_{i}∂A_{j}) = (∂)/(∂A_{i})⎡⎣(∂g^{T})/(∂A_{j})C^{ − 1}⎤⎦(N − g) − (∂g^{T})/(∂A_{j})C^{ − 1}(∂g)/(∂A_{i}).
The expected value of the first term is zero since the term in the bracket does not depend on the data and⟨N⟩ = g. The second term does not depend on the data so
F_{ij}
=
− ⟨(∂^{2}ℒ)/(∂A_{i}∂A_{j})⟩
=
(∂g^{T})/(∂A_{j})C^{ − 1}(∂g)/(∂A_{i})
The full matrix is
where I have used the notation of van den Bos
[2] Appendix D so
(∂g)/(∂A^{T}) = ⎡⎢⎢⎢⎣
(∂g_{1})/(∂A_{1})
(∂g_{1})/(∂A_{2})
⋯
(∂g_{2})/(∂A_{1})
(∂g_{2})/(∂A_{2})
⋯
⋯
⋯
⋯
⎤⎥⎥⎥⎦
and ⎛⎝(∂g)/(∂A^{T})⎞⎠^{T} = (∂g^{T})/(∂A).
If the number of bins is large so the energy width is small, we can use the Beer’s law exponential model of the expected value of the counts
g_{k}(A) = ⟨n_{k}⟩ = n_{k0}exp[ − μ_{k1}A_{1} − μ_{k2}A_{2}]
where
n_{k0} is the mean of the number of photons incident on the object in the
k’th energy bin. Differentiating,
(∂g_{k})/(∂A_{i}) = − μ_{ki}g_{k}.
We can express all the components as a matrix multiplication
(∂g)/(∂A^{T}) = − M^{T}C
and the Fisher matrix (
4↑) is
where I have used the fact that all covariance matrices are symmetrical so C^{T} = C.
The Fisher information matrix derived in my
previous post assuming log data is
F = M^{T}C^{ − 1}_{L}M (see Eq. 8 of my previous post, which I reproduce here)
Recall that the CRLB is the inverse of
F, so the Fisher information matrix from my previous post is the term in the parentheses of Eq.
6↑. With photon count data, the variance is the inverse of the number of counts and since
C_{L} is a diagonal matrix,
C^{ − 1}_{L} = C. Thus Eq.
5↑ has the same form as the Fisher information of constant covariance multivariate normal data, derived in my previous post but here it is derived without the assumptions that the covariance is constant and the mean count values are large. Therefore, we can use the same formula (
2↑) for the covariance of the Aspace data for any number of bins.
SNR for K bin PHA
To derive the SNR with Kbin data, we need to evaluate Eq.
3↑ with
δA = t_{f}⎡⎢⎣
1
− 1
⎤⎥⎦, M = ⎡⎢⎢⎢⎣
M_{11}
M_{12}
⋯
⋯
M_{K1}
M_{K2}
⎤⎥⎥⎥⎦, C = ⎡⎢⎢⎢⎣
⟨n_{1}⟩
0
⋯
0
⟨n_{K}⟩
⎤⎥⎥⎥⎦.
This is straightforward but tedious, so I used Maxima to automate the matrix computation. Maxima is an opensource symbolic math program that is available for download
here. A good introduction is available
here. I develop the code by writing a batch script with a text editor and then running it using the Maxima
batch command. For example, I copy this command to the clipboard and then paste it into the WxMaxima GUI. With modern computers, the evaluation is very fast and we can see the results almost immediately
batch("E:\\Projects\\Blog\\Posts\\P29D2kbin\\code\\SNRkbin.max")$
A Maxima script to compute the SNR with a 3 bin detector is
/*
SNRkbin.max
Maxima batch file to derive SNR^2 for K bin PHA
to invoke copy and paste into Maxima:
batch("E:\\Projects\\Blog\\Posts\\P29D2kbin\\code\\SNRkbin.max")$
REA 2012Feb14 11:38
COPYRIGHT: © Aprend Technology, 2012. All rights reserved.
*/
"clear all workspace variables"$
kill(all)$
"The SNR^2 USING 3 bin PHA"$
"the M matrix with general coeffs and its transpose"$
M:matrix([m11,m12],[m21,m22],[m31,m32])$
MT:transpose(M);
"the covariance matrix"$
C:diagmatrix(3,0);
"fill the diagonal"$
S:matrix([n1,n2,n3]);
for k:1 thru 3 do C[k,k]:S[1,k]$
"display c"$
C;
"The delta_A vector for basis materials feature and background"$
dA:matrix([1],[1])$
dAT:transpose(dA);
"the SNR^2"$
d2:ratsimp(dAT.MT.C.M.dA);
"factor the individual terms"$
d2simp:n1*(m12m11)^2 + n2*(m22m21)^2 + n3*(m32m31)^2;
"verify correct"$
ratsimp(d2d2simp);
The initial result was (we can use edit”copy as Latex” in the Maxima GUI to copy expressions and then paste them into a Latex or Lyx file)
(m32^{2} − 2 m31 m32 + m31^{2}) n3 + (m22^{2} − 2 m21 m22 + m21^{2}) n2 + (m12^{2} − 2 m11 m12 + m11^{2}) n1.
The program did not factor the quadratics in the parentheses but this is easily done. The SNR with 3 bins generalizes to
K bins as
As explained in my paper, the factor of
^{1}⁄_{2} is introduced to make the imaging task SNR comparable to the detection theory SNR, which assumes only one measurement. The notation
⟨μ⟩_{k} means the weighted mean of the function
μ(E) in the spectrum of bin
k. For example, if bin
k corresponds to energies from
E_{k} to
E_{k} + Δ_{E} and
S(E) is the spectrum incident on the detector, then the weighted mean is
⟨μ⟩_{k} = (^{Ek + ΔE}⌠⌡_{Ek}μ(E)S(E)dE)/(^{Ek + ΔE}⌠⌡_{Ek}S(E)dE)
Using the notation
δμ = μ_{2}(E) − μ_{1}(E), for the difference in attenuation coefficient functions, Eq.
7↑ is the same as Eq. 35 of
my paper
SNR^{2}_{K} = (t^{2}_{f})/(2)^{K}⎲⎳_{k = 1}⟨n_{k}⟩⟨δμ⟩^{2}_{k}.
Simulate Kbin SNR
I wrote Matlab code to compute the SNR as a function of the number of PHA bins. I normalized the result by dividing by the ideal, fullspectrum SNR. I also computed the results as a function of object thickness. Recall that the Avectors for this kind of object lie along a straight line in the Aplane. The results are shown in Fig.
2↓.
Conclusions
There are two interesting features of Fig.
2↑. First, as we expected, the SNR approaches the ideal result as the number of bins gets larger. This is shown mathematically in my paper using the mean value theorem. The second interesting aspect is that the performance with a small number of bins depends strongly on the object thickness. My interpretation of this result is that as the object thickness increases the spectral width also decreases due to beam hardening. Therefore, there is less energy information and it is easier for a low energyresolution detector to extract it. The improvement occurs for relatively thin objects. From the figure, the threebin SNR with a 5.5 cm object is about 93% of the ideal value.
Last edited May 12, 2012
© 2012 by Aprend Technology and Robert E. Alvarez
Linking is allowed but reposting or mirroring is expressly forbidden.
References
[1] Steven M. Kay: Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory . Prentice Hall PTR, 1993.
[2] Adriaan van den Bos: Parameter Estimation for Scientists and Engineers. WileyInterscience, 2007.