\documentclass[11pt]{article} \usepackage{palatino,graphicx} \pagestyle{headings} \def\hdr#1 {{\parindent 0pt \vskip -0.7in% EE369C Fall 2011-12 \hfill \\ Medical Image Reconstruction \hfill \\ \vskip 0.1in }} \topmargin=-0.5in \textheight=9in \evensidemargin 0.0in \oddsidemargin 0.0in \setlength{\textwidth}{6.5in} \begin{document} \ \hdr{1} \begin{center} {\Large\em EE369C: Assignment 6} \\[0.15in] {\em Due Thursday, Nov 17} \end{center} Your task for the homework this week is write a short description of the final project you would like to do. Submit a short description of the problem, what approaches you will look at, and what resources you need (references, data, etc). The entire description should only be a paragraph or so. Submit it the same was as you do your assignments on EE class. You can choose any project that is related to the subject of the class. This can be something related to your Ph.D. work, some problem you found interesting, or one of the problems listed below. The scope of the problem should be something beyond what you have done in the homework. However, it should also be restricted in scope, to fit the limited time available. One good choice is the implementation of one of the algorithms that we have discussed, that haven't appeared in the homework. The projects will be assessed based on the difficulty and importance of the problem, and on how well you do with your problem. You can do well either by doing a thorough job on a well understood problem, or by making progress on a difficult or new problem. The projects are due Dec 16th, which is the last day of finals week. Present your project by writing a 10-15 page report of your project. This is similar in length and scope to the assignment write-ups you have been doing. \section*{Possible Project Topics} \subsection*{Non-Cartesian Reconstruction} \paragraph*{1. Inverse Gridding} Given a 2DFT MR image, use inverse gridding to generate simulated data sets for other k-space trajectories, such as the rt-spiral and var-dens trajectories. Choose an appropriate oversampling ratio and gridding kernel, and show that the reconstruction of the simulated data looks like you would expect. This is a useful tool for several of the projects below, and is a good place to start. \paragraph*{2. Density Estimation\ \ } We can represent the gridding operation as a matrix multiplication \[ {\bf m} = G {\bf d} \] where ${\bf m}$ is the vectorized k-space data resampled on a Cartesian grid, and ${\bf d}$ is the acquired data. The matrix $G$ is large ($n^2\times m$, for an $n\times n$ image and $m$ data samples). However it is sparse. If the kernel is 4 samples wide, then only 16 samples in each column are non-zero. This means that using a matlab sparse matrix is very efficient. If you preallocate the sparse matrix to be large enough, this can be very fast. We can compute an improved density function by solving the problem of what data vector {\bf d} produces a uniform {\bf M} over the region of support. One solution to this problem is to use the iterative lsqr() matlab routine which uses a conjugate gradient algorithm to solve the least squares minimization problem. Use this approach to compute improved density estimates for the rt-spiral data set, and the var-dens data set. Note that the density will be a function of the particular kernel you choose. Also, you may want to introduce regularization. \paragraph*{3. Least Squares Reconstruction} In class we talked about how the gridding reconstruction is essentially a weighted least squares solution to the reconstruction problem. Show that directly solving the least squares problem is also now tractable using iterative solvers. This is one component of the non-Cartesian SENSE reconstruction in part 8, and relies on the inverse gridding operation of part 1. \subsection*{Off-Resonance Correction} For these projects you can use the data from Assignment 3, and also the high-resolution spiral coronary data set \verb+3T_Coronary.mat+ available in the \verb+data+ directory of the web site. A description of the data set is in \verb+3T Coronary.txt+. \paragraph*{4. Autofocus Metrics} In the homework assignment we minimized the magnitude of the imaginary component of the image as a focus metric. In the original Noll {\em et al.} paper, \begin{quote} ``Deblurring for non-2D Fourier transform magnetic resonance imaging.'' Noll DC, Pauly JM, Meyer CH, Nishimura DG, Macovski A. {\em Magn Reson Med.} 1992 Jun;25(2):319-33. PMID: 1614315 \end{quote} they used a fractional power of the magnitude of the imaginary component \[ f(x,y) = \sum_A |m(x,y)|^\alpha \] where $\alpha <1$, and $A$ is a local neighborhood of $x,y$. What is the effect of choosing different values of $\alpha$? What value should you choose? \paragraph*{5. Fast Multifrequency Reconstruction} In the homework assignment we simply used a nearest neighbor interpolation for the multifrequency reconstruction. Clearly we can do much better. This would allow for fewer reconstructions for a given reconstruction fidelity. Three papers that describe alternatives are \begin{quote} ``Fast conjugate phase image reconstruction based on a Chebyshev approximation to correct for B0 field inhomogeneity and concomitant gradients.'' Chen W, Sica CT, Meyer CH. {\em Magn Reson Med.} 2008 Nov;60(5):1104-11. PMID: 18956462 ''Multifrequency interpolation for fast off-resonance correction.'' Man LC, Pauly JM, Macovski A. {\em Magn Reson Med.} 1997 May;37(5):785-92. PMID: 9126954 ``Off-resonance correction of MR images.'' Schomberg H. {\em IEEE Trans Med Imaging.} 1999 Jun;18(6):481-95. PMID: 10463127 \end{quote} Implement one or more of these, and evaluate how much time can be saved in the reconstruction of the phantom data set of Assignment 3. \paragraph*{6. Joint Estimation of the Frequency Map and Image} A number of approaches have been proposed to jointly solve for the frequency map and the image reconstruction. A few of these are \begin{quote} ``Algebraic reconstruction for magnetic resonance imaging under B0 inhomogeneity,'' YM Kadah and X Hu, {\em IEEE TMI}, vol. 17,no. 3, pp. 362Ð370, 1998. %``Iterative reconstruction of single-shot spiral MRI with off-resonance.'' TB Harshbarger and DB Twieg, {\em IEEE TMI}, vol. 18, no. %3, pp. 196-205, 1999. ``Fast, iterative image reconstruction for MRI in the presence of field inhomogeneities.'' BP Sutton, DC Noll, and JA Fessler, {\em IEEE TMI}, vol. 22, no. 2, pp. 178-188, 2003 \end{quote} These work by estimating what object corresponds to the data that was acquired. This has advantages when there are sharp discontinuities in frequency that can't be resolved with a multifrequency reconstruction. Implement one of these, and reconstruct the data set from Assignment 3, and the 3T Coronary data set. % \item {\bf\em Fan Beam CT\ \ } %Next week we will be talking about backprojection and CT. In practice, almost all CT data is acquired using a fan beam. We have some fan beam data from a current commercial whole-body CT system. For this project, you will reconstruct the actual fan beam data. There are several alternative fan beam reconstruction algorithms. The filtered backprojection case over a 360 degrees of rotation will be covered in the homework. Some project options are: %\begin{enumerate} % \item Reconstruct the fan-beam data using rebinning, which we won't cover in the homework. % \item Modify the fan beam filtered backprojection to use only 180 degrees plus the fan angle. This is known as "halfscan" reconstruction, and is covered in the Kak and Slaney book. % \end{enumerate} \subsection*{Parallel Imaging} \paragraph*{7. Calibration Area/Kernel Size in SPIRiT} In Assignment 5 we used a calibration area of 24, and a kernel size of 5. What is the effect of different kernel sizes and calibration areas? For the data set we used, is there a better combination of the two? In general, how should we choose the calibration area and the kernel size? An article that looks at these issues for GRAPPA is \begin{quote} ``On optimality of parallel MRI reconstruction in k-space.'' Samsonov AA. {\em Magn Reson Med.} 2008 Jan;59(1):156-64. PMID: 18058935 \end{quote} This may give you some ideas what to look for. \paragraph{8. Non-Cartesian SENSE or SPIRiT} There are a number of ways to do parallel imaging with non-Cartesian data sets. The first method is described in the paper by Klaas Pr\"{u}ssmann, \begin{quote} ``Advances in sensitivity encoding with arbitrary k-space trajectories.'' Pr\"{u}ssmann KP, Weiger M, B\"{o}rnert P, B\"{o}siger P. {\em Magn Reson Med.} 2001 Oct;46(4):638-51. PMID: 11590639 \end{quote} This was the first well known paper to use a conjugate gradient algorithm for parallel MR image reconstruction, and has a nice, simple description of how to implement it. This is a good place to start. Another method is SPIRiT, as described in the Lustig paper you have already studied. This is a little more difficult to implement because it solves for all of the individual coil images, and that makes some of the operators more complex. However, being autocalibrating, it is preferred. A variable-density spiral four coil data set is available on the web site, \verb+spiral_sense.mat+. A description of the data is in \verb+spiral_sense.txt+. This is a simulation, so we have the actual sensitivities if you would like to use them. Implement either of these approaches. For both of them you will need the adjoint of the gridding operator, which is inverse gridding, so you will need to implement that first. The data set has 21 spiral interleaves, and four coils. Try reconstructions with acceleration factors of 2, 3, and 4 by subsampling in the interleaf dimension. \paragraph*{9. Null Space or Eigenvector Approaches to Parallel MRI} The SPIRiT algorithm is the first of a broader class of parallel imaging algorithms known as null space or eigenvector approaches. The basic idea of SPIRiT is that if $x$ is the true k-space data, then operating on it with the SPIRiT kernel should do nothing \[ G x = x. \] If we collect terms on the left \[ (G-I)x = 0 \] Hence the solution $x$ is in the null space of $(G-I)$. The operator $(G-I)$ is determined by calibration. The vector $x$ consists of known and unknown elements. The reconstruction problem then reduced to finding the unknown elements of $x$ such that $x$ is in the null space of $(G-I)$. Two approaches to this are \begin{quote} ``An Eigen-Vector Approach to AutoCalibrating Parallel MRI, Where SENSE Meets GRAPPA,'' Lustig {\em et al}, Proc. Soc. Mag. Reson. Med. 19 (2011), p 479. \end{quote} and \begin{quote} ``Parallel reconstruction using null operations.'', Zhang J, Liu C, Moseley ME. {\em Magn Reson Med. }2011 Nov;66(5):1241-53. Epub 2011 May 20. PMID: 21604290 \end{quote} One of the interesting features of these approaches is that you can solve for the implicit coil sensitivities used by the k-space algorithms. Implement either of these, and reconstruct the data from the homework. \subsection*{Compressed Sensing} \paragraph*{10. Parallel Imaging and Compressed Sensing} The SPIRiT paper describes how to introduce an additional ${ l}_1$ constraint to the basic SPIRiT algorithm. This allows a compressed sensing to be used to allow accelerations beyond what parallel imaging alone would support. This is known as $l_1$-SPIRiT. Implement $l_1$-SPIRiT to show that you can achieve higher acceleration factors than parallel imaging alone would allow. \paragraph*{11. Variable Density Sampling} In class we asserted that for compressed sensing and parallel imaging, the sampling density should be matched to the power spectrum of the image being sampled. In general, the power spectrum an image is of the form \[ P(k_r) = \frac{1}{k_r^\alpha} \] where $k_r$ is the radial dimension in k-space, and $\alpha$ is typically 2, and can range from 1 to 3. For the data set of Assignment 7, study the convergence rate and error when we choose different values of $\alpha$. How important is it to match the sampling density to the image power spectrum? How accurate should the estimate of the power spectrum be? \end{document}