\magnification=\magstep1
\centerline{Mathematics 372 -- Numerical Linear Algebra}
\centerline{Lab/Problem Set 11 -- Landweber Iteration and Computed Tomography}
\centerline{April 23, 25, 27, 2007}
\bigskip
\noindent
{\it Background}
\bigskip
We are now ready to work with one iterative method that applies for overdetermined
systems and that can be used for the image reconstruction problems from mathematical
SPECT tomography.   Recall that the image acquisition process can be modeled by 
the following equations.  Let ${\bf f}$ represent the true concentrations of the
radioactive tracer in a 2-dimensional axial slice
of the subject (e.g. a patient's torso).  We take a discretization (``pixelization'')
of this slice, so $f$ will be a real vector with $J$ components for some $J$
($f$ can be thought of as either a single, ``long,'' vector, or as a 2-dimensional
array of pixel values which can be plotted as an image). 

The process of data acquisition with multiple angular positions of the radiation
detector camera is modeled by the matrix-vector equation 
$${\bf g} = H{\bf f}\leqno(1)$$
where ${\bf g}$ is a (known) vector with $I$ components that represents the 
measured data by the detector, and the {\it projection matrix}
$H$ is $I \times J$.  ${\bf f}$ is the unknown here -- we want to solve
for it to reconstruct the image from the projection data ${\bf g}$.
Typically in the applications, $I > J$, so the system
above is {\it overdetermined}.  

The problem is to find the ``best fit'' $f$ for the data $g$, in the {\it least 
squares sense}, since (1) typically has no exact solution ${\bf f}$.  
Of course systems like this could be solved in theory either by solving the 
normal equations
$$(H^t H) {\bf f} = H^t {\bf g}\leqno(2)$$
directly, or by using a $QR$-factorization, or the SVD of $H$, as we have discussed.
However, in realistic applications, $I$ and $J$ are quite large (in the 1000's),
so these approaches can be infeasible.   Indeed, the MATLAB programs we will
use never actually write down the whole {\it projection matrix} $H$ or the 
whole {\it back-projection} matrix $H^t$ explicitly.  

Instead, we will use an iterative method called {\it Landweber iteration}.  
Given ${\bf f}^{(0)}$, we iterate 
$${\bf f}^{(k+1)} = {\bf f}^{(k)} + \alpha H^t \left({\bf g} - H {\bf f}^{(k)}\right), \quad 
k \ge 0\leqno(3)$$
for some appropriate choice of the parameter $0 < \alpha < 1$, called the {\it relaxation
parameter}.  Note that (3) can also be written in the standard fixed-point form
as follows:
$${\bf f}^{(k+1)} = (I - \alpha H^t H){\bf f}^{(k)} + \alpha H^t{\bf g}.\leqno(4)$$
The equation
(3) shows that Landweber can
can be thought of as an {\it iterative technique for solving the normal equations}
(2), since if $H^t {\bf g} = H^tH{\bf f}^{(k)}$, then ${\bf f}^{(k+1)} = {\bf f}^{(k)}$, 
and we are at a fixed point.
\bigskip
\noindent
{\it April 23 -- Understanding Projection and Reconstruction Visually}
\bigskip
We will be using a series of MATLAB programs written by Prof. Soares of our
department and a research student of his last year.  To obtain copies of them, enter
the following Linux commands in a Terminal window:
\bigskip
\leftline{\tt mkdir Tomo}
\leftline{\tt cp /home/fac/little/public\_html/NLA07/Tomo/* Tomo}
\leftline{\tt cd Tomo}
\leftline{\tt matlab \&}
\bigskip
\noindent
The first line creates a new folder in your directory, the second copies the program files
we will be using,  the third changes your current directory to the new folder, and the 
last launches MATLAB.

\bigskip
\noindent
{\it Note:}  To be safe, always get into MATLAB this way when working on this lab 
assignment.
\bigskip
\noindent
In the MATLAB command window, enter 
\bigskip
\leftline{\tt Landweber\_demo}
\bigskip
\noindent
This should launch a window containing three plots called Original Image, Projection Data, 
and Reconstructed Data.  

\bigskip
\item{$\bullet$} The default image shown on the left is a ``pixelized'' circular
disk.  
\item{$\bullet$} There are several choices for the color map used in displaying images -- I like
``Bone'' the best personally, but you can experiment to find one that ``works for you.''
\item{$\bullet$} In the left portion, you also get to specify a number of projection angles, and 
the way the images are displayed.  Experiment with these to see what they mean.
\item{$\bullet$} The Project button on the left does the projection step 
$${\bf f} \mapsto {\bf g} = H {\bf f}.$$
The resulting image can be hard to interpret, so don't spend a lot of time 
trying to understand it -- just realize that the process used to create it is
a version of the data acquisition process we discussed -- multiple angular
positions of the detector camera.  
\item{$\bullet$}  The center inputs let you specify the relaxation parameter $\alpha$ for 
Landweber iteration directly.  Or you can have the program estimate the maximum
$\alpha$ for which the iteration will converge for you.  If you press the Find Max Alpha
button, you will need to {\it wait} for the output (this computation can take a while,
depending on the size of the image, so be patient).  
\item{$\bullet$}  You set the Maximum Iterations 
number by typing in a value.  Try the value Maximum Iterations = 3 to start.
\item{$\bullet$} We will not be using the Counts feature at all -- just leave
it with the default value of 100000.
\item{$\bullet$}  The Add Noise button adds random ``noise'' to the projection 
data -- this is a fundamental aspect of using these techniques in the 
real world.  Gamma rays can be deflected by other structures in the body, etc.
and that can introduce this kind of noise in the data collected by the detector.  
You can press this button repeatedly to add ``lots of noise.''  
\item{$\bullet$}  When you press the Reconstruct button, the Landweber iteration
will be started and a message will be printed out in the MATLAB Command Window
when each iteration is begun.  The reconstructed image is displayed on the right.
\bigskip
\noindent
The default image is contained in a file called 
\bigskip
\leftline{\tt disk\_p32p32.mat}
\bigskip
\noindent
which is one of the files you copied.  There are several other files you 
should look at as well:
\bigskip
\leftline{\tt gaussian\_p64p64.mat}
\leftline{\tt gaussian\_offcenter\_p64p64.mat}
\leftline{\tt mcat\_p64p64.mat}
\bigskip
The last of these is a ``mathematical cardiac torso phantom'' -- an artificial
image representing a slice through a human torso showing the chambers of the
heart, spine, ribs, etc.   You can use the Load Image button on the left to 
get these.  (If you launch MATLAB from the directory containing these,
they should show up in the Pick a File dialog box.)  Experiment with all
of these and try to get a good intuitive feeling for how the reconstruction
quality is affected by things like the value of $\alpha$, the maximum number of
iterations used, whether noise is added, what ``chopping negatives''
means, and so forth.  
\bigskip
\noindent
{\it Theoretical and ``Conceptual'' Questions}
\bigskip 
\item{A)}  If you fix $\alpha$, does increasing the number of iterations 
{\it always} improve the reconstructed image quality?  If so, why?  If not, 
what else does this depend on (give a specific example where it does not 
improve things and explain what is ``going wrong'').
\item{B)}  If you fix the maximum number of iterations, 
how does the value of $\alpha$ used seem to affect
the reconstructed image quality?  
\item{C)}  Should the maximum value of $\alpha$ computed by the program
depend on {\it which image} you have loaded?  If not, what {\it does it depend on}?
(Note:  the maximum $\alpha$ is computed by a program called 
{\tt powermethod\_LW.m} that is part of the ``suite'' of programs
you copied.  You can examine the source code with an editor 
if you like.)  
\item{D)}  In fact, prove using our analysis of convergence of iterative
methods that Landweber iteration will converge for all $\alpha$
satisfying 
$$0 < \alpha < {2\over \sigma_1^2},$$
where $\sigma_1$ is the largest singular value of the projection matrix $H$.
(Since we usually do not {\it have} the actual projection matrix $H$, this is primarily
a theoretical result, not a way to determine what value of $\alpha$
to use.)  Can you see what the {\tt powermethod\_LW.m} function is
actually doing now? 
\item{E)}  You might have noticed that if we had the exact 
original object ${\bf f}$ and the exact projection data ${\bf g} = H{\bf f}$,
then the system $H{\bf x} = {\bf g}$, while overdetermined, would have
at least one {\it exact} solution ${\bf x} = {\bf f}$.  
So, strictly speaking, least squares is unnecessary.  However in the 
real world, there {\it will} be noise in the projected data.  We model 
the noise as an additional term added to ${\bf g}$, so we
are actually looking at $H{\bf x} = {\bf g} + \delta {\bf g}$.  Why will 
least squares usually actually be necessary in this sort of problem?
\item{F)}  What is ``chopping negatives'' and why is it reasonable
to do that in looking at the reconstructed image?  (It may help 
\item{G)}  Since the projection data ${\bf g}$ and the final reconstructed 
image ${\bf f}^{(maxiter)}$ can be seen as ``long'' vectors, one way to 
measure the {\it quality} of the reconstructed image would be
to use a vector norm and compute
$$\| {\bf g} - H{\bf f}^{(maxiter)} \|.$$
How would this idea depend on which vector norm we used? 
Which norm do you expect would be closest to saying ``the reconstructed
image looks like the original data'' intuitively?  For 
instance what would it mean if $\| {\bf g} - H{\bf f}^{(maxiter)} \|_\infty$
was small?  Similarly, what would it mean if 
$\| {\bf g} - H{\bf f}^{(maxiter)} \|_1$, or 
$\| {\bf g} - H{\bf f}^{(maxiter)} \|_2$ was small.
Which norm is Landweber iteration set up for?
\bigskip
\noindent
{\it April 25, 27 -- Further MATLAB investigations}
\bigskip
\item{A)}  A ``normalized'' version of the measure of reconstructed
image quality from conceptual question G above is called the 
{\it root mean-square error} (RMS).
Let 
$$\hat{{\bf g}} = H{\bf f}^{(maxiter)}$$
be the {\it estimated} projection data based on the final
reconstructed image.  Then we define:
$$RMS_{data} = \sqrt{{1\over I} \sum_{i = 1}^I (g_i - \hat{g}_i)^2} = 
{1\over \sqrt{I}} \| {\bf g} - \hat{{\bf g}}\|_2.$$
To study this in more detail, we can use another one of the programs
you copied at the start of the lab:  {\tt LWBRMSiter.m}.
\bigskip
\noindent
\item{} To run this you enter a command like this from the
MATLAB command window:
\bigskip
\centerline{\tt LWBRMSiter('mcat\_p64p64.mat',20)}
\bigskip
\noindent
\item{} The first argument is the filename for the image you want, the 
second is the maximum number of iterations.  The output will be
a plot showing how $RMS_{data}$ changes as $maxiter$ increases
from 1 to 20.  That is, the 20 in the input command is the ``maximum $maxiter$.''
The value of $\alpha$ returned by {\tt powermethod\_LW.m} is used here.
Generate plots for each of the 4 image files above and max $maxiter$ 20.
What is happening in each case?
\item{B)}  (This is a somewhat open-ended question.)  If you fix $maxiter$ 
to begin with, what is the best strategy for choosing $\alpha$ to give 
the smallest $RMS_{data}$ (intuitively, this would be the $\alpha$
with the fastest convergence rate)?   You may assume if you like that
$H^tH$ is invertible.  You can attack this either theoretically or
experimentally.  For the experimental way, for instance, you might 
create a modified version of the function from {\tt LWBRMSiter.m}
(call it something else) that lets you supply a value for $\alpha$ in 
addition to a $maxiter$ and use it to gather data about the question.
Note: you will have to think about how to set up the output
to help you answer this.  The output in the function supplied is 
geared to displaying a plot for $RMS_{data}$ after 1 iteration, after 2 iterations,
and so on up to a max $maxiter$ number of iterations.  

\bigskip
\noindent
{\it Note}
\bigskip
This is the last regular lab/problem set for the semester.  It will be due
in class on Monday, April 30.  In the course syllabus I said that the 
final lab project would count for 15\% in the 40\% of the course grade
based on the labs and problem sets.  But since we actually had two
more regular problem sets than I originally planned (including the last two on
eigenvalue and iterative methods), I am not going to do it that
way.  This will just be an ``ordinary'' problem set and it will carry
the same weight as the others.  
\bye

 
