Estimator for contrast agents-2 The iterative estimator
As its name implies, the maximum likelihood estimate is the value of the dependent variable that maximizes the likelihood given the measured data. One way to implement it is to use an iterative algorithm, which I discussed
here. In this post, I give a detailed a description of the code for an iterative estimator. The implementation is different than the one used in the previous post and is included as
AsolveIterFromSpectrum.m in the
code package .
Log-likelihood of Poisson distributed photon count data
Assuming that the PHA detector counts are independent and Poisson distributed, the likelihood function, which is the probability of a particular measurement
N considered to be a function of
A, is
The measurement vector N has components nk, which are the measurements with spectra k = 1…nspect. The expected values of the nk are λk(A) given by
where μ(E) is a vector whose components are the basis functions μj(E), j = 1…nbasis evaluated at energy E and the symbol ⋅ is the vector dot product.
The maximum likelihood estimate
ÂMLE is the value of
A that maximizes Eq.
1↑. It is convenient and equivalent to maximize the logarithm of the likelihood since the log is a monotonically increasing function
Notice that the measurement data samples
nk do not depend on
A so we can ignore the last term in summand of the log-likehood Eq.
3↑ during the maximization.
Code to compute log-likelihood and maximum likelihood estimate
The code for the iterative estimator
AsolveIterFromSpectrum.m is included with the package for this post. The following functions are used to compute the log-likelihood. Actually the negative is used so we can compute the maximum using Matlab’s
fminsearch function. The code is a straightforward implementation of Eq.
3↑. The
logPoissProbFuncLoc function calls the helper function
LambdasLoc to implement Eq.
2↑ to compute the
λk(A).
function logp = logPoissProbFuncLoc(A,solvedat)
% (negative of) exponent of Poisson pdf for given A.
% The measurement is counts for each PHA bin, passed in as solvedat.nspnt
lambdas = LambdasLoc(solvedat.specdat,A);
% use negative of log likelihood so get maximum from fminsearch
logp = -sum(solvedat.nspnt.*log(lambdas) -lambdas);
function lambdas = LambdasLoc(sp,As)
% the integral (number of counts) of the spectrum in each of the PHA bins
specnum_transmitted = sp.specnum(:).*exp(-sp.mus*As(:));
lambdas = zeros(sp.nbins,1);
for kbin = 1:sp.nbins
lambdas(kbin) = sum(specnum_transmitted(sp.idxs2bins{kbin}));
end
The log-likelihood for each PHA bin photon count measurement vector is maximized by the following code. The solvedat struct is pre-loaded with the incident spectrum and the basis functions attenuation coefficients. The counts for a particular measurement are loaded into the struct at the first line of code. The second line calls the fminsearch function. Notice the use of a Matlab anonymous function so the solvedat struct is passed to logPoissProbFuncLoc during the calls by fminsearch.
solvedat.nspnt = tdat.ns(kpnt,:)’;
% pass counts in as member of spectrum struct so
% spectrum and mus accessible to logProbFuncLoc
As=fminsearch(@(A) logPoissProbFuncLoc(A,solvedat),dat.As_initfit(kpnt,:),options);
dat.Astest(kpnt,:) = As’;
The initial guess
dat.As_initfit(kpnt,:) is either passed to the
AsolveIterFromSpectrum.m function by the calling program or computed as the linear maximum likelihood estimate.
Data required for estimator
Studying the code for the log-likehood function in the previous section, the data required to implement the iterative estimator are:
-
spectrum incident on the object, sp.specnum
-
PHA bin response functions, sp.idxs2bins
-
attenuation coefficients of basis materials, mus.
Schlomka et al.
[1] measured all of these for their experimental implementation.
Measuring the spectrum with clinical x-ray tubes is difficult because of the high photon count rates. Schlomka et al. measured the spectrum using a cooled germanium sensor with a high-energy resolution multi-channel analyzer. This was possible only because their experimental scanner used a micro-focus x-ray tube with a current of 0.1 milliAmperes (mA). A typical clinical CT tube current is hundreds of times larger so some method would be needed to reduce the count rate. Filtering, of course, is not feasible because it affects the spectrum so some experimenters have tried to use separation from the tube. Regardless of the method used to lower the count rate, the high resolution energy spectrum analyzer that Schlomka et al. used is usually not found in clinical institutions and is quite expensive.
The PHA bin response functions were measured by Schlomka et al. using a synchrotron accelerator at the German DESY laboratory as a nearly monochromatic source. This is an admirable achievement but it requires access to a major nuclear research laboratory and is clearly not practical for a commercial clinical x-ray system. Even if it could be contemplated to measure the responses during the manufacture, they could drift over time once the system is installed at a clinical institution therefor requiring periodic re-measurement.
The tabulated x-ray attenuation coefficient databases are quite accurate so using the computed attenuation coefficients would probably give useful results in terms of the virtual materials.
Last edited Mar 07, 2016
Copyright © 2016 by Robert E. Alvarez
Linking is allowed but reposting or mirroring is expressly forbidden.
References
[1] J P Schlomka, E Roessl, R Dorscheid, S Dill, G Martens, T Istel, C Bäumer, C Herrmann, R Steadman, G Zeitler, A Livne, R Proksa: “Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography”, Phys. Med. Biol., pp. 4031—4047, 2008.