Scan time optimization for post-injection PET scans

zy zyxw Scan Time Optimization for Post-injection PET Scans Hakan Erdogan and Jeffrey A. Fessler 4415 EECS Bldg., University of Michigan, Ann Arbor, MI 48109-2122 Abstract methods. Particularly we are interested in the optimum scan time fractions under a fixed total scan time constraint, which Previous methods for optimizing the scan times for PET would result in the smallest variance in a region of interest in transmission and emission scans under a total scan time the final emission image estimate. Previous studies of scan constraint were based on linear non-statistical methods and time optimization [2] were based on NEC criteria with multiple used noise equivalent counts (NEC) criteria. The scan acquisitions of emission and transmission data and focused on times determined by NEC analysis may be suboptimal when conventional reconstructions. The (co)variance approximations nonlinear statistical image reconstruction methods are used. developed here might also be useful for other purposes such For statistical image reconstruction, the predicted variance in as determining the weights in a weighted least-squares image selected regions of interest is an appropriate alternative to NEC reconstruction [3]. We analyze both the conventional and analysis. We propose a new method for optimizing the relative statistical reconstruction cases. We give approximate analytical scan times (fractions) based on analytical approximations to formulas for conventional and quadratic penalty attenuation the covariance of images reconstructed by both conventional map reconstructions and compare empirical results with the and penalized-likelihood methods. We perform simulations analytical predictions. Our analysis is based on Poisson to compare predicted standard deviations with empirical statistics and mathematical approximations [4]. ones. Results show that for statistical transmission image Let yE = [yF . . . and yT = [y: . . . yz]’ be emission reconstruction, the optimal fraction of the scan time devoted and post-injection transmission scan count vectors, and let to transmission scanning is shorter than for conventional p = [/11 . . .pp]’ and A = [A, . . .Ap]’ be attenuation map and transmission smoothing. emission image pixel value vectors respectively. zyxwv zyx zy zy YE]’ I. INTRODUCTION We define the survival probabilities as follows: cti(p) = where l i ( p ) represents the line integral along projection i of the attenuation map p. We also define the emission contamination count rate KC^ (A, p ) = Li ai ai( p ) p i (A). Here hi is the fraction of emission counts contaminating the transmission data (the portion in the transmission window for rotating rod sources), pi(A) represents the geometric projections of the true emission image A, and ai contains the detector efficiencies and a scaling factor that accounts for emission scan count rate. We assume that the emission scan measurements yE and the transmission scan measurements yT are independent Poisson measurements with corresponding means: e-’;(P), For PET reconstruction, one has to do two sets of scans, namely transmission and emission scans. One uses the attenuation correction information obtained from the former scan to aid in estimating the radiotracer emission image from the latter one. Conventional methods of reconstruction are based on linear processing of the transmission and emission data, multiplicative correction of attenuation factors in the sinogram domain followed by FBP to reconstruct the emission image. This approach ignores Poisson nature of the data. Recently, there is growing interest on reconstructheproject methods for attenuation correction in which one reconstructs the attenuation map and, after possibly some processing in Y& A) = TT(bi%(P) r? %(A, p ) ) , (1) the image domain, this map is reprojected to be used in the attenuation correction factors (ACF) computation. The use YP(A,p) = TE(Wi(P)Pi(X) $). (2) of statistical methods for reconstructing attenuation maps as well as emission images is becoming attractive in the medical Here, T~ and rE are transmission and emission scan A C p = l g i j p j and research community, especially due to faster computers and times respectively. l i ( p ) = [Gp]i = gijXj are geometric tomographic projections faster algorithms. In this paper, we reconstruct ACFs using both pi(A) = conventional and penalized-likelihood reconstructheproject of parameters p and A. b,, rT and r? are blank scan, (PL) methods for post-injection transmission scans. For brevity, transmission scan randoms and emission scan randoms count we reconstruct emission images with FBP only. Resolution rates respectively. We assume { b i } , { rT}, { ai}, { r?} and { g i j } matching is critical in attenuation correction, so we add are known constants throughout this work. a post-filtering step to statistical reconstructions to yield approximately Gaussian point spread functions which reduces 11. ACF ESTIMATION artifacts from point spread function mismatches. This post-filter Attenuation correction is a must for quantitatively accurate reduces the negative sidelobes from the point spread function of emission image reconstruction. We define attenuation penalized-likelihood reconstructions [ 11. correction factors (ACFS) -yi(p) = e ’ i ( ~ )= l / c t i ( p ) . This In this work, we study the effects of emission and is the multiplicative factor that corrects for the effects of transmission scan time duration on the variance of the attenuation in the emission data. We consider two different reconstructed emission image for different reconstruction ways of estimating the ACFs: 1) Conventional smoothing zyx zyxw + + + zyxwvutsr 0-7803-5021-9/99/$10.00 0 1999 IEEE. 1842 zyxwvutsrqp zyxwvuts zyxwvut zyxwvutsrq zyxwvutsr zyxwvut zyxw method and 2) Reconstructheproject penalized-likelihood (PL) method. In the non-statistical conventional method, we estimate the emission contamination by: and consequently, i. - smooth{ ki(yF/rE - $)} , (3) and we estimate the ACFs by reciprocating the survival probabilities, that is yi = l/&, where Si = smooth { (yT/rT - r: - ki)/bi} . (4) The smoothing operation is often used to reduce noise in the ACFs. We also use smoothing to reduce noise in the emission contamination estimate in (3). In a statistical reconstruction, one estimates the ACFs by = eli@) where ji is the attenuation map estimate computed by the reconstruction algorithm. The emission contamination estimate (3) is included in the model. The statistical reconstruction is considered in detail in section V. 111. EMISSIONIMAGERECONSTRUCTION For brevity, we consider here the conventional FBP method to reconstruct emission images. We define the attenuated emission projections function as To overcome this problem, we propose an approximation for the probability distribution function of the ACFs. We assume that Ti are lognormal distributed. A random variable is lognormal distributed if its logarithm is normally distributed. We believe this is a very accurate assumption because yi is an estimate of e l t ( f i ) and the projections of any random variable (here li(ji)) can be assumed Gaussian due to the Central Limit Theorem. This provides us extra information about the ACFs. With this assumption, one can compute the mean and variance of Ti’s directly in terms of mean and variance of S; in the conventional method and in terms of mean and variance of li in the statistical method. zyxwvut zyxw zi(X,p) = P i ( X ) a i ( p ) . A linear unbiased estimate of this function is 2. (9) The ACFs yi are not linearly related to variables with known covariances. In the conventional method, .i;. = 1/&. In the statistical method = e’*@). These are both nonlinear functions. Since the covariance of ti can be found exactly for the conventional method and the covariance of ji can be approximated for the statistical method, we can linearize these formulas around &i and & to get an estimate of the covariance of +. This linearization was the method used in [6] to estimate the variances of the ACFs. But, this linearization is not very accurate for especially the conventional method, because the function f(x) = 1/z cannot be closely approximated by a linear function especially when the denominator (survival probabilities) is close to zero and the variance of the denominator is high. - smooth{(yF/TE - T ? ) / E ~ } . Then, an estimate of the projections p i ( A) can be obtained by: pi = yjij. So, for the conventional method, we get: The emission image is reconstructed by standard FBP method. We use the ramp filter only because the estimate is already a smooth estimate of pi(X). Thus, and .i= FBPramp {fiji) . IV. EMISSION COVARIANCE ESTIMATES Even with the lognormality assumption, the covariance matrix of is not easy to compute directly. But, the diagonal of the matrix is known. So, we propose this approximation for the covariances: u~.u~j (6) (12) COV{Ti,+j} M u&;u&jc o v { S i , &j} The covariance of the emission image estimate vector obtained by the above procedure can be written as follows: cov { ;\}= PCOV { p } P’, where the matrix P represents the linear FBP operation with a ramp filter. We need to find the covariance of the random vector 6 = [ij+]L1. The computation of the exact covariance of this expression is computationally intensive and is not desirable. Instead, we prefer to evaluate this covariance as a separable sum of the covariances of the vectors i and +. For this purpose, we consider the Taylor series expansion of &?; in the neighborhood of E i r j where 2 and 7 are mean values of i and -i. respectively. Then: = u+;u+jp(&i,&j), (13) where p ( & , Sj) represents the correlation coefficient of the vector &. In matrix form: c o v {T} M DlCOV {ti} D1, where 1843 zyxwvut zyxwvut zyxwvutsr zyxwvu zyxwvu zyxwvutsrq zyxwvutsrq zyxwvuts We make sure that the diagonal of the covariance matrix of y matches the variances we get from the lognormal assumption. This formula assumes that the correlation coefficient of .i. is largely determined by the smoothing operator B and is the same as the correlation coefficient for &. Plugging in the approximation (9) for .i. and writing Zi 6ismooth{pi(Atrue)}, we get the following COV{@}R DECOV{ i }DE where and R + DTCOV A DT = D E V { smooth{pi(Atrue The mean and covariance of i can be found easily from the expression (5) since it is linearly related to yE. A simple analysis yields: Zj = smooth pi(Atrue)e-~s(~'rUe) and from ( 5 ) and (2): } { 1 {ai} B', COV{ i }= -BV TE A (15) + where qi = (Ejai(ptrue)pi(Atrue) $ ) / E ; , and B is a smoothing convolution matrix along the radial direction of the projection space. It was suggested in [5] that angular smoothing is not desirable in attenuation correction, so we smooth only in radial direction. For conventional ACF computation, ignoring the noise in the emission contamination estimate, the covariance of ti can be found from (4) and (1). 1 COV{&} = -BV {si} BT, TT (16) where si = ( b i a i ( p t r u e ) + TT + & ) / b f . Here, B is the same smoothing matrix as in (15). The same operator B is used to obtain both i and ti to avoid artifacts from resolution mismatch [7, 81. We used Gaussian smoothing as suggested in [7] which avoids any artifacts in the reconstructed image. The mean of emission contamination can be determined from (3) as N ii = B [ I C i ~ i a i ( ~ ~ ~ ~ ~ ) pThe ; ( variance X ~ ~ ~ of ~ )&i ] ~ can = ~ be. found from (16) as A U$i = Si/TT B&. k Using (4), one can find the mean values of & as 6=B[(~i(p~~~~)]zl. (17) The variance of the sum over a region of interest in the emission image can be found from (6), (14), (15) and (16) as V. PENALIZED-LIKELIHOOD ATTENUATION RECONSTRUCTION While conventional method of ACF computation has been used for some time, reconstructheproject methods have gained some interest recently. In a statistical reconstructheproject method for ACF computation, an attenuation map estimate ji is found from noisy transmission data by maximizing the penalized-likelihood objective function @ ( p ;yT) = L ( p ; yT) P R ( p ) ,where L ( p , yT) is the log-likelihood function and R ( p ) is a regularizing roughness penalty function. After estimating the attenuation map ji, we estimate the ACFs by: 9; = e l i ( f i ) , where l i ( j i ) = [GfiIi is the geometric projection of the attenuation map estimate fi. If one uses FBP for emission reconstruction, then i should be smoothed to yield similar resolution with the -i. [9] in order to reduce resolution mismatch artifacts. A. Resolution Penalized likelihood (PL) or penalized weighted least squares (PWLS) methods are very attractive image reconstruction methods due to their superb noise reduction properties. The variance weighting in PWLS method reduces the variance of the estimates as compared to penalized unweighted least squares (PULS) or FBP reconstructions, because it makes use of the statistical information in the measurements. However, attenuation maps reconstructed with PL or PWLS methods have non-uniform resolution [1] even with a quadratic penalty. This non-uniform resolution is caused by the variance weighting in PWLS (or PL) method and hence does not exist in a PULS reconstruction. Due to this non-uniform resolution, ACF computation by PL method from a real transmission scan causes resolution mismatch between the emission data and reconstructed ACFs. This mismatch reveals itself as artifacts in the final reconstructed emission image. Fessler's certainty based penalty [9] yields more uniform resolution in terms of the average FWHM of the point spread function over the image. But, it still has non-uniform resolution in that the psf is not circularly symmetric but the level contours look like ellipses whose orientation are image dependent and space-variant. Stayman and Fessler have recently proposed a new modification to the quadratic penalty [lo] which yields more circularly symmetric uniform resolution properties. We used this modification in our reconstructions. This modification makes the resolution properties of the PL method close to PULS method. Quadratic PULS method was shown to be essentially equivalent to FBP method with the following constrained leastsquares (CLS) filter defined in spatial frequency domain by (equation (50)in [9]) sinc(ku) / sinc(u) Fp(u;P) N with vE = qi(wF)2and vT = Ezl S ~ ( W ? ) ~and , where U is a vector of ones in the region of interest and zeros elsewhere. We define the vectors WE A = B ' D ~ P ' ~wT , B'DTP'U. = sinc(k.u)2 + +3 zyxw ' U E [0,0.5] (19) where U denotes spatial frequency, k is the ratio of the detector strip width to the pixel size of the system model, and c is a constant dependent on system geometry. This CLS filter has high negative sidelobes in the space domain. The filters that 1044 zyxwvuts zyxwvuts zyxwvutsrq zyxwvutsr zyxwvuts zyxwvuts zyxwvuts zy zyxwvutsrqpo smooth the ACFs and emission data have to be matched. So, the emission data should be blurred with the same filter (19). But, due to high negative sidelobes of filter in (19), after dividing the appropriately blurred emission data to computed survival probabilities from reconstructions, we get artifacts especially for higher blurring amounts (higher Ps) around the boundaries of the image. So, we conclude that the results in [7] does only hold for Gaussian smoothing. estimate of the attenuation map p. We again ignore the noise in the emission contamination estimate and use the mean value for it in our approximations. The formula yields: where To overcome this problem, we first reconstruct a higher resolution image using a smaller ,f3 value than desired and then we filter the projections with the'following filter: Here R is the Hessian of the penalty function and includes the modified penalty weights [ 101 and FT = rT E * . + where Fg(u; w) is the desired Gaussian filter with desired FWHM w. Now, the emission data is also filtered with the Gaussian shaped filter Fg(u;w). This approach reduces artifacts and yields acceptable images. The ACF computation in this case is done as follows: In this case, the variance of the sum over a region can be predicted with a formula similar to (18). The emission part of the formula is now wE = B'V eli P ' u and qi remains same. The transmission part changes a lot due to statistical method as term wT should be changed to: { -} ji = argmax@(,u;yT), P f = BzGji, eii, (20) =Ui = where B2 is the convolution matrix corresponding to above. F2( U ) { -} B. Covariance Approximations The covariance formula in (9) is still valid in PL transmission reconstruction. We use the following first order Taylor series expansion for the ACFs: -.. - ;U = el; M eft + e'*(ii - li), (21) r where = BzGD is the mean projection vector where G , = arg max@(p;jjT) is the image reconstructed with noiseless P20 data. ji is a very good approximation for the mean of fi [4]. Note that, we do not use the lognormality assumption here, because we believe that the above approximation is accurate enough and lognormal assumption leads to much more computation. From (21) and (20), zyx The most computationally intensive part' in this computation is the part where H-lv* should be computed for U* = G'B',V Eie" P ' u . This operation can be performed by solving the equation Ha: = v* using iterative methods such as conjugate gradient. Also, we assume the mean for =U is now, ? = e 1;. These variance predictions are useful, because they do not require hundreds of empirical reconstructions of data [4]. However, they require knowing the true parameters and noiseless sinograms. For real data, these are not known, but one can still get a good approximation of variances by replacing the true parameters by their noisy counterparts [4]. Finally, the optimal time fraction for the emission scan can be found by minimizing the variance in (18) with respect to the emission scan time when total scan time is fixed. For the QPL method, the simple analysis yields E Topt To find the covariance of the implicitly defined estimator ji, we use the formulas introduced in [4]. The general form of penalized-likelihood estimates is ,u is the parameter vector fi = arg max@(,u;yT), where - - totalvE VE-UT . Note that for the conventional method, the above formula is invalid because the vT term is not independent from the scan time duration rT. J ! and yT is the measurement vector. This defines an implicit function fi = h(yT). A first order Taylor expansion of the equation V @ ( p ;yT) = 0 around (a,yT) yields the following approximation [4]: VI. RESULTS We have done a series of simulations to test the proposed variance predictions and to find the optimal scan times under a total scan time constraint. We used two 2-D images (22) corresponding to attenuation map and emission image to Cov {ji} M QCov { yT} Q', generate noisy transmission and emission data with 150000 and where Q = [-VZo@(/i,GT)]-' V1l@(p,yT). We use this 50000 counts per minute respectively. The true images are formula to evaluate the covariance of the penalized-likelihood shown in Figure 1. The transmission scan had 5% randoms 1845 zyxwvutsr sum of pixels inside the heart region zyxwvut and an emission contamination of 5%. Emission scan had 10% randoms. Randoms rates were assumed to be constant. The total scan time was 20 minutes. To obtain empirical standard deviations, 300 realizations were generated for each scan time distribution. The emission images were reconstructed with FBP with a smoothing filter that yields about 9 mm FWHM resolution in the image domain. ACFs were computed using the conventional and quadratic penalized-likelihood statistical methods. The resolutions for the ACFs were matched for these two methods. The standard deviations of the sums over the heart region in the reconstructed emission images were found empirically and predicted analytically using the derived formulas. The results are shown in Figure 2 as a plot of standard deviation estimates versus emission scan time fraction. The statistical method not only reduces the overall variance, but also yields a larger optimum emission scan time fraction (about 40%) as compared to the conventional method (about 30%). The standard deviation is reduced by about 1520% in the statistical method as compared to the conventional method. The predictions seem to match the empirical data for the statistical reconstruction, but the predictions for the conventional method seem to underestimate the standard deviations. We conjecture that, the approximations used in deriving the variance formulas causes the mismatch. We are currently working on improving our approximations. Due to highly nonlinear processing of data however, it is likely that there will be some discrepancy between predicted and empirical standard deviation estimates. zyxwvutsrqp zyxwv I 1 wnlulon s u n tlnn hction Figure 2: Standarddeviation of the sum over the heart region estimates versus emission scan time fraction for conventional and statistical ACF computations. VIII. REFERENCES J. A. Fessler and W. L. Rogers, “Spatial resolution properties of penalized-likelihood image reconstruction methods: Spaceinvarianttomographs,”IEEE Tr.Im. Proc., vol. 5 , no. 9, pp. 134658, September 1996. T. Beyer, P. E. Kinahan, and D. W. Townsend, “Optimization of transmission and emission scan duration in 3D whole-body PET,” IEEE Tr. NUC.Sci., vol. 44,no. 6, pp. 2400-7, December 1997. C. Comtat, P. E. Kinahan, M. Defrise, C. Michel, and D. W. Townsend, “Fast reconstruction of 3D PET data with accurate statistical modeling,” IEEE Tr.NUC.Sci., vol. 45, no. 3, pp. 10839, June 1998. J. A. Fessler, “Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): Applications to tomography,” IEEE Tr. Im. Proc., vol. 5 , no. 3, pp. 493-506, March 1996. M. E. Daube-Witherspoon and R. E. Carson, “Investigation of angular smoothing of PET data,” in Proc. IEEE NUC.Sci. Symp. Med. Im. Con$, volume 3, pp. 1557-61,1996. C.Steams and D. C.Wack, “A noise equivalent counts approach to transmission imaging and source design,” IEEE Tr. Med. Im., vol. 12, no. 2, pp. 287-292, June 1993. A. Chatziioannou and M. Dahlbom, “Detailed investigation of transmission and emission data smoothing protocols and their effects on emission images,” IEEE TI:NUC.Sci., vol. 43, no. 1, pp. 2904, February 1996. H.Erdogan and J. A. Fessler, “Statistical image reconstruction methods for simultaneous emissiodtransmission PET scans,” in Pmc. IEEE NUC.Sci. Symp. Med. Im. Con$, volume 3, pp. 157983, 1996. J. A. Fessler, “Resolution properties of regularized image reconstruction methods,” Technical Report 297, Comm. and Sign. Proc. Lab., Dept. of EECS, Univ. of Michigan, Ann Arbor, MI, 48109-2122,August 1995. W. Stayman and J. A. Fessler, “Spatially-variant roughness penalty design for uniform resolution in penalized-likelihood image reconstruction,” in Proc. IEEE Intl. Con$ on Image Processing, volume 2, pp. 685-9,1998. VII. CONCLUSION zyx We presented new approximate formulas for covariances of reconstructed emission images with conventional and statistical ACF computation for post-injection scans. These formulas can be used to predict the variance of the sum over a region of interest in the final reconstructed emission image instead of expensive empirical reconstructions. These formulas can also be used to determine optimal scan times devoted to emission and transmission scans under a total scan time constraint. Results show that, statistical ACF computation not only reduces the overall standard deviation but also yields higher optimum emission scan time fraction than the conventional method. The covariance approximations for statistical method work well for quadratic penalties, but not for non-quadratic penalties. We plan to extend our work to cover non-quadratic penalties in the future. attenuation map emission Image zyxwvutsr Figure 1: The attenuation map and emission image used to generate simulation data. This work was supported in part by NIH grants CA-607 11 and CA-54362, and by the Whitaker Foundation. First author was supported by TUBITAK-NATO Science Fellowship for doctoral studies. 1846