MEDICAL IMAGING EXPERIMENT医学成像实验.doc

MEDICAL IMAGING EXPERIMENTPART I0. Introduction the goals of this experimentThe goal of this experiment is to introduce the field of medical imaging and to implement basic techniques of iterative reconstruction. After doing this experiment, you should have an idea about inverse problems arising in medical tomography and about numerical optimization methods that can be used in order to solve these problems.The first part of the experiment is dedicated to introduction into tomography (particularly, the notion of projection and Radon transform), and basics of optimization a theory which tells how to algorithmically minimize a given convex function of many variables. We will implement a simple iterative algorithm for reconstruction of an object from its projection. The second part is dedicated to the relation between the Radon and the Fourier transforms, and the frequency-domain reconstruction techniques arising from it. We will implement an iterative inverse non-uniform fast Fourier transform (NUFFT) and use it for tomographic reconstruction. 00. General remarksThe main bibliographical source is the book by Kak and Slaney, Principles of computerized tomography, available online on . Chapters 1 and 3 referenced in this experiment are attached.The programming language for this experiment is MATLAB. - theoretical questions requiring a written answer.: - computer questions requiring writing a code in MATLAB.2 - computer questions requiring to show a figure or a plot.1. Theoretical background tomography and Radon transform1. Read Chapter 1 (pp. 1-4) from Kak&Slaney. What is the main difference between computerized tomography (CT) and conventional X-ray imaging ? What is the reason for the word “computerized” in the term CT ?2. Read Chapter 3, part 3.1 (pp. 1-8) from Kak&Slaney. We will deal with parallel projections only. What is a projection ? Give a mathematical expression for a projection at angle q of an object f(x,y) What is the Radon transform ? Find the Radon transform of . Sometimes the collection of projections of an object is called a sinogram. What is, in your opinion, the reason for this name ? Explain why it is necessary to take projections in the range of angles q between 0 and p and not 2p.3. In practical applications, we will be interested in a discrete version of the Radon transform. The object f(x,y) will be represented as an image, i.e. a matrix Fmn, where m,n=1, N are pixel coordinates (w.l.o.g., the matrix is assumed to be square). Projections will be taken at K discrete evenly spaced angles qk between 0 and p; a projection at angle qk will denoted by a vector Plk (l=1,L , k=1,K) . The Radon transform will be therefore a discrete linear operator . Write an expression for a discrete approximation of the Radon transform. Approximate integrals as Riemann sums. Denote the step by and . Denote the dimensions of the matrix Fmn by . Draw a schematic diagram of a discrete projection. What is the maximum number of pixels (elements of Fmn) that are summed along a line of every projection ? What should be the value of L in terms of N ? Explain. Hint: look at projection at angle 450.4. Now assume that the image is parsed into a column-stack, i.e. represented as a column vector containing successively concatenated columns of the matrix Fmn. We will denote the obtained vector by.The Radon transform is a linear mapping which can be represented by an matrix R. What can you say about the nature of the matrix R ? Does it contain many non-zero elements ? What are the dimensions of the matrix R for a image, and projections taken every 10 ? What problems would you expect in treating a matrix of such dimensions ? How does the nature of the matrix is helpful in this regard ? Assume that you are given a black box, which, given an image x (as a column stack), is able to compute its Radon transform, i.e. y=Rx. Suggest a way to find the elements of R. Hint: use the fact that in every finite-dimensional linear space L there exists the standard basis , such that.Think how to express k-th column of R.5. The problem of tomographic reconstruction can be expressed as the problem of solving the set of linear equations, given in the matrix form as,where x is the unknown object (in column-stack representation) and y are the measured projection (also column stack). Problems of such kind are termed as inverse problems. Mathematically speaking, these problems are always ill-posed. Write an analytical solution of the inverse problem. Pay attention: R is not square ! What is the obvious difficulty in the proposed solution ? Next section will be dedicated to an attempt to overcome this difficulty.2. Optimization theory1. We will give a brief introduction to optimization theory. This theory will allow create iterative algorithms to produce a solution of the inverse problem. The main result in optimization is the following: Let be a convex C2 function, i.e. and D be a convex set, i.e.Then, possesses a global minimum, which coincides with its local minimum. Which of the following is a convex function ? Explain by drawing. Mark the local and the global minima of the functions.xa.xb.xc.xd.2. The meaning of this result is that in order to find the global minimum of a convex function, it is enough to be able to find its local minimum. We will define the gradient of the function as the vector of partial derivatives:and the Hessian as the matrix of second-order derivatives: Assuming that is a C2 function (i.e. has continuous derivatives of 2nd order), what type of a matrix is the Hessian ? Let be the minimizer of . Recall the results from Calculus and write what should the gradient of the function in the point obey. 3. We will be interested in solving linear equations of the type .The inverse problem can be formulated as an optimization problem: given the observations y, we would like to find the object x, such that the difference between y and Rx are minimal. I.e. (*) Write the function in (*) explicitly as a quadratic function, i.e. Prove that is convex. Assume that is a positive semidefinite matrix, i.e.Bonus: prove that is positive semidefinite. Find the gradient . Bonus: find the Hessian . Using the result of paragraph 2 (condition on the gradient at the minimum), find an analytical solution to the inverse problem. Is it equal to the solution in section 1 ?3. Part I lab work - simulations1. Assume images of size 3232, angles 0:2:179. : Use the result of 1.4. to find the Radon transform matrix R. Instructions: since dealing with a very large matrix, you will not be able to store it “as is”. We will take advantage of the matrix sparsity. a. Compute the dimensions and of the matrix R: run MATLAB function radon on an arbitrary image of size 3232 to see the dimensions of the outcome. b. Estimate the number of nonzero elements in the matrix. c. Allocate a sparse matrix with this number of nonzero elements (advised to take a bit larger number) using spalloc.d. Create the matrix according to 1.4.2 Use spy to display the pattern of the matrix R.: Generate a simulated brain image of size using phantom (use defaults). Compute the Radon transform by using the function radon and by parsing the image into a vector and computing Rx. Compare the results. 2. Recall that we are interested to minimize the cost function in order to reconstruct the image from the observation y.: Write a MATLAB function f,df,d2f = func(x,) which will compute the values of at point x. Use nargin to allow different number of output arguments. Use the results from 2.3.: Reconstruct the image from the data y supplied in dataset1. Use MATLAB function fminunc for unconstrained minimization of in the following way:OPT = optimset(Display,iter,DerivativeCheck,on, MaxIter,100,Gradobj,on,Hessian,off,LargeScale,on);xo = fminunc(func,x0,OPT, ); (See HELP).: Compare to reconstruction using iradon (inverse Radon transform using a frequency method). 2 Display the results.3. Now we will get some impression about what is termed ill-posedness. In the above setting (images of size 3232, angles 0:2:179), is the problem over-determined or under-determined ? (i.e. are there more variables than data or vice versa ?).2 Plot the singular values of R on logarithmic scale (use svd(full(R). This may take a while). Do you expect a problem to solve the system Rx = y as you suggested in 1.5 ? Explain. Now assume angles 0:10:179, images of size 3232. Construct the matrix R for this case. Here we have 5 times less projections than previously. Is now the problem over-determined or under-determined now ?.2 Plot the singular values of R on logarithmic scale for this case. What is the difference ? Do you expect now a problem to solve the system Rx = y as you suggested in 1.5 ? Explain. In which of the cases (the first or the second) do you think the problem is ill-posed ?: Add Gaussian noise with variance 0.01 to the data y in dataset2 (use randn). Reconstruct the image. What do you observe ? Explain in light of the previous conclusions. Give a definition to ill-posedness.4. One of the ways to cope with ill-posedness is called Tikhonov regularization. In this method, the original problem is replaced by where is some parameter that is chosen depending on the problem and the data. From the first order optimality condition (), write an analytical solution for the minimizer of (i.e. the solution of the new optimization problem). Compare to the result in 1.5 (or 2.3.). Explain in general terms how the regularization solves the problem you noticed in 3.3.: Write a MATLAB function f,df,d2f = func_r(x,) which will compute the values of at point x. : Perform reconstruction of the noisy data in 3.3 (dataset2) with regularization. Select . 2 Show the reconstructed images obtained with different values of . Explain the results. What is the problem with the proposed regularization ?PART II1. Theoretical background Fourier slice theorem1. Read Chapter 3, part 3.2 (pp. 56-60) from Kak&Slaney. We will deal with parallel projections only.What is the relation between the Radon transform, the 1D Fourier transform and the 2D Fourier transform ? Given the Radon transform , suggest a frequency-domain algorithm for reconstruction of .2. Gridding1. We will derive a nave and simple method called gridding for fast computation of non-uniform Fourier transforms (NUFT). For simplicity, lets consider the 1D case. Assume we are given a 1D discrete signal . Its discrete time Fourier transform (DTFT) is a periodic function of continuous variable, given by.A standard (uniform) Fourier transform can be thought of as sampling on evenly-spaced grid . A NUFT is obtained by arbitrary sampling at frequencies (not necessarily uniform). We will denote the values of the NUFT by , and will write in matrix-vector notation:. Assume that we are able to compute the fast Fourier transform (FFT) on the grid (q is called the oversampling factor). Show how to approximately compute the forward NUFT using the oversampled FFT and interpolation. Estimate the computational complexity of the algorithm. Show how to approximately compute the inverse NUFT using the oversampled FFT and interpolation. Estimate the computational complexity of the algorithm.2. Now assume the signal f is 2D. Generalize the previous notation and results to the 2D case.3. The Fessler-Sutton optimal NUFFT1. The problem of approximate forward NUFT computation is the choice of interpolation coefficients. We will derive here optimal interpolation in min-max sense as proposed by Fessler and Sutton. Again for simplicity, we consider the 1D case.Let us denote the q-oversampled FFT by . The forward NUFFT is obtained by interpolating the values of to the non-uniform grid: ,where denotes interpolation operation, which makes use of p neighboring uniform samples (i.e. each row of has only p non-zero elements). Fessler and Sutton proposed finding by minimization of the worst-case error at each point of the non-uniform grid. In other words, the interpolation coefficients for point k of the non-uniform grid (the vector of non-zero elements of the k-th row of ) are obtained as the solution to the following minimization problem:where denotes the p rows of used for the computation of k-th point of the non-uniform grid (corresponding to the interpolation coefficients). Find an analytic expression for . Hint: the constraint can be replaced by . Substitute the expression you found and express . Assume that the interpolation coefficients are pre-computed. Estimate the complexity of the Fessler-Sutton forward NUFFT algorithm.2. Suggest how to perform iterative computation of the inverse NUFFT using the forward NUFFT. Write the cost function that should be minimized and its gradient. Estimate the computation complexity of the cost function and the gradient computation. 4. Part I