Nonlinear Dual Reconstruction of SPECT Activity and Attenuation Images

In single photon emission computed tomography (SPECT), accurate attenuation maps are needed to perform essential attenuation compensation for high quality radioactivity estimation. Formulating the SPECT activity and attenuation reconstruction tasks as coupled signal estimation and system parameter identification problems, where the activity distribution and the attenuation parameter are treated as random variables with known prior statistics, we present a nonlinear dual reconstruction scheme based on the unscented Kalman filtering (UKF) principles. In this effort, the dynamic changes of the organ radioactivity distribution are described through state space evolution equations, while the photon-counting SPECT projection data are measured through the observation equations. Activity distribution is then estimated with sub-optimal fixed attenuation parameters, followed by attenuation map reconstruction given these activity estimates. Such coupled estimation processes are iteratively repeated as necessary until convergence. The results obtained from Monte Carlo simulated data, physical phantom, and real SPECT scans demonstrate the improved performance of the proposed method both from visual inspection of the images and a quantitative evaluation, compared to the widely used EM-ML algorithms. The dual estimation framework has the potential to be useful for estimating the attenuation map from emission data only and thus benefit the radioactivity reconstruction.


Introduction
Single photon emission computed tomography (SPECT) has become an indispensable tool in clinical trials and medical practice. Attenuation correction in SPECT has significant values to better understand the physiological processes associated with the disease (i.e. cancers, heart diseases) and to provide improvement in patient diagnosis and treatment. However, even with ample efforts devoted to the development of various attenuation estimation and compensation techniques, it remains one of the key open issues in SPECT imaging [1][2][3][4][5][6][7][8][9][10][11].
Tissue attenuation map is usually estimated based on transmission data by scanning the patient with a rotating external radionuclide source [3,5,6,12], or obtained from X-ray computed tomography (CT) system [10,11,[13][14][15][16]. However, transmission based attenuation correction clearly increases the patient's dose, and requires maintaining additional radioactive sources. Further, if multiple imaging sessions are needed, it may be difficult for some patients to tolerate for a longer scan at one time, and leads to coregistration problem in emission image reconstruction, especially for deformed tissues and organs. In addition to the added equipment cost and the well-known beam hardening problem, similar registration issue exists for CT attenuation data.
It has been the goal of many recent research efforts to simultaneously estimate the activity and attenuation distributions from emission data. Iterative statistical methods have been extensively studied, with the main incentive being that they explicitly take into account the specific SPECT data statistics. Some of the most notable works include the use of differential attenuation method [17], gradient ascent [18], Tikhonov regularization [19], and expectation maximization (EM) [20][21][22][23]. With the SPECT imaging and measurement processes in state space representation, a recent work has adopted the extended Kalman filtering (EKF) procedures to linearize the augmented state representation to provide the joint estimates in the minimummean-square-error sense [24]. Such EKF based framework, however, has several potential drawbacks. First, the derivation of the Jacobian matrices, the linear approximations to the nonlinear functions, often leads to filter instability. Furthermore, it has been shown that estimation bias may originate from a coupling between the state variables and the model parameters, which suggests that if the system parameters can be separated from the system state variables, the precision of the estimation could be improved.
Following the same spirit while addressing the limitations, we present a robust unscented Kalman filter framework for the joint estimation of SPECT activity and attenuation parameters. Instead of linearizing using the Jacobian matrices and thus overcoming the associated shortcomings, our effort deals with the nonlinear estimation process by using a deterministic sampling approach to capture the mean and covariance estimates [25]. In addition, two iteratively coupled filters are used to sequentially estimate the activity variables (with fixed attenuation estimates from the last iteration) and the attenuation map (with fixed activity values from the previous estimation). Furthermore, such framework can explicitly recognize the uncertainties in the measurement data and the model structure, thus has the potential to produce more robust reconstructions.

SPECT Emission Scan Model
In SPECT imaging, once radio-tracers are injected into the subjects, they are delivered to the tissues/organs by the blood flow and participate in the related physiologic/metabolic processes. The general mathematical model for the emission scan can be stated as where column vector x(t)~fx j Dj~1,:::,Ng is the radiopharmaceutical concentration of the object, e o is the background noises (e.g. scatter events), and y is the emission measurement data that is acquired by a rotating detector head around the patient at each angle and represented lexicographically as a column vector y(t)~fy i Di~1,:::,Mg. Here, i indicates different projection defined by rotating angle and different detector bin, and M is the total number of projections. G represents the emission system matrix including attenuation effects, and its element indicates the probability that a photon emitted from pixel gets detected in a specific detector bin. Let the attenuation map be given by column vector m~fm k Dk~1,:::,Ng, we may explicitly account for the attenuation effects in the system matrix G by factoring it as Here, The symbol ' : ' means the dot product of two matrices. ½L ijk~l ijk represents the length of the ray through voxel k with respect to a photon emitted from voxel i being detected in projection bin j. The photon survival probability considering attenuation can now be represented by a matrix A with elements ½A ij~a ij as The term D, with elements ½D ij~d ij , is the photon-detection probability without taking into account the attenuation effects.
The mathematical model of the SPECT emission measurement (Eqn. (1)) can now be rewritten as or in discrete form In the following sections, we will present the dual estimation method which works by alternating between using the UKF statefilter to estimate the radioactivity based on fixed attenuation parameters, and the UKF parameter-filter to estimate the attenuation coefficients given previously estimated activity states. Although two separate state-space representations are constructed for the activity and attenuation estimation problems, both of them use the same emission scan model (Eqn. (4) or (5)) as the measurement equation.

Activity Estimation: UKF State Filter
State Space Representation for Radioactivity. In emission tomography, the goal is to reconstruct the radioactivity distribution x(t) from the measurement data y(t). The system equation of the SPECT imaging system, which describes the radioactivity evolution of the pixels, can be written in the form of with initial activity x 0 and system noise u s that accounts for the statistical uncertainty of the imaging model. In general, Eqn. (6) represents the dynamic changes of the state variable x(t), and it reduces to the conventional static reconstruction problem when the transition matrix S is an identity. The associated measurement equation, which describes the observations provided by the imaging data y(t), is given by Eqn. (4). And now we have the following state space representation for radioactivity distribution: where u s and e o , with covariance matrices Q s and R o , model the uncertainties of the imaging system and the measurement data respectively. Two important observations on the nonlinear SPECT imaging system represented by Eqns. (7) and (8). First, since the noises in emission sinogram are typically Poisson distributed, it is difficult to perform standard estimation on such non-Gaussian system. By applying the Anscombe transformation [26], however, the Poisson noise could be converted into a Gaussian one, and various Kalman filtering techniques become viable options. Secondly, the EKF, probably the most widely used estimation algorithm for nonlinear systems, has several drawbacks. EKF requires the system is almost linear on the time scale of the updates, and is difficult to implement and to set proper parameters due to the cross-talk between state and parameter. To overcome such limitations, we have adopted the unscented transformation (UT) to accurately propagate mean and covariance information through nonlinear transformations [25], with little additional computational cost.
UKF State-Filter. For the SPECT imaging system given by Eqns. (7) and (8), the purpose of the UKF state-filtering is to search for the optimal estimates of radio-tracer variable x(t), given fixed attenuation values m(t). Since detailed discussion of the unscented Kalman filtering is certainly beyond the scope of this paper [25], we present here a more algorithmic description while ignoring some theoretical considerations.
The unscented transformation is a method for calculating the statistics of a random variable. Given a random variable x (dimension L) with mean x x and covariance P x through a nonlinear function y~f(x). To compute the statistics of y, we form a sigma point matrix X and their weights according to the following: 2(Lz ) ,i~1,:::,2L where~a 2 (Lzk){L. The parameter a[½0,1 determines the spread of the sigma points, b §0 is used to incorporate any prior knowledge about the distribution of x, and k is a scaling parameter and is often set to be zero. Once having the above definitions, the UKF state-filter is initialized withx x(0)~x x 0 and covariance matrix P x(0) , the state estimates and their error covariance matrices are computed sequentially until convergence: 1. Calculate the sigma point weights as in (9); 2. Project the state variable x(t) ahead: 3. Project the error covariance P x(t) ahead: 4. Calculate the sigma points as defined in Eqn. (9): 5. Filtering of the measurement equations:   6. Compute the Kalman Gain: 7. Update the estimate with the measurement: 8. Update the error covariance: Here, the previously stored information in the prediction step is combined with the new information coming from the next measurement y(t) and the Kalman gain matrix to refinex x(t) and P x(t) in the correction step. The covariance of the measurement error R o and system error Q s is assumed to be known and set to time-invariant.

Attenuation Estimation: UKF Parameter Filter
Once the UKF state-filter converges, it is followed by the estimation of a coupled UKF parameter-filter aiming to recover the attenuation map of the object being imaged, given the estimated radioactivity map. The system equations for the parameter-filter are m(tz1)~m(t)zu p ð20Þ Here, u p is the process noise with covariance matrix Q p , and we assume that the attenuation parameter vector m is temporally  constant. Initializing the unscented parameter filter withm m(0)~m m 0 and covariance matrix P m(0) , the parameter-filter follows similar recursion steps as the state-filter (Eqns. (9)-(17)) until convergence, given the sigma point calculation scheme of Eqn. 9. The coupled radioactivity state and attenuation parameter estimation processes are iteratively repeated as necessary, until stable results are achieved. The final optimal estimates then become the reconstructed SPECT activity and attenuation maps.

Validation with Synthetic Data
Synthetic Zubal phantom is used to quantitatively evaluate the accuracy and robustness of the framework, where a simplified human thorax with two lungs is represented by 32|32 attenuation and activity (Fig. 1). Specifically, while both lungs have average attenuation coefficient of 0.04/cm, attenuation of the left lung is nonuniform but that of the right lung is uniform. SPECT projection data have been generated for a parallel beam geometry and 90 views uniformly spaced over 360 degrees. To generate realistic data, simulations in our study are performed using toolbox GATE [27], which can provide a relatively accurate reference for the assessment of new image reconstruction algorithms. And we have performed two studies, one with mean activities of 50 counts/pixel (low count), and the other with 200 counts/pixel (high count).
For these two sets of synthetic projection data, the radioactivity maps are reconstructed using four reconstruction methods: FBP  [28], EM-ML [29], EKF [24], and dual UKF (The attenuation maps are also reconstructed for the EKF and UKF methods). For the UKF framework, the state (activity map) is estimated by the UKF state-filter, where, for the first iteration, the attenuation parameters come from the initial guess. Otherwise, we use the values estimated in the last iteration from the UKF parameterfilter. Consequently, these state estimates are used by the UKF parameter-filter to recover the attenuation coefficients. This process is executed iteratively until it meets the convergence criterion, which is defined using two consecutive normalized errors x(tz1) and x(t) through Ex(tz1){x(t)Evz with z being a small constant, and x defines the normalized error between the estimated and the exact value (J~x for the state-filter process, J~m for the parameter-filter process) with where J t i is the estimated value, J r i is the corresponding true value, and i indicates the pixel.
A detailed statistical analysis on the estimation results against the ground truth phantom map is performed. Let N p andx x be the total number of pixels in the region of interest(ROI), e.g. body anatomy, and the final reconstruction results respectively, and x tr be the ground truth, we have the following error definitions: The results for the reconstruction of the attenuation maps are shown in Fig. 2 for the low count case. Visually, the UKFreconstructed attenuation image clearly shows all relevant anatomical structures, where the lungs are easily seen with clear shape and size. The RMSE values in ROI of body anatomy using the UKF method are 0:0190 (low counts) and 0:0187 (high  counts), which, though not as impressive, are still somewhat smaller than the corresponding values of 0:0204 (low counts) and 0:0199 (high counts) obtained by the EKF method (note that the true average attenuation coefficient is 0.04). The worse performance of the EKF strategy could be caused by the first-order Taylor approximations of state transition that provided insufficiently accurate representations.
The results for the reconstruction of the activity maps are presented in Fig. 3, and the quantitative tabulation of the reconstruction accuracy is in Table 1. These figures and results illustrate that traditional EM-ML and FBP methods, with unknown attenuation map settings, produce some noticeable errors. The dual UKF estimation framework, on the other hand, consistently yields the best quality radioactivity estimates for both high and low count data. Same conclusion can be drawn from the visual examples of the selected horizontal profiles, as shown in Fig. 4.

Reconstruction from Physical Phantom Scanning Data
The second data set used for validation has been on a real cylinder phantom. The dimension of the phantom is 200 mm (diameter) | 290 mm (depth). A Teflon rod and two hollow PMMA cylinders with diameters of 50mm are inserted in the phantom's volume, as shown in Fig. 1. The phantom is filled with 99m Tc concentration with a total radioactivity of 20mCi (100kBq/ cc) and the two hollow cylinder rods are filled with air and pure water respectively. The phantom was scanned with a Siemens ECam duet ECT scanner by two detector head rotating at total 64 angle position around 180 degree and the acquiring time at each position is 30 seconds. The final sinogram data has 64|128 projections for each slice. Once again, FBP, EM-ML, EKF and dual UKF strategies have been used to reconstruct activity maps from the measurement data, as shown in Fig. 5. Visually, it is evident that the UKF method produces the best reconstruction results, especially for the three cold areas. The results of EKF and UKF reconstructed attenuation map are shown in Fig. 6, with RMSE values of 0:0007 for UKF and 0:0044 for EKF.

Reconstruction from Real Patient Scanning Data
The dual reconstruction strategy has also been evaluated on clinical studies, where the patients are undergoing 99m Tc sestamibi stress tests. Using a Siemens ECam duet scanner, all projections are acquired over 120 angles covering a circular 360 degrees acquisition orbit in a continuous step-and-shoot mode. With an acquisition time of 16s/frame, the total photon counts for each slice are 148761. The dual reconstruction of activity and attenuation maps are shown in Fig. 7, and the clinically standard FBP reconstruction (without attenuation correction) is also shown for comparison. Since the transmission data is not available from this imaging site, we can only make qualitative visual inspection of the images. It is quite clear, however, that the estimated attenuation map agrees with general knowledge of the imaged area, and the UKF estimated activity map exhibits improved contrast between heart and soft tissue.
There is usually a considerable increase in computation for improved performance. The computational load increases when moving from the EM to the UKF. However, as the UKF gives a better approximation in time update step, the UKF estimate is able to converge quite faster comparing to the EM. Furthermore, this proposed approach runs efficiently on graphics processing units(GPUs) since large amounts of computations are done in matrix forms. Further investigations on the implementation with GPUs are underway.

Conclusions
A dual UKF strategy has been derived for joint reconstruction of the attenuation map and activity distribution solely from SPECT emission sinograms. Constructing the state transition of the activity distribution through state space evolution equations and the photon-counting measurements through observation equations, we rely on the unscented Kalman filter principles to first generate estimates of activity maps with sub-optimal attenuation parameter estimates, and then recover the attenuation maps given these activity estimates. These coupled iterative steps are repeated as necessary until convergence. Simulated and physical phantoms, as well as real patient data, are used to evaluate the proposed strategy.