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


\def\hdr {{\parindent 0pt  \vskip -0.7in%
EE369C  Fall 2020-21 \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

\begin{center}
{\Large\em Assignment 7} \\
{\em Due Nov. 19} \\[0.15in]
\end{center}

\paragraph*{Assignments, and the Class}

This will be the last assignment for the class.  If you do a project, you can skip assignment 6 and 7.  Note that the due date is in two weeks, which is the last day of class.

\paragraph{Projection Reconstruction for CT}

For this assignment we will reconstruct fan beam data acquired from a commercial CT system.  The data is from Sarah Patch when she worked at the GE Global Research Center (she is now a Prof. of Physics at Marquette Univ.).  It was acquired on a multi-detector system.

We will do this with a rebinning reconstruction.  In order to do this we will need a parallel beam back projection algorithm, which we will do first.  

\paragraph{1) Parallel Beam Backprojection}

We will start with parallel beam back projection, and simulated parallel beam data of the famous Shepp-Logan phantom.  The phantom and the sinogram data look like this
\begin{center}
 \includegraphics[width=2 in]{slp.pdf}  \hfil
\includegraphics[width=1.25 in]{pd.pdf} 
\end{center}
The data is in the matlab file \verb+sl_phantom.mat+, and has the variables \verb+slp+ which is an image of the phantom, and \verb+pd+ with is the projection data.  The projections are the columns of the matrix (the transpose is shown above).  There are 256 samples per projection and 402 projections over 180 degrees.

To reconstruct the projection data you need to first rho filter each column, and then backproject the result.  Write a matlab function for rho filtering.  Then use the \verb+back_projection.m+ function from the web site to do the back projection.  The image should look pretty close to the image above.  

Include the reconstructed image, and also a reconstruction without doing the rho filter.  


\paragraph{2) Display the Fan Beam Sinogram}

Next we will look at the fan beam data.  There are 888 samples (detectors) over an arc of about 55$^\circ$.  A total of 984 projections are acquired over a full $2\pi$ rotation. The data is in the matlab file \verb+fan_beam_data.mat+ available on the web site via the link.  The data is the variable \verb+d+.  The fan angle in radians is in the variable \verb+fan_angle+. 

The data is quite pretty.  Make an image of the sinogram, and include it. You may have to adjust the window to make it show up well.


\paragraph{3) Parallel Beam Reconstruction of the Fan Beam Data}

Before we get to the fan beam algorithm, it is interesting to try the parallel beam reconstruction directly.  

First rho filter the sinogram, and display the data, again windowing down so that you can see most of the data (there are a few bright areas that will be clipped).  Note that the SNR is noticeably worse in some areas than others.  Compare the low SNR areas to the sinogram you plotted above.  Why is the SNR lower in these areas?

Next, backproject the rho filtered data. Since we have $2\pi$ of projection angles, choose any interval of $\pi$, and perform the backprojection.  The data will support a reconstruction of 512 samples. Display the result.  You should be able to tell there is something in the data, but it won't be pretty.

\paragraph{4) Rebinning}

Next is the rebinning operation.  What we will do is compute a new parallel projection data set for each $\theta$ that we have.  In this case we have about six $\gamma$ samples for each $\theta$.  First, define
\begin{verbatim}
>> [nsamples nprojections] = size(d);
>> dtheta = 2*pi/nprojections;
>> dgamma = fan_angle/nsamples;
\end{verbatim}
For each projection angle $\theta$ we will go through the fan beam data to produce a parallel projection.  The projection angle $\theta_p$ for  the $n^{th}$ fan beam set and the $m^{th}$ ray is
\[
\theta_p = n \Delta \theta + m \Delta \gamma
\]
where we assume $m$ goes from \verb=-nsamples/2+1= to \verb+nsamples/2+.   Recall from the notes that the rebinning operation essentially resamples the fan beam sinogram along diagonal lines.  This is the equation of that diagonal line.  The first term is the constant offset for a parallel projection we are creating, and the second term is the linear offset with fan angle.

To calculate the linear offset, consider the projection at angle $\theta_p = 0$.  Then
\begin{eqnarray*}
 n \Delta \theta  & =  & - m \Delta \gamma \\
 n &  = & -\frac{\Delta \gamma}{\Delta \theta} m
 \end{eqnarray*}
 We can then use this to compute any parallel projection by adding the projection angle constant offset.  If we define \verb+ndth+ for $m$,  the parallel projection data for a sample \verb+ns+ and parallel projection \verb+np+ is given by
 \begin{verbatim}
 >> dp(ns,np) = d(ns,np+ndth(ns));
 \end{verbatim}
 where \verb+ndth+ should be rounded to an integer to keep matlab from complaining. The sinogram extends over $2\pi$, so you need to take care of the end cases, where the index will wrap around (use the modulo function \verb+mod()+). This is essentially a nearest neighbor interpolation.
 
Show the rebinned sinogram.  It should look similar to the sinogram above.  The main difference is that the trace of particular bright objects now look like sinusoids, before they were distorted depending on where they were in the field of view.

Make sure you check for discontinuities or gaps.  This will help you debug your code.


\paragraph{5) Backprojection}

There is one more correction we should do, the correction for the non-uniform spacing of the parallel projections, but we'll put that off for a moment.

Use your parallel beam backprojection algorithm to reconstruct the rebinned data.  Again, use any set of $\pi$ projections.  Show the reconstruction.  This should look pretty nice.

\paragraph{6) Sample Density Correction}

Finally, we should correct the rebinned parallel beam data for the non-uniform beam spacing.  This is a fairly small effect for this data, but it still helps.

The fan-beam rays are at
\[
m_f = \sin(m \Delta \gamma)/(\Delta \gamma).
\]
measured in beam widths, and \verb+m=[-ns/2:ns/2-1]+. If \verb+mx = max(abs(mf))+, what we want is to resample the data uniformly over the range \verb+[-mx:(mx-1)]+, where the \verb+-1+ is to make the number of samples even.  Use the matlab routine \verb+interp1()+ to do the resampling. For a projection \verb+dp(:,jj)+ the corrected data \verb+dc(:,jj)+ is
\begin{verbatim}
>>  dc(:,jj) = interp1(mf,dp(:,jj),[-mx:(mx-1)]);
\end{verbatim}
By default, \verb+interp1()+ uses linear interpolation, which is probably good enough here.  Note that the number of samples in the corrected projections will be smaller that what you started with.

Reconstruct the result, and compare to your previous reconstruction.  There should be fewer artifacts, and the resolution and definition of structures towards the edge of the FOV should be better.


\end{document}