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

(1) Pr[N(A)] = ^nspect∏_k = 1(e^{− λ_k(A)}λ^n_k_k(A))/(n_k!)

The measurement vector N has components n_k, which are the measurements with spectra k = 1…nspect. The expected values of the n_k are λ_k(A) given by

(2) λ_k(A) = ⟨n_k⟩ = ⌠⌡S_k(E)e^{− A⋅μ(E)}dE, k = 1, …, nspect

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 Â_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

(3) ℒ(A) = ^nspect⎲⎳_k = 1n_klog(λ_k(A)) − λ_k(A) − log(n_k!).

Notice that the measurement data samples n_k 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

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 Bä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.