Positron Emission Tomography Image Reconstruction (Section 1.7.5 of Linear and Nonlinear Programming by Griva-Nash-Sofer, SIAM, 2009) Positron emission tomography (PET) is a medical imaging technique that helps diagnose disease and assess the effect of treatment. Unlike other imaging techniques such as X-rays or CT-scans that directly study the anatomical structure of an organ, PET studies the physiology (blood flow or level of metabolism) of the organ. Metabolic activity is an important tool in diagnosis: cancerous cells have high metabolism or high activity, while tumor cells damaged by irradiation have low metabolism or low activity. Alzheimer's disease is indicated by regions of reduced activity in the brain, and coronary tissue damage is indicated by regions of reduced activity in the heart. In a PET scan the patient is injected with a radioactively labeled compound (most commonly glucose, but sometimes water or ammonia) that is selected for its tendency to be absorbed in the organ of interest. Once the compound settles, it starts emitting radioactive emissions that are counted by the PET scanner. The level of emissions is proportional to the amount of drug absorbed, or in turn, to the level of cell activity. The emissions are counted using a PET scanner that surrounds the body. Based on the emissions counts obtained in the scanner, the goal is to determine the level of emissions from within the organ, and hence the level of metabolic activity. The output of the reconstruction is typically presented in a color image that reflects the different activity levels in the organ. We describe the physics of PET in further detail. As the radioisotope decays, it emits positrons. Each positron annihilates with an electron, and produces two photons which move in nearly opposite directions, each hitting a tiny photo-detector within the scanner at almost the same time. Any nearsimultaneous detection of an event by two such detectors defines a coincidence event along a coincidence line. The number of coincidence events yj detected along each of the possible coincidence lines j is the input to the image reconstruction. Consider the situation depicted in the figure below, where a grid of boxes or voxels has been imposed over the emitting object (for simplicity, the figure is depicted in two dimensions; the concept is readily extended to three dimensions). Given a set of measurements yj along the coincidence lines j = 1,...,N we seek to estimate xi, i = 1, ...,n the expected number of counts emitted from voxel i, where n is the number of voxels in the grid. PET Most reconstruction methods are based on a technique known as filtered back projection. Although this technique yields fast reconstructions, the quality of the image can be poor in situations where the amount of radioactive substance used must be small. Under such situations it is necessary to use a statistical model of the emission process to determine the most likely image that fits the data. The approach is via the maximum likelihood estimation technique. The radioactive emissions from voxels i = 1,... ,n are assumed to be statistically independent random variables that follow a Poisson distribution with mean Xi. Denote by Ci,j the probability that an emission emanating from voxel i will hit detector pair (coincidence line) j. The n x N matrix C = {Ci,j}depends on the geometry of the scanner and on the tissue being scanned, and is assumed to be known. Using these assumptions one can show that the emissions emanating from voxel i and hitting detector pair j are also independent Poisson variables with mean rate Ci , j xi and the total emissions received by the detector pairs j = 1, ... , N are independent Poisson distributed variables with mean rate C x . i, j i i Let q = CeN where eN is a vector of 1's. The vector q denotes the sum of the columns of C (which need not be 1). It is computationally easier if we write the optimization model using the logarithm of the likelihood function. If we ignore a constant term, the resulting logarithm is f ML = qT x + y j log ( CT x ) j j (See the Exercises.) Since the emission level is nonnegative, the final reconstruction problem becomes qT x + y j log ( CT x ) max f ML = j j s.t. x 0 The size of the problem can be enormous. If one wishes to reconstruct, say, a volume of, say, 5 cubic cm at a resolution of half a millimeter, then the size of the grid would be 100 by 100 by 100, corresponding to n = 100,000 variables. Problems of this size (or even larger) are not uncommon. The size of the data is also huge. The scanner may have thousands of photo-detectors and since any pair of these can define a coincidence line, the number of coincidence lines N can be on the order of millions. Since every function evaluation requires the computation of a matrix product CTx, and the matrix C is large, the function evaluations are time consuming. The efficient solution of such large problems often requires understanding of their structure. By structure we mean special characteristics of the function, its gradient, and Hessian. Often structure is associated with the sparsity pattern of the Hessian, that is the number of zeros, and possibly their location. The special structure of fML and its derivatives can be used in designing effective methods for solving the problem. Here we will just give the formulas for the derivatives. Defining y = CT x , we can write the gradient and Hessian of the objective function, respectively, as 1 y f ML ( x ) =q + CY 2CT 2 f (x ) = CYY ML diag(y ) . The matrix C itself is sparse, and only a small fraction of its where Y diag( = = y ) and Y are of course also sparse. Even so, the Hessian entries are nonzero. The diagonal matrices Y and Y 2 f ML ( x ) is dense; almost all of its entries are likely to be nonzero. A key challenge in the design of effective algorithms is to exploit the sparsity of C. Tasks to be accomplished in the case study a) Show that the likelihood that the total emissions received by all detector pairs is y given that the expected number of emissions from all voxels is x is given by: yj Ci , j xi i e y Ci , j xi e j ( y j ) i = , where y j = = P( y | x ) y y ! ! j j j j yj C x i, j i i (Use hints given at the end of this document.) Prove the final expression for the maximum likelihood estimation objective function fM L Hint: Take the logarithm of the likelihood and omit the constant term that does not depend on x. b) Derive the formulas for the gradient and Hessian matrix of fML. c) Show that the Hessian of fML may be dense, even when its matrix factors are sparse. Do this by considering C = (I I en) and y= y= e 2 n+1 where I is the identity matrix, and ek is a vector of ones of size k. Show that every element of the Hessian is nonzero. d) The purpose of this problem is to write a program in MATLAB to solve a PET image reconstruction problem. Choose the algorithm that you think will be most efficient to solve this case, based on the information that you have up to this point. (i) Write a MATLAB script of your chosen solution method and test it on a problem with n = 9 variables corresponding to a 3 x 3 grid, and with N = 33 detector pairs. The data are C= (B B B), where B is a sparse nx(n + 2) matrix with the following nonzero entries: Bi,i =a, Bi,i+1 = b, Bi,i+2 = a, i = 1, ... , n, where a = 0.18, b = 0.017, and yT =( 0 0 1 19 27 30 40 50 35 15 1 0 0 1 7 20 38 56 55 38 20 7 1 0 1 3 17 38 40 20 7 1 0). (ii) Test your software on a problem with n = 1080 variables corresponding to a 36 x 30 grid, and with N = 2164 detector pairs. The data are C = (B 2B), where B is defined as in part (i) with the parameter values a = 0.15 and b = 0.05. The vector y can be downloaded in text format from the Web page for this book (http://www.siam.org/books/ot108). Display the values of the first row of the reconstructed image. (iii) Identify the image you obtained in (ii). You will need software for displaying intensity images. Hints for part (a): Require some basic background in stochastic methods First let Zij be the number of events emitted from voxel i and detected at coincidence line j, and let Yj be the total emissions received by detector pair j, for j = 1, ... , N. Use the assumptions given in the main section to prove that {Zij} are independent Poisson variables with mean Ci , j xi and that {Yj } are independent Poisson distributed variables with mean rates y j = C likelihood function may be written as yj Ci , j xi i e Ci , j xi i P( y | x ) = ! y j j x . Then show that the i, j