\documentclass[11pt]{article}
\usepackage{palatino,graphicx}
\pagestyle{headings}

\def\hdr {{\parindent 0pt  \vskip -0.7in%
EE369C  Fall 2020-21 \hfill \\
Medical Image Reconstruction \hfill  \\
\vskip 0.1in }}

\frenchspacing

\topmargin=-0.5in
\textheight=9in
\evensidemargin 0.0in
\oddsidemargin 0.0in
\setlength{\textwidth}{6.5in}

\begin{document}
\
\hdr


\begin{center}
{\Large\em EE369C:  Assignment 4} \\[0.15in]
{\em Due Thursday, Oct. 22}
\end{center}


The problems this week will be concerned with 2DFT SENSE reconstruction. The data is an axial brain scan with an eight channel head coil.  It is available here:
\begin{verbatim}
http://ee369c.stanford.edu/data/brain_8ch.mat
\end{verbatim}
The sensitivities were estimated from the data. The data and coil sensitivities are in the following variables
\begin{itemize}
\item \verb+im(160,220,8)+  the complex, fully encoded images for each  channel
\item \verb+map(160,220,8)+  the estimated complex sensitivities for each channel
\end{itemize}
This data and the maps were generated by Miki Lustig.

\paragraph*{1. Displaying the Data}  We have lots of channels, so you will want to display images as an array.  One way to do this is use the \verb+montage+ function which is in the image toolbox.  Unfortunately, \verb+montage+ wants a 4D array.  You can do this with \verb+reshape+ like this:
\begin{verbatim}
>> montage(reshape(abs(im),160,220,1,8), [0 imax]);
\end{verbatim}
where we are displaying the absolute value of \verb+im+ over a range of 0 to \verb+imax+.  Note that the phase encoding direction of the data is the up/down direction, since it is usually taken to be the shortest dimension.  If we were displaying this image for clinical purposes, we would rotate it by 90 degrees, to get it into the standard orientation.

Show the coil data, and the image data using and appropriate range to see the brain (the subcutaneous fat around the edge of the head should be saturated).  Note the variation in coil sensitivities.

\paragraph*{2. Multicoil Reconstruction}

The next task is to compute a reconstruction of the data using the full multicoil reconstruction we described in class, and the square root of the sum of squares reconstruction.  Each of these should only require about one line of matlab code!  A useful command here is \verb+sum(x,nd)+ which sums an array \verb+x+ along the \verb+nd+ dimension.  

Display both reconstructions.  Ideally the multicoil reconstruction should have slightly better SNR, and be more uniform, although it is hard to tell here.  

 \paragraph*{3. $g$-Factor Maps}
The most important characteristic of a SENSE acquisition is the geometry factor $g$.  This tells you how well conditioned the reconstruction problem will be. 

Write an matlab function that determines the $g$-factor when undersampling in the x and/or y dimensions:
\begin{verbatim}
g = gfactor(map,Rx, Ry)
\end{verbatim}
where \verb+map+ are the coil sensitivities, and \verb+Rx+ and \verb+Ry+ are the acceleration factors in the x and y dimensions.  For convenience, assume that \verb+Rx+ and \verb+Ry+ divide evenly into the image dimensions.

Since the coil sensitivities were estimated from the data, the sensitivity maps only exist over the head itself.  This means that a given voxel can be aliased with a variable number of other voxels, and this can be confusing to keep straight. 

One option is doing the calculation for each source voxel independently.  This takes longer, but is easier to keep straight.  Make the sensitivity of the voxel you are interested in the first column in the \verb+C+ matrix.  The aliased voxels are separated in the image by multiples of \verb+Nx/Rx+ and \verb+Ny/Ry+ if the image is \verb+Nx+ by \verb+Ny+ in size.  Append the sensitivities of any aliased voxels that have non-zero sensitivities to the \verb+C+ matrix.  Then when you compute the g factor, just look at the (1,1) element that corresponds to the source voxel you are interested in. That way you don't have to explicitly keep track of how many other voxels are aliased, and where they are.  

Compute and display the g-factor maps for \verb+Rx+ for 2, 3, and 4, for \verb+Ry+ for 2, 3, and 4, and for \verb+Rx=Ry=2+.  For the factor of 3, reduce the size of the coil data by one to make it evenly divide by 3. Display the results on a scale of \verb+[0 5]+.  Display 0 if there is no sensitivity map at a voxel.

\paragraph{4.  Generate Undersampled Images}

Write a matlab function that takes the fully sampled data, and generates the undersampled aliased images.
\begin{verbatim}
ima = undersample(im,Rx,Ry);
\end{verbatim}
This will take the image data, do a 2D DFT, zero out lines in the frequency domain, and then inverse transform to produce the aliased images.  For a factor of 4, you will zero out three out of every four lines.  For \verb+Rx=Ry=2+ you will zero out every other row and column.

Test your routine on the $Rx=Ry=2$ case, and display your result.

\paragraph{5. SENSE Reconstruction}  

Finally, write a matlab function that takes the coil sensitivity maps and the aliased, undersampled images from the previous part, and computes the SENSE reconstruction of the image:
\begin{verbatim}
im = sense(ima, map, Rx, Ry)
\end{verbatim}
The program logic is identical to the gfactor function you wrote above. You just have to add the calculation of the image value for each voxel.

Reconstruct and display images for \verb+Rx+ of 2 and 4, \verb+Ry+ of 2 and 4, and \verb+Rx=Ry=2+.  As a test, try reconstructing the \verb+Rx=Ry=1+ case.  This should give you exactly the same thing as the multicoil reconstruction for part 2. The full multicoil reconstruction is often called the $R=1$ SENSE reconstruction. At low SNR's this does much better than the square root of sum of squares reconstruction.





\end{document}