\documentclass[oneside,titlepage]{report}
\usepackage{a4}
\usepackage{aeguill}
\usepackage{makeidx}
\usepackage[T1]{fontenc}
\usepackage[english]{babel}
\usepackage{graphicx, fancyheadings, tabularx}
\usepackage{moreverb, latexsym, amssymb, amsfonts, multicol}


%\newenvironment{code}
%{\begin{list}{}{\setlength{\leftmargin}{1em}}\item\scriptsize\bfseries}
%{\end{list}}

%\newenvironment{code}
%{\begin{list}{}{\setlength{\leftmargin}{1em}}\item\tiny\bfseries}
%{\end{list}}

\bibliographystyle{acm}

\makeindex

\author{Jonny~Reichwald}
\title{Evaluation of Automatized Methods for Coregistration and Interpolation of PET and MR Brain Images in Theory and Practice\thanks{Supervisors: Per~Karlsson and Fredrik~Georgsson}}
\date{30th April 2004}

\begin{document}
\maketitle
\pagestyle{empty}
\newpage
\mbox{}
\pagenumbering{arabic}
\begin{abstract}

There exist several ways to study the brain \textit{in vivo} (within a living organism). These include Magnetic Resonance Imaging (MRI), functional Magnetic Resonance Imaging (fMRI) and Positron Emission Tomography (PET). With the PET method, it is possible to see where different neuroreceptors are located.

Unfortunately, the resolution of the PET camera is very low, which makes it difficult to study small regions of interest in the brain. There are several methods to enhance the image quality, one of them is to coregister the image against an MR image, which can be thought of as laying the images over each other and trying to make them overlap as well as possible. When this is done, it is possible to draw regions of interest on the MR image, and then use them to see how much activity it is in the same region in the PET image.

At the Karolinska Hospital the process of drawing regions of interest is done manually, a very time consuming process which is suspectible to possible bias from the person drawing the regions. The work presented in this thesis will concentrate on evaluating different methods that can be used for the coregistration process, how different interpolation schemes affect the resulting data and compare available programs for this. Especially different interpolation schemes are looked upon, with the surprising result that the most commonly used, trilinear, is the worst performer.

\end{abstract}
\newpage
\mbox{}
\newpage
\pagestyle{headings}
\setcounter{page}{1}
\renewcommand{\thepage}{\roman{page}}
\tableofcontents
\newpage
\renewcommand{\thepage}{\arabic{page}}
\pagestyle{headings}
\setcounter{page}{1}
\chapter{Introduction}
\section{History}


The Positron Emission Tomography (PET) technique was first developed in the early 70's. In 1974 a joint project between the Karolinska Institute, the Royal Institute of Technology and Stockholm University began, and in the beginning of the 80's it was the most outstanding PET centre in the world.
In the beginning, PET was of particular use for studying Cerebral Blood Flow (CBF), and the metabolism in the brain. These areas are nowadays mostly studied with functional Magnetic Resonance Imaging (fMRI). Instead, PET has become an important tool for studying the neurotransmitter system of the brain.

\section{Applications}


PET imaging is used to study where in the brain different neuroreceptors are located and in which quantity they exist. 
With this data, it is possible to see how it correlates with different hypothesises. 
It can be used to see how effective selective serotonin reuptake inhibitors (SSRI) are \cite{pindolol2001}, how different personality traits correlates with different amounts of neuroreceptors, for example how novelty seeking correlates with D2 receptors \cite{suhara2001}, how schizophrenia physically affects the brain compared to a healthy individual \cite{karlsson2002} or how aging affects neuroreceptor concentrations \cite{meltzer}. It can also aid in the search for new drugs,  and psychological correlations with the aid of statistical analysis of changes in the brain.

\section{Goal of this work}

To enhance the data produced by PET imaging, it is preferrable to coregister PET images with MR images, which makes it possible to extract activity information from the PET data in regions anatomically defined on the MR data. Previously, the methods for doing this automatically has not been thoroughly tested at the Karolinska Hospital. The work presented in this thesis describes the different methods in chapter \ref{theory}, compares them with real world tests in chapter \ref{comparison}, and describes the results in chapter \ref{results}. A closer look at what happens with the image voxel values when processing the images is done in parallel.

\chapter{Theory}
\label{theory}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/synapse.eps}
\caption{A: axon, B: neurotransmitter 1, C: neurotransmitter 2, D: synaptic cleft, E: synapse for neurotransmitter 2, F: dendrite.}
\label{fig:synapse}
\end{figure}

The brain consist of around a hundred billion nerve cells, or neurons. These are in turn made of four regions, the cell body, the dendrites, the axons and the terminal. 
The axons from one neuron is connected to the dendrites from other neurons, but there is a small gap between them. To communicate over this gap the axons emit neurotransmitters from the terminal which is the end of the axon.
These are received by the receptors, synapses, on other neurons dendrites (see fig. \ref{fig:synapse}). There exists around 300 known neurotransmitters, of which the more common ones are dopamine, adrenaline, noradrenaline, serotonine and acetylcholine \cite{hubinfaq}. 
When conducting a PET study, a synthetic neurotransmitter, or ligand, is manufactured and then labeled with an radioactive agent. The radioligand will bind to free receptors in the brain. 


\begin{table}
\begin{center}
\begin{tabular}{llr@{.}l}
\multicolumn{3}{c}{Radioisotopes} \\
\hline
Isotope & Name & Half-life (in seconds) \cite{toi} \\
\hline
\hline
$[^{11}C]$ & Carbon-11 & 1223&0\\
$[^{15}O]$ & Oxygen-15 & 122&24\\
$[^{13}N]$ & Nitrogen-13 & 597&9\\
$[^{18}F]$ & Flourine-18 & 6586&2
\end{tabular}
\end{center}
\caption{Some of the more common labeling agents used in PET studies.}
\label{isotopes}
\end{table}

Because of the radioligands short half-life, it must be manufactured close to the PET facility to keep the radioactivity at sufficiently high levels. 
Therefore the radioligand is manufactured locally, in a cyclotron, which is a kind of particle accelerator with which neutrons are collided with the agent producing the radioactive isotopes.

When the radioligand decays it emits a positron ($e^{+}$), which will lose momentum while traveling through the brain tissue, and eventually it collides with an electron. This results in the annihilation of both and two gamma ray photons with an energy of $511KeV$ are emitted, in almost opposite directions (fig. \ref{decay}).
\begin{figure}
\centering
\includegraphics[width=7cm]{figures/decay.eps}
\caption{2D illustration of the PET camera. A: a radioactive decay occurs which emits a positron. B: the positron collides with an electron and emits two gamma ray photons. C: The photons are detected by photomultipliers.}
\label{decay}
\end{figure}

The gamma rays are detected by the tomograph with the aid of small devices able to convert photons to electrons. These devices are called photomultipliers, and they also amplify the current making it possible to detect single photons.
When a photomultiplier detects a decay, a time window is open for a limited amount of time, for a corresponding photomultiplier to detect a photon in some other area. When both has made a detection, the time and position are stored.

The photon signals are not perfect, there is always noise, which is reduced by doing Monte Carlo simulations of the expected decay levels, and data considered to be noise is discarded. The rest of the data is saved in sinograms for later reconstruction. The data files are called sinograms because of how they look, a point that is not in the center of the tomograph will look like a sine wave in the acquisition file if it is stored as an orderered matrix of all the volumes of response (VOR) and viewed with the same ordering.

\section{Reconstruction}

\begin{figure}
\centering
\includegraphics[width=7cm]{figures/reconstruction.eps}
\caption{Illustration of the backprojection reconstruction algorithm in 2D. A: The decay emits photon pairs which are normal distributed around the angle $180$ degrees, sometimes they are emitted in other relative angles, which is depicted in the Line of Response (LOR) C. B: The ideal case of a decay resulting in the LOR D.}
\label{fig:reconstruction}
\end{figure}

%VISA ETT SINOGRAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAM


With the sinograms as input the final images are created with a reconstruction\index{reconstruction} algorithm, the most commonly used is the filtered back projection. By increasing the values for the voxels intersecting the VORs, an image is created (fig. \ref{fig:reconstruction}.) The final spatial resolution of the PET image is limited by several factors such as detector size, positron range of the isotope of interest, undersampling of the signal in the linear or angular directions for the image reconstruction process, patient motion and Compton scatter \cite{levin1999}.
An extensive description of the filtered back projection algorithm in 3D is given in \cite{bendriem1998}. 

\section{Transform Models}
\label{transmodel}
An archetypical form can be transformed into shapes that differs greatly, as D'Arcy Thompson showed by taking the shape of a fish and changing it into another species by transforming its shape \cite{thompson}. Probably there is no transform for each brain into an archetypical standard brain, but some attempts have been made to make a standard atlas brain to which other brains are transformed \cite{talairach1988}. The process of transforming an image into the shape of another is called coregistration, or registration.

To be able to register images, one of the images must be transformed to the coordinate system of the other image, either by a global or local transformation. A global transformation affects each part of the image in the same way contrary to a local transformation which can affect only a small part of the image. 

\begin{figure}
\centering
\includegraphics[width=8cm]{figures/transformations.eps}
\caption{Examples of global 2D transformations. A is the original image, B is a rigid body transformation, C an affine transformation, D a projective transformation and E is a curved transformation.}
\label{fig:trans}
\end{figure}

\label{trans}
There are several methods to transform\index{transform models} an image, the simplest one is the rigid body transformation\index{rigid body transformation}, which rotates and translates an image, but preserves internal distances and angles. It has six parameters in $\mathbb{R}^3$, three to specify the rotation and three tho specify the translation. With intrasubject, intramodal image registrations, the head can be assumed to be a rigid body since it do not tend to change shape. Therefore a rigid body transform is sufficient for those scenarios.
For intrasubject image registrations with different modalities prior adjustments may be needed, such as scaling and shearing to compensate for warping effects of the MR camera. A rigid body transformation with scaling and shearing is called an affine transformation\index{affine transformation}. An affine transformation preserves colinearity and ratios of distance.

By writing a vector $[x\ y\ z]^T$ as $[Hx\ Hy\ Hz\ H]^T$ it is called a homogenous vector. The factor $H$ is usually set to the identity. By describing points as homogenous vectors it is possible to express compositions of linear transformations as multiplications and inversions of their corresponding matrices.
Both rigid body and affine transformations can be expressed as a composition of linear transformations of the form $T_B \circ T_A$. For the transformations $T_A$ and $T_B$ this can be expressed as $(T_B \circ T_A)(u) = T_B(T_A(u))$. This only holds if $T_A: \mathbb{R}^n \mapsto \mathbb{R}^k$ and $T_B: \mathbb{R}^k \mapsto \mathbb{R}^m$, but since $u$ always stays in the same space, this is not an issue. To reduce the number of matrices with which $u$ has to be transformed with, the matrix representation of the original linear transformations can be multiplied with each other since $(T_B \circ T_A)(u) = T_B(T_A(u)) = B(Au) = (BA)u$. With the knowledge that compositions of transformations is equivalent with the multiplication of their matrices we can express the rigid body transformation matrix as

\begin{equation}
\label{rigid}
T = 
\left(
\begin{array}{cccc}
cos \phi cos \theta & cos\phi sin\theta & -sin \theta & dx\\
sin \psi sin \phi cos \theta - cos \psi sin \theta & cos \psi cos \theta + sin \psi sin \phi sin \theta & sin\psi cos \phi & dy\\
sin \psi sin \theta + cos \psi sin \phi cos \theta & cos \psi sin \phi sin \theta - sin \psi cos \theta& cos \psi \cos \phi & dz\\
0 & 0 & 0 & 1
\end{array} 
\right)
\end{equation}

in $\mathbb{R}^3$. $\theta$ is rotation about the z-axis, $\phi$ rotation about y-axis and $\psi$ is rotation around the x-axis. $dx$, $dy$ and $dz$ specifies the translation in the different planes. With the composition of additional transformation matrices affine transformations can be written in the same form but with 12 parameters instead of 6\footnote{An affine transformation is said to have 12 Degrees of Freedom (DOF) while a rigid body transformation has 6DOF.}\cite{anton1994}.

There exist more advanced transforms which are used when an image shall be transformed to another space, such as projective and curved transformations, which can be seen in fig. \ref{fig:trans} \cite{ashburner2000, woods1993}.
These are useful when the brain must be warped to fit an atlas brain such as the one described by Talairach and Tournoux \cite{talairach1988}, or the MNI brain atlas \cite{evans1993}. In SPM \cite{spm} this is done with a affine transformation followed by an estimation of nonlinear transforms made from linear combinations of 3-dimensional Discrete Fourier Transform (DCT) basis functions \cite{cox1999}. The number of parameters now rise to 1176 with the default settings.

\section{Alignment}

There are regions in the brain that do not have a perfect anatomical definition, the definitions may be vague and rely on other data. One such region is the dorsolateral prefrontal cortex\footnote{See appendix \ref{anadef} if unsure about the anatomical definitions mentioned here.}, which can be defined in the coronal plane as the grey matter with posterior limit at the first slice frontal of the corpus callosum and frontal limit at the most frontal part of cortex. Medially, it goes from the cingulate gyrus laterally down to the last sulcus adjacent to orbitofrontal cortex \cite{raz1997}. The plane defined by being placed frontally of corpus callosum will change orientation depending on the orientation of the image. This will in turn affect the definition of the dorsolateral prefrontal cortex area, and therefore it may be important to keep the brain aligned.

For study types where different brains are transformed to an atlas space, to be compared statistically against other brains it is of great importance that the brains have a good alignment with the target brain before any nonrigid transormations take place, otherwise it becomes warped more than required which will distort the image data more than neccessary.


\begin{figure}
\centering
\includegraphics[width=8cm]{figures/acpc.eps}
\caption{The markers show the posterior commisure (PC) and the anterior commisure (AC).}
\label{acpc}
\end{figure}

To avoid problems like the ones mentioned above, it is preferable to always align the heads to a similar orientation. By using a line which can be imagined to go through the anterior commisure and the posterior commisure -- the AC-PC line, and then using this line as the y-axis (pointing forward), intersecting the x-axis in AC, the number of possible undefined movements of the brain is reduced to two. Further, by using the plane that bisects the brain into two hemispheres and intersects the AC-PC line, and using this as the definition of the zy-plane, the brain is AC-PC aligned and it is not possible to move or rotate the brain in any direction without violating the AC-PC alignment (see fig. \ref{acpc}).

The MR images are usually not exactly aligned to the AC-PC plane, and the first post-processing step is to align it. This is done manually by translating and rotating the image until the AC and PC is clearly seen in the transverse view, and the coronal view is straightened up to make the interhemispherical fissure as vertical as possible. The image is then resliced to the new coordinate system.



\section{Coregistration}

Coregistration\index{coregistration} can be done intrasubjectively or intersubjectively, that is to say, between images from the same person or between different persons. Furthermore, it can be done intramodally or intermodally, the intramodal case is when the images are from the same modality, say two MR images, and the intermodal case is when the modality differs, for example MR-PET.

Coregistration can be done either manually or automatically, or by both methods by first applying an automated method and then adjusting the image manually until the examiner is satisfied. 

The automatic methods often convolve the image with a Gaussian kernel before coregistration to avoid local minima/maxima and to speed up the search.

The manual method is not usable if the image is to be transformed by an affine or even more complicated transform, since it quickly becomes unsustainable for the examiner to keep track of all the parameters.

The images that are being coregistered is referred to as the target image $u$ and the source image $v$, the target image is not transformed and serves only as a model to which the source image is transformed to.

\subsection{Automatic Coregistration Algorithms}

There have been several attempts to develop an automatic coregistration routine, but it is not an easy task when the images are from different modalities. Below are the most successful algorithms described. Other approaches for PET-MR coregistration has been tried, such as surface-based methods \cite{pelizzari1989} and chamfer matching \cite{jiang1992}.

\subsubsection{The AIR and SPM99 Coregistration Algorithm}

For intramodal images that are to be coregistered with a rigid body transform, the parameter estimation can be done by searching for a minimization of square error. The same method can be used for warping transforms, but sensible restraints are required to avoid large local changes. Intermodal images can not be subject to neither of these approaches, at least not directly since the images represent different  tissues. 

In AIR (a program described in section \ref{airdesc}), doing a PET-MR coregistration, the operator is supposed to first strip the skull, scalp and meninges from the MR image, and thereafter applying the parameter minimizing program \cite{woods1993}.

For the early incarnations of the SPM package (see section \ref{spmdesc} for program description) this was solved by automatically segmenting the grey tissue, white tissue, cerebral brain fluid and skull with a clustering algorithm which can exploit the information differences from $T1$ and $T2$ weighted MR images \cite{ashburner1997,hartigan1975}. The segmented brain can then be used in any of the methods mentioned above.

\subsubsection{Coregistration by Maximization of Mutual Information}
Each image voxel $v_i$ represent a symbol in a set $ S = \{v_1, v_2, \ldots, v_m\}$ where $m$ is limited by the number of voxels. Each symbol has an associated probability $P_i$ which can be expressed as 
\begin{equation}
\label{eq:voxelprob}
P_i = \frac{1}{m} \sum_{k=1}^{m}f(v_k)
\end{equation}
\begin{equation}
\label{eq:bin}
f(\cdot) = \left\{
\begin{array}{cc}
1 & ,P(v_k) = P(v_i)\\
0 & ,P(v_k) \neq P(v_i)
\end{array} \right.
\end{equation}

The entropy, the unpredictability of a sequence of symbols drawn from $S$, of these voxels is

\begin{equation}
H = - \sum_{i=1}^m P_i\ log_2\ P_i
\end{equation}

where $H$ is the entropy measured in bits and $P_i$ is the associated probabilites for the symbols $v_i$. A bit is the amount of information it is required to disambiguate an uncertainty which can be resolved by a yes/no answer \cite{duda2001, shannon1948}. 

Eqn. \ref{eq:bin} can be modified so it equals $1$ for a range instead of an exact value, thus becoming 

\begin{equation}
\label{eq:bin2}
f(\cdot) = \left\{
\begin{array}{cl}
1 & , P(v_i) - \epsilon \leq P(v_k) \leq P(v_i) + \epsilon\\
0 & , else
\end{array} \right.
\end{equation}

where $\epsilon$ defines the half width of the window where $P(v_i)$ and $P(v_k)$ is considered equal. This makes real calculations simpler.


The entropy $H$ for voxels from the different images can be used in a dissimilarity function, derived from the Kullback-Leibler distance. The Kullback-Leibler distance does not measure distance in the conventional way, like the Minkowski metric ( eqn. \ref{eqn:minkowski}), but is based on the differences in statistical distributions \cite{duda2001}.

For the two distributions $P(x)$ and $Q(x)$, their distance can be expressed as

\begin{equation}
\label{eqn:kullback}
D_{KL}(P(x), Q(x)) = \sum_x Q(x)\ log_2\ \frac{P(x)}{Q(x)}
\end{equation}

where $D_{KL}$ is the Kullback-Leibler distance. This is not directly useful for registration of images of different modalities, since they have different variables and the Kullback-Leibler distance only measures distances between distributions for the same variable $x$ \cite{wells}.


By introducing a second variable $y$, we can get two distributions, $P(x)$ and $Q(y)$, with different variables.
The concept of mutual information is the amount of uncertainty that can be removed from one variable with knowledge about the other variable.

\begin{equation}
\label{eq:mi}
I(P; Q) = H(P) - H(P|Q) = \sum_{x,y} R(x,y)\ log_2\ \frac{R(x,y)}{P(x)Q(y)} 
\end{equation}

$I$ is the mutual information and $R(x,y)$ is the joint probability of $x$ and $y$. 

\begin{figure}
\centering
\includegraphics[width=10cm]{figures/pet_coreg.eps}
\caption{Joint histogram of two PET images. A: the images are unregistered. B: registered images. The source image is the same as the target image, but slightly rotated. Usually there is one MR and one PET image, but to show how the intensity values align to the diagonal in the histogram when registered it was done with images with the same modality.}
\label{fig:joint}
\end{figure}

The joint probabilites can be calculated by making a joint distribution and then normalizing it. The joint distribution can be calculated by making a joint histogram of the images to be registered. The joint histograms before and after coregistration can be seen in fig. \ref{fig:joint}. The pixel intensity at a point $(x,y)$ in the joint histograms is given by the number of voxels in the area overlapped by both images with a certain value given by the x-axis in one of the images, and the number of voxels with a value given by the y-axis in the other image. The brighter the pixel, the more voxels with those intensities. The bin size for the histogram is choosen somewhat arbitrarily, but earlier work has suggested that an optimal bin size is $W = 3.49\rho N^{-1/3}$, where $W$ is the width of the histogram bin, $\rho$ is the standard deviation of the distribution and N is the number of available samples \cite{scott1979}. The calculated values are shown in table \ref{tbl:scott}. The number of bins in SPM2 and PMOD (further information about PMOD can be found in section \ref{pmoddesc}) is $256$, which probably is a good number since it is preferrable to err on the low side, since the statistics will become more reliable and the execution time will decrease \cite{tsao2003}.


By using the histogram values in eqn. \ref{eq:mi}, a measure of the mutual information is given, which then can be a dissimilarity measure to minimize with any minimization technique with the transformation matrix as parameter.

The method of coregistration by maximization of mutual information was first independently suggested by Viola \cite{viola, wells} and Collignon \cite{collignon}, and has been further improved by Studholme et al. \cite{studholme1998}.


\begin{table}
\begin{center}
\begin{tabular}{cccccc}
\multicolumn{6}{c}{Bin sizes} \\
\hline
Modality & $\rho$ & N & W & max/min  & $n$ bins\\
\hline
\hline
MR & 400.3901 & 10223616 & 6.4383 & 7054/0 & 1095 \\ 
PET & 72.2049 & 770048 & 2.7493 & 1377/-95 & 1412
\end{tabular}
\end{center}
\caption{Optimal bin size calculated for an MR and a PET image.}
\label{tbl:scott}
\end{table}

%PET
%num voxels: 128 128 47 = 770048
%stddev: 72.7049
%varians: 5286.003392 
%mean: 36.3453
%
%MR:
%voxels: 256 256 126 = 10223616
%varians: 160312.203075
%stddev: 400.3901
%mean: 263.9119
%
%
%

\subsection{Manual Coregistration}

There are several methods for manual coregistration, but they all share the attribute that they need an operator to carry out the procedure. These procedures are all very simple, but nonetheless still powerful enough to fulfill the requirements of the coregistration procedure.

Common for all landmark methods is that points are marked on the different images which corresponds to each other, then a transformation is calculated from these points and then applied to the target image.

\subsubsection{Landmarks}
By marking different parts of an image and then marking the corresponding parts in a different image it is possible to calculate the registration parameters by finding a suitable transform for the markers -- $T: u \mapsto v$ which transforms the markers in the first image $u$ to the corresponding markers in the second image $v$. If a linear transform cannot be found, a least squares, or any other optimization method can be applied to find a transformation that best minimizes the distance between the corresponding markers. The optimization procedure tend to be very cheap since the number of landmarks usually are very few, and the calculations are done very fast on any modern computer.
Another method is to let $T$ be a transformation of a higher degree, which is useful for registration against atlases, but this requires many landmarks to make it accurate and is not easily done.

\paragraph{Intrinsic landmarks}
\label{ana_land}
The skull tissue, eyes and other internal landmarks can be used for the coregistration. The problem here is that the PET images are not anatomical, but functional and what is seen on the PET image is not easily translated into something anatomical. Still, it is possible to tell where the brain tissue ends in the PET images with many radioligands, and for some ligands which have a high affinity for certain anatomically well-defined regions, those regions can be used as landmarks.
A downside with intrinsic landmarks is that they must be placed manually which needs user interaction and is very time consuming. 


\paragraph{Extrinsic Landmarks}
Extrinsic landmark registration relies on the ability to place external objects close to the patient. These objects shall be clearly visible in the image, and be attached to the patient so they stay in the same place on all images \cite{bergstroem}. 
One method is to put markers on a helmet worn by the patient (see section \ref{helmet}), which is then seen on the images and can then serve as landmarks which has to be matched. The method used for external landmarks is the same as in the intrinsic landmark method, but since the markers retain their internal distances and angles to each other no advanced parameter estimation techniques are required, although scaling can become neccessary if the images are in different modalities.

\subsubsection{Image Overlay}
\label{overlay}
The source image is viewed over the target image, and can be translated and rotated until the examiner is satisfied. This method will only work with rigid body transformations which makes it unsuitable for intersubject registrations, but for intrasubject registrations it is very good. Any search for the transformation parameters becomes unneccesary, since they already are determined by the examiner.


\section{Interpolation Methods}
\label{interpmethods}
The image data is a discrete representation of continuous data, which can bee seen as data points in a lattice of sample points. Voxel values for an image $u$ is stored as an array $u(x_i,y_j,z_k)$, where $(x_i,y_j,z_k)$ is the index.  $(x_i,y_j,z_k) \in G$, where $G$ is the dimensions of the image $u$. When a transformation $T: u \mapsto v$ is applied, the data points from $u$ will almost certainly be misaligned with the lattice points in $v$. A method to calculate approximate values must then be used.


\begin{eqnarray}
\label{eqn:inv}
Tu & = & v \\
TIu & = & Iv \\
TTT^{-1}u & = & TT^{-1}v \\
TT^{-1}u & = & T^{-1}v \\
Iu & = & T^{-1}v \\
\label{eqn:inv2}
u & = & T^{-1}v
\end{eqnarray}


Since there are no values in the target space to use for interpolation, the inverse $T^{-1}$ can be used which transforms from target space to source space, which is shown in eqn. \ref{eqn:inv} -- \ref{eqn:inv2} \cite{hamre1998}. Now voxel locations $v(x,y,z)$ can be transformed to the source image space by the inverse transformation $T^{-1}v(x,y,z) = u(x,y,z)$, where $(x,y,z) \in \mathbb{R}^3$. These voxels are probably not well aligned either, but now a convolution can be applied to the surrounding data points which gives a value to the given point. This new value is then used as the interpolated value in the new image  at coordinates $v(x,y,z)$. 

In a very general form, interpolation can be described by
\begin{equation}
\label{eqn:interp}
f(x) = \sum_{k \in G} f_k \varphi(x - k)
\end{equation}

$f(x)$ is the interpolated value at any point in the image space and $f_k$ the sample values. $\varphi(x-k)$ decides which weight each sample shall have -- the convolution kernel, or basis function. The difference between the interpolation algorithms mentioned in this text lies in the choice of basis function.


Interpolation is not only useful for registration of images, but is used when the image is upsampled or downsampled too.
Depending on the choice of coregistration method, it is sometimes required that the images are of the same modality.
The PET images at the Karolinska Hospital has a voxel size of $2.202 \times 2.202 \times 3.125\ mm$. The MR images can have different resolutions depending on the type of scan, down to  $1.016 \times 1.016 \times 1.000\ mm$. 

Coregistration methods which operate directly on voxel values need to have the images in the same modality and they can use any of the following interpolation methods for use in their dissimilarity functions. Although they tend to use a simple interpolation method for this task since they need to recalculate the voxel values for each iteration.

There exist other methods than the convolution based methods mentioned here, such as morphology based methods \cite{grevera1996, herman1992}, but they are not included in this text since there is no program known to the author that is currently using any of them.


An in-depth discussion about interpolation is given in \cite{thevenaz2000}.

\subsection{Nearest Neighbour Interpolation}
\label{nn}
\begin{figure}
\centering
\includegraphics[width=6cm]{figures/square.eps}
\caption{Nearest Neighbour basis function.}
\label{fig:nn}
\end{figure}


Nearest neighbour is a method that is very simple, and therefore fast. This makes it appropriate when fast calculations are needed, and accuracy is of a smaller importance. 

By taking the value of the datapoint nearest to the calculated value $u(x,y,z)$, a crude but quick interpolation has been done. To measure which one is the closest, a metric must be used.
\begin{equation}
\label{eqn:minkowski}
D_M(\mathbf{a}, \mathbf{b}) = \left(\sum^3_{i=1}|a_i - b_i|^k\right)^{1/k}
\end{equation}
$D_M$ is the Minkowski distance, $\mathbf{a} \in u$ and $\mathbf{b} \in v$. The most common case is when $k = 2$, which is the Euclidean distance. With $k = 1$, the Manhattan distance is used, which can be computationally even faster, but as usual at the cost of precision.

By substituting $\varphi$ in eqn. \ref{eqn:interp} with a function $NN$ that returns $1$ when $k$ is the closest neighbour, we have constructed a function that can be used for interpolation. A graph of the of the function $NN$ is shown in fig. \ref{fig:nn}.

It is a common belief that nearest neighbour interpolation never should be used for final data, only as a fast method for intermediate images used in optimization. Later on it is shown that this belief may be ungrounded.

\subsection{Trilinear Interpolation}
\label{linear}
Trilinear filtering is the most intuitive and most commonly used method for interpolation. It is a bit more sofisticated than the nearest neighbour method, but still fast to calculate. The basis function for trilinear filtering is shown in fig. \ref{fig:linear}.

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/linear.eps}
\caption{Linear interpolation basis function.}
\label{fig:linear}
\end{figure}

We construct the basis function by taking the Euclidian distance $D$ ($k=2$ in eqn. \ref{eqn:minkowski}) to the currently indexed data point, and subtract it from the distance between samples. The result is used as a weight for the contribution of that data point to the final data.

\begin{equation}
\label{eqn:trilinear}
\varphi = 1 - D_k
\end{equation}

$D_k$ the distance to currently indexed data point from our current position to which we wish to interpolate a value. It is assumed that the distance between data points is 1 in this equation, but other values work as well if they are first normalized.

\subsection{Partial Volume Interpolation}

This method is similar to trilinear interpolation, and is used in conjunction with the mutual information algorithm, it updates the image histogram in a more smooth manner -- for each neighbour, it increases its histogram entry by the weight for that neighbour (eqn. \ref{eqn:trilinear}). This makes the histograms change gradually and no sudden changes can occur which reduces the risk of being trapped in a local maxima when searching for the global maximum of mutual information.

\subsection{Windowed Sinc}

The sinc function is defined by

\begin{equation}
sinc(x) = \left\{
\begin{array}{ll}
\frac{sin(\pi x)}{\pi x} & ,x \neq 0 \\
1 & ,x=0 \\
\end{array} \right.
\end{equation}

\begin{figure}
\centering
\includegraphics[width=7cm]{figures/sinc.eps}
\caption{Sinc function for $-3 \pi \leq x \leq 3 \pi$}
\label{fig:sinc}
\end{figure}

and it is displayed in fig. \ref{fig:sinc}. It is an interesting function when working with signal sampling due to its box-shape in Fourier space, which makes it suitable for image reconstruction from Fourier space signals, for example MR signals.
The sinc function oscillates into infinity and when used as basis function it is said to require infinite support, but the images are spatially limited and can only give limited support. Furthermore, the calculations become burdensome when each voxel value in the whole image has to be included in the calculations. To avoid this, it is possible to introduce a window function that limits the sinc function to an interval:

\begin{equation}
\label{eqn:window}
h(x) =  	\omega(x)sinc(x)
\end{equation}

where $\omega(x)$ is the window function with band limitation given by $m$:

\begin{equation}
\omega(x) = \left\{
\begin{array}{ll}
\omega(x) & ,0 \leq |x| \le m\\
0 & ,m \leq |x| \\
\end{array} \right.
\end{equation}

When trying to estimate a signal that may be discontinuous with the sinc function (or any other continuous signal), ringing will occur. This is because it is not possible to reconstruct a discontinuous signal with only continuous functions, and at the points of discontinuity in the original signal, the estimated signal will ripple. This is called the Gibbs phenomenon. It is possible to reduce it, for example with the Lanczos window described by:

\begin{equation}
\label{eqn:lanczos}
\omega(x) = \frac{(\pi x)}{m}
\end{equation}

There exist a plethora of more or less exotic window functions, for a more in-depth discussion of their properties, see \cite{harris1978, meijering2001}.


\subsection{Splines}


The discrete case of eqn. \ref{eqn:interp} is described by the equation
\begin{equation}
\label{eqn:interpdiscrete}
f_{k_0} = \sum_{k \in G} f_k p_{k_0 - k}
\end{equation}
with $k_0 = x$, and $p_{k_0 - k} = \varphi(x-k)$. 
 
For the methods mentioned earlier, only the basis function has been described, but to make it possible to use other basis functions we must understand that eqn. \ref{eqn:interp} only holds true if the basis function vanishes for all integer arguments except at the origin -- the function is an interpolant. A property that holds true for linear, nearest neighbour and sinc interpolation \cite{unser2002}.

For basis functions with noninterpolant properties,
an alternative form of eqn. \ref{eqn:interp} can be expressed as
\begin{equation}
\label{eqn:interpdiscrete2}
f_{k_0} = \sum_{k \in G} c_k p_{k_0 - k}
\end{equation}

The difference is the introduction of $c_k$ instead of $f_k$, where $c_k$ is coefficients. 
Eqn. \ref{eqn:interpdiscrete2} makes it possible to perform the interpolation in two steps -- first the calculations of the coefficients, and then calculate the values $f(k_0)$ from $c_k$.
If $p$ is known \textit{a priori}, eqn. \ref{eqn:interpdiscrete2} can be viewed as a linear system of equations with $c_k$ as the unknown variables. This can be expressed as $c = P^{-1}f$ and to solve it, $P$ must be inverted. If $p$ is interpolating, then $P = I$, and it is trivial to solve the system. 
For the interpolation methods mentioned earlier, with interpolating basis functions, $P = I$ yields the equation $c = If = f$, and the need for eqn. \ref{eqn:interpdiscrete2} is nonexistant.


Higher degree splines are unfortunately not interpolating, but inversion of $P$ can still be avoided, by recognizing that discrete convolution is equivalent to \ref{eqn:interpdiscrete2}, and can be written
\begin{equation}
\label{eqn:conv}
f_{k_0} = (c * p)k_0
\end{equation}
By using the inverse $(p)^{-1}$ to convolve both sides of \ref{eqn:conv} we get
\begin{equation}
\label{eqn:conv2}
c_{k_0} = ((p)^{-1} * f)_{k_0}
\end{equation}
which can be used to calculate all coefficients.

Splines can be described by the (n+1)-fold convolution of a spline of degree 0 -- $\beta^0$-spline, which is almost identical with the nearest neighbour basis function. From this it follows that the linear basis function is equivalent with $\beta^1 = \beta^0 * \beta^0$. Higher degrees are given by:

\begin{equation}
\label{eqn:splinen}
\beta^n = \beta^0_1 * \beta^0_2 * ... * \beta^0_{n+1}
\end{equation}

Splines of degree $n>1$ are not interpolating and therefore the coefficients in eqn. \ref{eqn:interpdiscrete2} must be calculated. By using eqn. \ref{eqn:conv2}, this can be done by finding the inverse convolution $(b^n)^{-1}$ for the $\beta$-spline basis function $b^n$. Since convolution can be viewed as a kind of filter, the coefficients can be found by filtering the image data with the inverse convolution kernel, and a recursive implementation is given in \cite{thevenaz2000}. 
With the coefficients given by $(b^n)^{-1}$ and basis function $\beta^n$, we can write \ref{eqn:interp} as

\begin{eqnarray}
\label{eqn:cardinalspline}
f(x) & = & \sum_{k \in G}  \left((b^n)^{-1} * f\right) (k)\beta^n(x - k) \\
	  & = & \sum_{k \in G} f(k) \eta^n(x-k)
\end{eqnarray}

where $\eta$ is the cardinal spline
\begin{equation}
\label{cardinaldef}
\eta^n(x) = \sum_{k \in G} (b^n)^{-1} (k) \beta^n(x-k)
\end{equation}

Since \ref{cardinaldef} takes the image data values as parameters instead of coefficients it must have the interpolation property. This is true, and it looks very much like the sinc kernel for higher degrees. As $n \rightarrow \infty $, $\eta^n$ converges to the sinc function \cite{unser1999}. This obliterates the practical need for the sinc function and provides a relatively computationally cheap method for a theoretically good algorithm.


%\cite{unser1993}
%\cite{meijering2001}

\section{Locating activity differences}

The goal with PET imaging is to see if any significant difference of neurotransmitter levels or neurotransmitter bindings is occuring in different parts of the brain. There are two methods for doing this, either with the definition of regions of interest (ROIs) or by statistical parametric mapping.

\subsection{Regions of Interest}

When doing studies with PET, it is hard to see where in the brain certain anatomical areas are, and to do quantitative studies for the bindings in a certain region requires that it is somehow delimited.

A region of interest (ROI) is an area which contains an anatomical structure. This structure can manually be defined by drawing it on an MR image in several planes, which together will define a volume of interest (VOI), or the ROIs can be defined on an atlas brain to which the source image is transformed, and with that it will have ROIs defined. The difference between a ROI and a VOI is often neglected in casual talk, and the two should be seen as interchangeable when mentioned in this text.

Values extracted from ROIs are very stable since they are based on the average radioactivity in a large volume. It is important to remember that the ROIs cannot be defined for regions smaller than twice the sampling frequency, since the sampling frequency is twice as high as the resolution of the camera, to fulfill the Nyquist sampling criterion \cite{nyquist1928}.

It can be of interest of how the voxels in the borders of the VOI is handled. If they are partially outside the volume, should they be included, excluded, interpolated or excluded/included depending on how much of their area is outside? This is an open question which has not been further investigated.

\subsection{Parametric Mapping}

Parametric Mapping requires that the images all are aligned to a standard brain. By looking at the same voxel in all the test subjects a hypothesis for significant change between the reference group and the test group can be tested. This results in an image of the standard brain where significant changes are seen \cite{worsley1996}.

This method allows for the discovery of novel areas where activity changes, and does not require any manual work, at the cost of much noise. There are attempts to reduce the noise, for example by preprocessing the images with a wavelet filter \cite{cselenyi2001},  but this method is still not as stable as the ROI method. SPM is especially developed with parametric mapping in mind.


\section{Time Activity Curves}

After the images are aligned and coregistered, and the ROIs are drawn, it is time to extract the data. This is done by analyzing the voxels which reside inside the VOIs. Their value is averaged for each time frame as a separate value for each volume. The data is then used as parameters in the modelling stage.

\section{Kinetic Modelling}

The time activity curves do not tell how many neuroreceptors there are in the region relative to the surrounding regions. To find out, we must know what the activity is in reference tissue where the concentration of neuroreceptors of current interest are negligible. As reference, cerebellum is commonly used.

The radioligand is injected in tracer amounts, which results in the receptors never being fully occupied by the radioligand. The ligand can be bound to nonspecific receptors and specific receptors. The later is the receptor targeted by the ligand, which is of interest for the study. Because of the tracer amounts, there is a constant flux between the free ligands and the specific/nonspecific bound ligands. It is possible to set up a model that describes the kinetics for the ligand flux, but the parameters are unknown. By assuming that the amount of nonspecific binding is the same in the reference region as in the region of interest, it is possible to estimate the parameters of this model by inserting the time activity values and time \cite{gunn1997}. The end result is a measure called binding potential (BP), which describes the affinity of the ligand for a certain region. 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55



\chapter{Comparison of the methods}
\label{comparison}
In this chapter the programs available for doing coregistrattion, AC-PC alignment and other image manipulations will be compared. Some of them uses the same theoretical methods, but still they produce different output as will be shown.

Tests have been performed to measure speed and accuracy of foremost interpolation methods.

\section{Current methods in use at the Karolinska Hospital}

Table \ref{capabilities} shows the different capabilities of all the programs that where evaluated. No program except PMOD can do all the required tasks and the data has to be moved between different platforms and converted between different formats. Conversion is very time consuming, especially if there is problem with any of the programs.

A typical usage scenario can look like:

\begin{itemize}
\item Convert the PET images from the ECAT6 format to the Analyze format, and the MR images from General Electric format to Analyze format.
\item Extract frame duration times from the original PET images since Analyze is unable to store that.
\item Edit the MR images in ImageJ to remove edge slices which often exhibit aliasing effects.
\item Do AC-PC alignment and coregistration in SPM.
\item Convert PET images from Analyze to the VFF format.
\item Draw ROIs in Hba with great care, lest Hba crashes, and extract Time Activity Curves (TAC).
\item Postprocess statistics and do modelling in Excel with the aid of formulaes describing the BP.
\end{itemize}

Several of the above mentioned steps consist of file conversion. This is hard to avoid if different programs are to be used. All the available programs that can convert images are listed in table \ref{tbl:convert}. Older MR images are in a format used by General Electric in their older MR cameras and there exist a small program to convert these images to the ECAT6 format. New MR images are in the DICOM format. PET data was until recently in ECAT6 format only, and still is occationally. Nonetheless, the tools to convert from the old formats are required when older scans are to be included in a study.

The featureset of the programs used for drawing ROIs differ. An important one is their ability to draw ROIs in different projections. The projections supported by each program are listed in table \ref{draw}.

\begin{table}[b]
\begin{minipage}{\textwidth}
\renewcommand\thefootnote{\it\alph{footnote}}
\begin{center}
\begin{tabular}{llccccc}
\multicolumn{7}{c}{Program capabilites} \\
\hline
Name & Platforms & Alignment & Coreg. & ROI & TAC & Stats. \\
\hline
\hline
Hba & Irix & - & - & x & x & - \\
ECAT6 & Solaris & - & - & x & x & - \\ 
ImageJ & All with Java & x & - & - & - & - \\
mars & All with Matlab & x & x & - & - & - \\
SPM & All with Matlab & - & x & - & - & - \\ 
KaleidaGraph & Windows, MacOS & - & - & - & - & x \\
Excel & Windows, MacOS & - & - & - & - & x \\
PMOD & All with Java & x & x & x & x & x \\
ECAT7 & Solaris & - & - & x & x & -
\end{tabular}
\end{center}
\end{minipage}
\caption{Matrix of programs for different tasks, autumn 2003.}
\label{capabilities}
\end{table}


\begin{table}[b]
\begin{minipage}{\textwidth}
\begin{center}
\begin{tabular}{llcccc}
\multicolumn{6}{c}{File conversion programs} \\
\hline
Name & Analyze & ECAT6 & ECAT7 & VFF & DICOM \\
\hline
\hline
Szolt's programs\footnote{Szolt's programs are several Matlab-scripts, and are not part of any package, and hence the strange name. These scripts can only convert images in certain directions. These are from ECAT6 to VFF, ECAT7 to VFF, Analyze to VFF and VFF to Analyze.} & rw & r & r & rw & - \\
(X)MedCon & rw & rw & r & - & rw \\ 
ImageJ & rw & - & - & - & r \\
MRICro & rw & r & r & r & r \\
SPM & rw & - & - & - & r \\
Vesa's programs & rw & rw & rw & - & - \\
\end{tabular}
\end{center}
\end{minipage}
\caption{The programs available for conversion between different formats. There also exist a program to convert files from a proprietary format used by General Electric in their older MR cameras to the ECAT6 format.}
\label{tbl:convert}
\end{table}


\begin{table}[t]
\begin{minipage}{\textwidth}
\begin{center}
\begin{tabular}{lll}
\multicolumn{3}{c}{ROI drawing programs} \\
\hline
Name &  Projections & Other \\
\hline
\hline
PMOD & All with Java support & Transverse, saggital and coronal\footnote{Requires reloading of the images if the projection is to be changed and viewing ROIs from different projections at the same time is not possible.} \\
Hba &  Transverse, saggital and coronal & Can interpolate missing ROIs\\
ECAT6 & Transverse \\ 
ECAT7 & Transverse \\
MRIcro & Transverse \\
\end{tabular}
\end{center}
\end{minipage}
\caption{Features of the ROI drawing programs.}
\label{draw}
\end{table}

%\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Descriptions of evaluated programs}

Below are descriptions of the programs involved in this evaluation, and their primary use are described along with which kind of license they are under.

\subsection{ImageJ}

\begin{figure}
\centering
\includegraphics[width=10cm]{figures/imagej.eps}
\caption{ImageJ with the scale tool and histogram of the image.}
\label{fig:imagej}
\end{figure}

ImageJ is a simple Java program that works under MacOS 9/X, Windows and Linux (fig. \ref{fig:imagej}). It is mainly developed by National Institutes of Health (NIH), but there are several other contributors and the program lies in the public domain, free to modify by anyone who wants \cite{imagej}.

The default download is in itself not very useful, but it has a library of plugins of considerable size, that makes it simple to extend its functionality. With plugins for reading Analyze images, and more plugins for reordering the image frames and realigning the image, it becomes very versatile, but it lacks some of the more sofisticated tools that are required for the tasks given. Because of this, ImageJ can foremost be seen as a file conversion tool between different resolutions, unwanted slice remover and a simple tool for AC-PC alignment. When used for reslicing, it uses trilinear interpolation.

\subsection{SPM99/2}
\label{spmdesc}

Statistical Parametric Mapping (SPM) is a Matlab program with optimized functions written in C, running under Windows, MacOS and Linux, and it is probably easily portable to other platforms as well. SPM was theoretically conceived and first implemented by Karl Friston. Its main purpose is to do statistical parametric analyzis of PET and fMRI data, but the coregistration and reslicing parts are very useful in themselves. Nowadays the program is maintained and developed by members and collaborators of the Wellcome department of Imaging Neuroscience \cite{spm}. SPM99 was released in 1999 and parts of it is still used in other programs, but SPM2, released 2003, is the latest and most developed version. SPM is released under the GNU General Public License \cite{gpl}.

The program is nonintuitive when first introduced to it, but the design makes more sense as experience with the program increases. One of its quirks is to keep an accompanying \texttt{.mat} file with each Analyze file it has been working with. This file contains a transformation matrix, and information if the image is flipped\footnote{See appendix \ref{anadef} for more information about image orientations.}. This is useful if an image needs to be rotated, and a reslice of the whole volume is unneccessary. When a \texttt{.mat} file is accompanying an image, SPM do not try to figure out if the image shall be flipped, but if it is not there the default setting is to always flip images. This can easily result in an image that is flipped when it should not be, and vice versa. When viewing an image SPM automatically loads the transformation from the \texttt{.mat} file and applies it to the image. This can give the false impression that an image already is aligned or coregistered while it in fact is not.
Problems like the one mentioned above and the general approach of not telling the user much of what is going on gives SPM a steep learning curve, but when the user is familiar with the processes it becomes an advantage, that it does not interrupt the workflow with confirmations for each action taken.
SPM is designed with large studies in mind (e.g. fMRI studies), and if data is present in an uniform way, it is easy to do large batch jobs, which is useful for time consuming tasks such as the coregistration.
It is easy to extend SPM since it is written in Matlab, and for testing new methods and writing small test programs it is by far the most superior program.


\subsection{Hba}

Human Brain Atlas\footnote{This program was originally used to map structures to an atlas, and therefore it is named Human Brain Atlas. It also exist a Monkey Brain Atlas which is identical except for the atlas files.}, or Hba, is a program for drawing ROIs and extracting TACs from the data. It was originally developed by the MR centre at the Karolinska Hospital. The licensing is unknown and altough the source is available, it is not free.

It is a very fragile program which is prone to crashing, especially if it is used a bit differently from what the original developers designed it for. Furthermore, it only works on Irix, and only with a color palette with 256 colors, probably due to optimizations for older SGI graphic cards. At the Karolinska Hospital there are only two SGI machines which can be a bottleneck\footnote{Currently methods to make other X servers use Hba are being tested, with some success. This will remove the bottleneck of only having two physical workstations.}. Hba can only handle VFF files which makes it essential to have tools for converting to this format available if this program is to be used. 

Despite all the drawbacks of this program, it actually is one of the better programs available -- it is fast (even on SGI machines from the late 90's), it is possible to view and draw ROIs in coronal, saggital and transverse projections, it does not require that the PET and MR images have the same modality as long as the voxel sizes are given in the image files, the drawing tools are easy to use. Altogether, it is a program designed to be easy and fast to use.
The final VOI is constructed with Bézier curves with control points from the ROI drawings. This makes it possible to skip slices when drawing -- a ROI for the skipped frame can be interpolated.
In Hba it is possible to draw ROIs in one modality and then applying them to an image in another modality, which in practice results in the possibility to draw on high resolution MR images, while extracting the TACs from lower resolution PET images, which saves both storage space and time, since the processing of large images takes more time. This ability, combined with the ability to interpolate missing ROIs makes it possible to draw on every third MR image with preserved precision, since it is common for MR slices to be $1 mm$ thick, while PET slices usually are $3 mm$.


\subsection{Mars}

Mars is implemented in Matlab, and targetted at a generic x86 Linux machine. It also works under MacOS X, and after the installation of standard GNU tools, also under Windows. It can act as a simplifying layer for AIR and the coregistration part of SPM99. Neither AIR nor SPM99 is platform independent, so a suitable platform with support for these programs must be chosen, if they are going to be used together with Mars \cite{mars}. Mars is being developed as a coregistration part of an EU project called PVEOut, mainly by people at the Copenhagen University Hospital. Licensing information is not available.

The co\-registration methods available in Mars are inter\-active image overlay, intrinsic landmarks (\ref{ana_land}), SPM99 and AIR. The most useful method is the image overlay method for manual coregistration and AC-PC alignment. The SPM99 and AIR coregistration algorithms are inferiour to the mutual information method and are therefore not of any use. Currently Mars is of limited use since it does not provide anything that other packages cannot solve, but it is still under development and it may be useful in a later incarnation.

\subsection{AIR}
\label{airdesc}

AIR consist of a bundle of command line tools written by Roger P. Woods. It is written in ANSI C and compiles on most platforms. To obtain a copy of the source code\footnote{The AIR package is only available as source code.} it is required that the downloader belong to the academia. The sole purpose of AIR is to coregister images.

The user interface, or rather lack of, makes AIR very hard to use, especially since it depends on many parameters. The program requires the MR images to be skull-stripped prior to registration.

\subsection{Vesa's programs}

These are conversion programs for most formats except DICOM and VFF. They are written by Vesa Oikonen at the PET centre in Turku, Helsinki and are available under the GNU Library General Public License \cite{lgpl}. They are command line based, but for the given task it is adequate. The programs should compile on most Unices, and has been tested under Linux and MacOS.

\subsection{ECAT6}

ECAT6 is a program developed by Siemens and it is made especially for use together with their PET cameras. It runs only under Solaris, and is a proprietary program with its own proprietary format. Due to time constraints, it has not been thoroughly tested. ECAT6 can be used to draw ROIs and extract TACs, and has been heavily used earlier.

\subsection{ECAT7}

ECAT7 is a newer version of ECAT6, and like ECAT6, it only runs under Solaris. Currently, ECAT7 is not installed on any workstation, and is not tested.

\subsection{PMOD}
\label{pmoddesc}

\begin{figure}
\centering
\includegraphics[width=8cm]{figures/pmod_fusion.eps}
\caption{PMOD fusion tool in action with MR image top left and and summation PET image to the top right. As can be seen in the lower right window where the images are overlayed they are not yet aligned.}
\label{pmodfusion}
\end{figure}


PMOD is a commersial program written in Java, it has been tested under MacOS X, Linux and Windows. The licensing fees are variable depending on the number of licenses and if the buyer is an academic institution or if a corporate license is required. A single license for academic use with the modules required for work at the Karolinska Hospital costs US \$4800 \cite{pmod}.

PMOD is a program that is made to be able to do all the steps necessary in PET image analyzis. It is built as a set of different modules, of which each can do a certain step. 
It can coregister images with mutual information, the SPM99 method, least squares and manual overlay methods.
It is possible to draw ROIs and the tools for the ROI drawing are good. Directly after the ROI drawing, it is possible to extract the statistical data and do kinetic analyzis of the data in another PMOD module. The integration of the different modules makes the workflow smooth and the user is not disrupted by minor annoyances such as file format conversion between the different steps.

PMOD is easy to learn, with easily accessible options (see fig. \ref{pmodfusion}). Furthermore it can load images in most formats directly, without any conversion tools. Although it is important to remember that it assumes that the radioactivity is measured in $kBq/cc$ if it is not specified in the image header. It is possible to change metric during the loading phase, but this will be converted to $kBq/cc$. When an image is saved, it will be saved with intensity values in $kBq/cc$ regardless of the original format, which can be confusing.

The GUI relies heavily on the use of a mouse and it cannot be scripted which is unfortunate when batch processing of images is needed. It is also slow when displaying images, especially if interpolation of the image is turned on.

A drawback for PMOD that is easy to overlook is that it is a closed source program. If anything is broken, it will stay broken until a patch is available. Options that would be usable cannot be added, but has to be requested and hopefully included in the next version.
Further drawbacks of its closed source nature are that it is not possible to see how different methods are implemented, the only way to find out is to mail the developers and ask them. Even then, it is not possible to see how the low-level details are handled, for example if they introduce some randomness into the mutual information coregistration process to reduce artefacts \cite{tsao2003}.
As mentioned earlier, PMOD is written in Java, with the effect that it in our experience requires at least 1.5GB of RAM to be useful, and preferrably even more.
Although it seems that Java cannot make use of more than approximately 1.9GB memory for heap allocation which makes it neccessary to restart PMOD after heavy usage, since there is no way to tell the Java Virtual Machine to run the garbage collector.
The choice of language could also be a possible reason for the experienced slowness, but the PMOD team says that the next release will focus on solving the memory and speed issues.

PMOD can only do trilinear interpolation when reslicing, with no option to change it.

\subsection{(X)MedCon}
(X)MedCon is a conversion program written by Erik Nolf, and is available under the GNU General Public License \cite{gpl,xmedcon}.

It can load and save in the most commonly used formats, and is easy to use. Unfortunately, the produced images are useless because it distorts the voxel values when converting from ECAT6/7 to Analyze, which is the conversion that is most useful. This reduces the uses of the program to only viewing images, a task that is possible in other programs.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Alignment and Head Movement}


There are two aspects when trying to minimize error caused by head movement. The first is to minimize head movement during the scan, and the second is to try to correct the errors which has occured afterwards. The ``afterward'' approach can be done by comparing successive frames and aligning them internally. Another method is to register the head movement during the scan and adjusting the images accordingly.

\begin{figure}
\centering
\includegraphics[width=9cm]{figures/helmet.eps}
\caption{Plastic helmet made for each patient to fixate the head position during the PET scan.}
\label{fig:helmet}
\end{figure}
\label{helmet}
To keep the subject from moving the head a special plastic helmet is made, individually for each subject \cite{bergstroem}. The helmet is fixated in the correct position by using dental markers molded to fit each subjects teeth (see fig. \ref{fig:helmet}). This helmet is also supposed to keep the head in the same position in both the PET camera and the MR camera. Unfortunately, it seems like the device keeps the head tilted around 15 degrees differently in the transaxial dimension. This is easily corrected later, but could be avoided with better calibration of the gear that holds it in place.

No method for registering head movement during the scan is currently in use, and no coregistration is done for each individual timeframe. This is because the possible movement when wearing the helmet is believed to be negligible.

For alignment, especially to the AC-PC line, SPM, ImageJ, PMOD and Mars all do a good job. Easiest to use is Mars.


\section{Reslicing}
\label{reslicing}
Reslicing occurs when a transform model is applied to an image, which requires the image to be interpolated.

The programs support different methods for the interpolation step. To evaluate some of them, the nearest neighbour, trilinear and spline interpolation methods has been tested by rotating a volumetric image in the coronal plane. It was rotated $10$ times with a step size of $\frac{2\pi}{10}$ radians. With this step size the resulting interpolated image was orientated as the original, and ROIs defined on the original could be reused without introducing any error by rotating the ROIs or by having a different set of ROIs. No translation was done. The time activity curves for the resulting images can be seen in appendix \ref{app:interp}, in figures \ref{fig:intcer}, \ref{fig:intlput}, \ref{fig:intrput}, \ref{fig:intlstria}, \ref{fig:intrstria}, \ref{fig:intlcau} and \ref{fig:intrcau}. By observing that the activity for the different interpolation methods tend to stay internally ordered under the ellapsed time, it has been assumed that it is possible to sum all values into a single value. By subtracting these values from the original data, a relative error has been calculated which can be seen in fig. \ref{fig:interror}.

The radioligands for the subjects are listed in table \ref{tbl:ligands}.

Volumes of the different regions are listed in table \ref{intvol}.

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/interp_error/anqv.eps}
\includegraphics[width=6cm]{figures/interp_error/erkr.eps}
\includegraphics[width=6cm]{figures/interp_error/olgi.eps}
\includegraphics[width=6cm]{figures/interp_error/olma.eps}
\caption{Interpolation error (lower is better, execution time increases to the right). The total amount of activity has been calculated for each of the 7 regions, after being rotated by the different interpolation methods. These values have then been subtracted from the total activity in the original data, and finally normalized. These values are plotted against the normalized interpolation execution time. The dashed line is only a hint to show which data points that come from the same region and do not contain any information.}
\label{fig:interror}
\end{figure}

\begin{table}
\begin{center}
\begin{tabular}{lll}
\multicolumn{3}{c}{Radioligands used} \\
\hline
Subject & Radioligand & Traces \\
\hline
\hline
1 & $[^{11}C]PE2I$ & Dopamine transport, DAT \\
2 & $[^{11}C]MADAM$ & 5HT transport, SET \\
3 & $[^{11}C]SCH23390$ & D1 receptor \\
4 & $[^{11}C]raclopride$ &  D2 receptor \\
\end{tabular}
\end{center}
\caption{The PET images in this text where made with these radioligands.}
\label{tbl:ligands}
\end{table}

\begin{table}
\begin{center}
\begin{tabular}{llrr}
\multicolumn{4}{c}{ROI Volumes} \\
\hline
&& \multicolumn{2}{c}{Volume ($mm^3$)}\\
\hline
Region & Subject & Hba & PMOD \\
\hline
\hline
Cerebellum & 1 & 18535.9 & 13252.6	\\
Cerebellum & 2 & 15146.4 & 18467.8	\\
Cerebellum & 3 & 19857.4 & 13252.6	\\
Cerebellum & 4 & 22952.4 & 10466.7	\\
Left Caudate Nucleus & 1 & 3148.11 & 3018.9	\\
Left Caudate Nucleus & 2 & 3820.83 & 3716.4	\\
Left Caudate Nucleus & 3 & 4300.15 & 3627.7	\\
Left Caudate Nucleus & 4 & 4053.64 & 4179.6	\\
Left Putamen & 1 & 2299.03 & 2546.1	\\
Left Putamen & 2 & 3848.22 & 3188.3	\\
Left Putamen & 3 & 3560.63 & 3627.7	\\
Left Putamen & 4 & 3327.82 & 3398.7	\\
Left Ventral Striatum & 1 & 783.761 & 403.5	\\
Left Ventral Striatum & 2 & 931.242 & 721.0	\\
Left Ventral Striatum & 3 & 1122.97 & 699.3	\\
Left Ventral Striatum & 4 & 821.684 & 442.5	\\
Left Dorsolateral Prefrontal Cortex & 1 & - & 22642.9	\\
Left Dorsolateral Prefrontal Cortex & 2 & - & 21963.6	\\
Left Dorsolateral Prefrontal Cortex & 3 & - & 24032.7	\\
Left Dorsolateral Prefrontal Cortex & 4 & - & 20612.3	\\
Right Caudate Nucleus & 1 & 3213.42 & 2845.6	\\
Right Caudate Nucleus & 2 & 3588.02 & 2727.2	\\
Right Caudate Nucleus & 3 & 3985.17 & 3449.3	\\
Right Caudate Nucleus & 4 & 3861.92 & 3740.1	\\
Right Putamen & 1 & 2050.84 & 1987.8	\\
Right Putamen & 2 & 3451.07 & 2931.5	\\
Right Putamen & 3 & 3861.92 & 3618.4	\\
Right Putamen & 4 & 3519.55 & 3139.8	\\
Right Ventral Striatum & 1 & 509.444 & 325.5	\\
Right Ventral Striatum & 2 & 890.158 & 479.6	\\
Right Ventral Striatum & 3 & 821.684 & 633.3	\\
Right Ventral Striatum & 4 & 657.347 & 366.8	\\
Right Dorsolateral Prefrontal Cortex & 1 & - & 21932.4	\\
Right Dorsolateral Prefrontal Cortex & 2 & - & 21670.6	\\
Right Dorsolateral Prefrontal Cortex & 3 & - & 24213.2	\\
Right Dorsolateral Prefrontal Cortex & 4 & - & 19897.5	\\
\end{tabular}
\end{center}
\caption{The volumes of the different regions during the tests. The volumes for the regions in the interpolation tests are the Hba volumes.}
\label{intvol}
\end{table}


From the tests we can assume that $\beta^2$-spline interpolation performs as good as the splines of higher degree, but at less than half the time as the $\beta^7$-spline. Both trilinear and nearest neighbour interpolation performs worse in most cases, but sometimes one or both of them performs better than the spline methods. Most interesting is that trilinear interpolation overall is the worst performer. 



\section{Coregistration}

If it existed an universal method to measure how good the coregistrations were done, this method should be used as the similarity measure in a new coregistration method. This method would then become the best since it can always make the best coregistrations.
Since no such method exist, the methods has to be compared to a given ``golden standard''. In this case the golden standard is an occular inspection by the examiner, who can use the overlay method (\ref{overlay}) for comparisons. It is in no way said that image overlay is the best method to inspect coregistration results, but there are no other usable metrics known to the author. The results of an ordinary coregistration is shown in figure \ref{fig:coreg}.

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/before_coreg.eps}
\includegraphics[width=6cm]{figures/after_coreg.eps}
\caption{This is a typical example of an coregistration, this time done with the MI method. As can be seen, the difference in not that large and small movements in either direction can be hard to classify as good or bad.}
\label{fig:coreg}
\end{figure}

Several methods has been tested, but none could surpass mutual information coregistration.
Both SPM2 and PMOD uses Mutual Information for coregistration, and to see if there is any difference between these two, an extreme case has been tested since they performed similarly under usual conditions. A PET image was translated, rotated and cut in half. This is a case which simulates an atrophy that probably is mortal, and a PET scan in this condition is most unlikely, but it pushes the coregistration methods to their limits. The results are shown in fig \ref{fig:misalignedcoreg}. Surprisingly, PMOD fails to register the PET image that is cut in half, even though it is using the same method as SPM2. Under normal circumstances PMOD and SPM2 producre very similar images.

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/half_brain.eps}
\includegraphics[width=6cm]{figures/half_brain_pmod.eps}
\includegraphics[width=6cm]{figures/half_brain_spm.eps}
\includegraphics[width=6cm]{figures/half_brain_air.eps}
\caption{The top left image is the original data -- severely misaligned and cut in half, which was given to AIR, SPM2 and PMOD. The top right image is the MI coregistration by PMOD. It is interesting to see how it actually identifies cortex, but places it in the wrong place. The bottom left is the result from SPM2. The bottom right image is done with AIR, but without skull stripping. It was not able to perform a coregistration.}
\label{fig:misalignedcoreg}
\end{figure}


\section{ROI drawing}

The time activity curves from ROIs drawn in different programs and different brain alignments can be seen in appendix \ref{app:roi}, figures \ref{fig:clput} -- \ref{fig:crdlpfc}. No data is available for dorsolateral prefrontal cortex with Hba because of the authors inability to draw them without crashing the program.

\chapter{Results}
\label{results}

In this chapter the methods that are preferrable will be discussed, and why they are viewed as the best choice at the time being. If a method is described as fast, it is completed well below one second and does not delay and anger the user.

\section{AC-PC Alignment}

All programs that can do this do it equally good, since it is a manual procedure. The programs that can be used are Mars, PMOD, SPM and ImageJ. There are two parameters that affects the choice here -- usability and interpolation method.

Mars is by far the easiest program to use for this due to its quickness and large view of the images. The quickness is partly due to the convinient keyboard short-cuts it uses to rotate and translate the image. It also shows the selected projection in full screen mode which makes it easy to discern small differences, and contrast settings are easily modified.

PMOD can do alignment but for unknown reasons it is incredibly slow. It is easy to use but it relies on the use of a mouse and requires the operator to click on small buttons repeatedly to do small controlled transformations. For large studies, the operator will be fatigued. It takes on average 10 seconds for a single rotation/translation step regardless of its size on a Pentium4 at 2.8GHz with 2GB memory at 800MHz. 

SPM can do AC-PC alignment, even though the original design did not take this into consideration. When manually entering in the transormation matrix values it shows the result and the transformation matrix can then be saved and used for alignment. It is fast, but not particularly user friendly. The transformation can be saved as parameters only in a separate file, which is applied later which reduces the number of interpolations of the image data.

ImageJ can with the aid of plugins do alignment, and even though it is written in Java just as PMOD, it is considerably faster. 


Overall, the choice of program if AC-PC alignment is needed is best decided by user preferences, or availability. There is no need to use a particular one because of some superior feature.


\section{Coregistration}

The manual coregistration procedures are very good for validation purposes, but not suitable for use in real studies since they are suspectible to operator bias and do not work in a consistent way and they are extremely time consuming.

The AIR coregistration method requires too much preprocessing of the images to save any time compared to the manual methods. The only advantage then is that the method used for coregistration will be the same for all subjects and will not depend on the user. AIR is not recommended for use because of the cumbersome command line interface, its sub-optimal performance on images poorly aligned prior to coregistration, the need for MR images to be skull-stripped and the inability to view images.

Mars can only do coregistration manually, and even though it makes this easy, it suffers from the drawbacks for manual methods mentioned above. 

SPM2 can coregister with the mutual information algorithm, and the image in fig. \ref{fig:misalignedcoreg} was calculated in 140 seconds. The program asks for the number of subjects and batch processing of an arbitrary number of images is possible.

PMOD can use the same method as the ones used in AIR and in SPM99, but mutual information is the default and is the method with the best results. It took PMOD 110 seconds to calculate the image in fig. \ref{fig:misalignedcoreg}, which is faster than SPM2 but then the image loading and resampling times of the images were not included. The lack of batch processing or scripting support makes the coregistration process in PMOD slower overall, since the user is required to repeat the process for each subject. PMOD will eventually run out of memory when this task is repeated, which is probably due to the nature of Java where there is no method for the programmer to free unused memory, and the garbage collector is not doing what it should fast enough.

To conclude, concerning the algorithm, the mutual information method is superiour to all other methods and should be used.
The best program for coregistration from a user perspective is SPM2 since it supports batch processing, and SPM2 is also the best performer, at least in extreme cases as in fig. \ref{fig:misalignedcoreg}. PMOD perform similarly to SPM2, and the choice is between these two, since they are the only ones that use the mutual information algorithm.

\section{Interpolation}

Nearest neighbour is not widely available and is only available in PMOD as an interpolant for the automatic coregistration algorithms to make them run faster, because it is regarded as an inferior method for final data.

Trilinear interpolation is the most common method, and can be used in all tested programs. 

Spline interpolation is currently only available in SPM2.

Sinc interpolation has some issues since it cannot have infinite support, which could be solved by carefully choosing a window function. But choosing this is not trivial and sinc interpolation is not neccessarily best when it is not ideal cases of information that is to be represented. Because of the problems associated with sinc interpolation, it is not in common use, and has not been tested.

In the tests performed in section \ref{reslicing} the interpolation algorithm with least error where the spline algorithm of degree 2 and higher. Internally for the spline algorithms, they performed very similar. Therefore the best choice, if calculation times are to be considered, is the $\beta^2$-spline interpolation. This claim is strengthened by the results in \cite{meijering2001}, although the results presented there are synthetic and not tested on real data as is the case in this test.

A surprising result is that nearest neighbour interpolation in most cases perform better than trilinear interpolation, something that was not theoretically expected. A possible reason is that ROIs in most cases contain higher pixel values than surrounding data, and as trilinear interpolation aquires the interpolated data linearly from all neighbours, the data spills out from the ROI to the surrounding area. Nearest neighbour interpolation has a very abrupt cutoff instead which sets the interpolated voxel values either to a value from within the VOI or to one from outside the VOI, which over large numbers of pixels may result in an average that is closer to the original values.

The errors in the tests in this text are exaggerated since the images where interpolated 10 times. Although, due to the smoothing effect the error cannot be said to be 10 times smaller when only one interpolation is performed\footnote{This is not true for nearest neighbour interpolation, since it is not an intensity interpolator.}. The images where coregistered beforehand with trilinear interpolation when reslicing, which makes them slightly presmoothed, which further reduces the 10 times factor.

In conclusion, of the tested algorithms, any high degree spline interpolation algorithm is the best choice, with the surprising result that nearest neighbour interpolation is the second best. Trilinear interpolation shows the least desirable results.

\section{Drawing Regions of Interest}

The ROIs can take different shapes in different programs and on different computers for several reasons. Simple things as the contrast and $\gamma$-ramp can change the position of perceived contours. The light conditions and type of monitor can also be a factor. In an attempt to see if the differences beteeen programs are large, the same ROIs has been drawn in both Hba and PMOD, and to see if AC-PC alignment has any effect, both AC-PC aligned and unaligned images has been used for ROI drawing in PMOD. The result of this is presented in figures \ref{fig:clput} -- \ref{fig:crdlpfc}.

From the results it can definitely be concluded that AC-PC alignment does not have any influence on the resulting values for dorsolateral prefrontal cortex, the only region of the drawn regions that theoretically should be affected by the AC-PC alignment. The curves are very similar in all cases except for the left dorsolateral prefrontal cortex in subject 1, but the difference is still smaller than for other areas.

Regarding the other regions, most results are similar with some fluctuations. The ROIs are constant through time but the TACs are not particularly smooth, which suggests that the errors are small enough to be influenced by noise in the image data. This suggests that the choice of program do not significantly influence the final result.

Subject 1 has the most varying result. The reason for this could possibly be explained by looking at the radioligand used, $PE2I$, which has a very low specific/non-specific affinity ratio which makes it hard to use.

This leaves the choice of ROI drawing program to the user. Hba is definitely the fastest of them all, especially if its ability to interpolate between slices is exploited. PMOD is very good with a user friendly user interface, with the drawback that all slices need to be drawn upon.

%ECAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAT6


\section{Conclusion}

There seem to be two viable paths for future studies. One is to use PMOD for all required steps, and the other is to use the combination of SPM2 and Hba.

The advantages with PMOD is the reduced amount of time spent handling files and converting them and its ease of use to get started with. A downside is, as mentioned earlier, it is a commersial program with closed source, which makes it impossible to correct bugs and look at its internal workings. PMOD is probably the program of choice for studies of moderate size where the lack of scriptability is compensated by its ease of use, and when there are no doubts regarding the methods.

Advantages with SPM2 and Hba over PMOD are that they are faster to work with, they are free and generally seems to have a better design when it comes to large studies, even though they are not the easiest to start with. 
They still require some glue between them in the form of file conversion programs and the frequent crashings of Hba is an issue. But since they are separate programs, Hba can easily be replaced by any other ROI drawing program and there is no risk of vendor lock-in as is the case with PMOD. For development work SPM2 is clearly superiour.

It is not necessary to do AC-PC alignment, since no difference where discernible in the areas that theoretically should be affected. Other smaller areas where more affected, which suggests that ROI definition variations and image noise are factors that should be reduced first.

Interpolation should be done with splines of degree 2 or higher, which only SPM2 can do. If spline interpolation is not available, nearest neighbour interpolation should be used for quantitative studies, while trilinear interpolation can be used for viewing purposes, since it looks better than nearest neighbour.


\chapter{Discussion}

\section{Limitations}

No rigorous statistical analysis has been done on the resulting data. This is because all the PET images used where done with different ligands, and since different ligands behave in different ways, it would not be meaningful. To produce any useful statistical data, more images produced with the same ligand must be used.

With the knowledge that interpolation affects real PET images differently from test images, more interpolation algorithms could have been included.

\section{Further work}

To truly see why the interpolated values differ the way they do, a possibility would be to draw ROIs around the original ROIs. With this data it would be possible to see how intensity values are influenced by surrounding intensities and how these values spills in or out.

It would be interesting to see how much the interpolated values would deviate from the original with a smaller number of interpolations. Only one is not possible, since the ROIs drawn on the original cannot be reused. A good number would be two.

The interpolation, AC-PC alignment and coregistration all serve to make the BP values more reliable, but those have not been calculated. To validate the results for final BP values, this must be done.

The file format Analyze which is in widespread use has several disadvantages. It does not contain all the information that is required to do a study and must therefore be complemented with more files containing the missing information. There exist a new Analyze format which may solve this, and other formats such as MINC \cite{minc}. Developing a new format is tempting, but then the gargantuan task to make all developers accept this as the \textit{de facto} file format standard remains.


\chapter{Acknowledgments}

I would like to thank my supervisor at the Karolinska Hospital, Per Karlsson, for taking the time to teach me how PET works, and for his help in showing me how the current programs work. Stefan Pauli for explaining much of the physics that makes PET a possibility. My supervisor at Umeå University -- Fredrik Georgsson for giving me help and feedback for this report.
Jacqueline Borg for setting me in contact with the PET group. Lars Farde for making my stay at the Karolinska Hospital a possibility. Vesa Oikonen for adding features and fixing bugs in his programs almost before I had them reported. The PMOD team for taking the time to answer my questions. 
\appendix

\chapter[Interpolation results]{Interpolation data for different interpolation algorithms}
\label{app:interp}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac/cer_anqv.eps}
\includegraphics[width=6cm]{figures/tac/cer_erkr.eps}
\includegraphics[width=6cm]{figures/tac/cer_olgi.eps}
\includegraphics[width=6cm]{figures/tac/cer_olma.eps}
\caption{Interpolation of cerebellum in four subjects}
\label{fig:intcer}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac/l_dput_anqv.eps}
\includegraphics[width=6cm]{figures/tac/l_dput_erkr.eps}
\includegraphics[width=6cm]{figures/tac/l_dput_olgi.eps}
\includegraphics[width=6cm]{figures/tac/l_dput_olma.eps}
\caption{Interpolation of left putamen in four subjects}
\label{fig:intlput}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac/r_dput_anqv.eps}
\includegraphics[width=6cm]{figures/tac/r_dput_erkr.eps}
\includegraphics[width=6cm]{figures/tac/r_dput_olgi.eps}
\includegraphics[width=6cm]{figures/tac/r_dput_olma.eps}
\caption{Interpolation of right putamen in four subjects}
\label{fig:intrput}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac/l_vstria_anqv.eps}
\includegraphics[width=6cm]{figures/tac/l_vstria_erkr.eps}
\includegraphics[width=6cm]{figures/tac/l_vstria_olgi.eps}
\includegraphics[width=6cm]{figures/tac/l_vstria_olma.eps}
\caption{Interpolation of left striatum in four subjects}
\label{fig:intlstria}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac/r_vstria_anqv.eps}
\includegraphics[width=6cm]{figures/tac/r_vstria_erkr.eps}
\includegraphics[width=6cm]{figures/tac/r_vstria_olgi.eps}
\includegraphics[width=6cm]{figures/tac/r_vstria_olma.eps}
\caption{Interpolation of right striatum in four subjects}
\label{fig:intrstria}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac/l_cau_anqv.eps}
\includegraphics[width=6cm]{figures/tac/l_cau_erkr.eps}
\includegraphics[width=6cm]{figures/tac/l_cau_olgi.eps}
\includegraphics[width=6cm]{figures/tac/l_cau_olma.eps}
\caption{Interpolation of left caudate nucleus in four subjects}
\label{fig:intlcau}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac/r_cau_anqv.eps}
\includegraphics[width=6cm]{figures/tac/r_cau_erkr.eps}
\includegraphics[width=6cm]{figures/tac/r_cau_olgi.eps}
\includegraphics[width=6cm]{figures/tac/r_cau_olma.eps}
\caption{Interpolation of right caudate nucleus in four subjects}
\label{fig:intrcau}
\end{figure}


\chapter{ROI data for different programs}
\label{app:roi}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac_comp/l_dput_anqv.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_dput_erkr.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_dput_olgi.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_dput_olma.eps}
\caption{Comparison of time activity curves for left putamen in four subjects with different programs.}
\label{fig:clput}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac_comp/r_dput_anqv.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_dput_erkr.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_dput_olgi.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_dput_olma.eps}
\caption{Comparison of time activity curves for right putamen in four subjects with different programs.}
\label{fig:crput}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac_comp/l_vstria_anqv.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_vstria_erkr.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_vstria_olgi.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_vstria_olma.eps}
\caption{Comparison of time activity curves for left striatum in four subjects with different programs.}
\label{fig:clstr}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac_comp/r_vstria_anqv.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_vstria_erkr.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_vstria_olgi.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_vstria_olma.eps}
\caption{Comparison of time activity curves for right striatum in four subjects with different programs.}
\label{fig:crstr}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac_comp/l_cau_anqv.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_cau_erkr.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_cau_olgi.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_cau_olma.eps}
\caption{Comparison of time activity curves for left caudate nucleus in four subjects with different programs.}
\label{fig:clcau}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac_comp/r_cau_anqv.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_cau_erkr.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_cau_olgi.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_cau_olma.eps}
\caption{Comparison of time activity curves for right caudate nucleus in four subjects with different programs.}
\label{fig:crcau}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac_comp/l_dlpfc_anqv.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_dlpfc_erkr.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_dlpfc_olgi.eps}
\includegraphics[width=6cm]{figures/tac_comp/l_dlpfc_olma.eps}
\caption{Comparison of time activity curves for left dorsolateral prefrontal cortex in four subjects with different programs.}
\label{fig:cldlpfc}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=6cm]{figures/tac_comp/r_dlpfc_anqv.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_dlpfc_erkr.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_dlpfc_olgi.eps}
\includegraphics[width=6cm]{figures/tac_comp/r_dlpfc_olma.eps}
\caption{Comparison of time activity curves for right dorsolateral prefrontal cortex in four subjects with different programs.}
\label{fig:crdlpfc}
\end{figure}


\chapter[Anatomical Definitions]{Anatomical Definitions and Image Orientations}
\label{anadef}

An image, either MR or PET, is by convention viewed ``from below'' in the \textit{transverse} plane (fig. \ref{projs}).

The side where the spine is located is the \textit{dorsal} side, while the opposite side is called \textit{ventral}. In human anatomy, the \textit{anterior} side is the front with nipples and navel and is equivalent to the \textit{ventral} side, while the back is called the \textit{posterior} side, which in humans is equivalent with the \textit{dorsal} side.
Structures close to the midline is called \textit{medial}, and towards the sides they are called \textit{lateral}.
When viewed \textit{coronally}, the image can either be flipped or not. If it is flipped, the right side is right on the patient too (the image is looking at us). Unflipped images are the other way around. Radiologists tend to prefer unflipped images while neurologists prefer them flipped. This is unfortunate since it is not obvious if the image is flipped or not.

In fig. \ref{mr_pet_voi} and \ref{dlpfc_cer} the anatomical regions used in this text are marked.


\begin{figure}
\centering
\includegraphics[width=7cm, height=5cm]{figures/jore_saggital_124.eps}
\includegraphics[width=7cm]{figures/jore_transverse_128.eps}
\includegraphics[width=7cm]{figures/jore_coronal_124.eps}
\caption{Projections commonly used. Top: saggital, middle: transverse, bottom: coronal.}
\label{projs}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=8cm]{figures/coronal_voi_mr.eps}
\includegraphics[width=8cm]{figures/coronal_voi_pet.eps}
\caption{ROIs first drawn coronally on an MR image and then translated to a PET image in PMOD. A: right putamen, B: right dorsal caudate nucleus, C: right ventral striatum, D: left caudate nucleus, E: left putamen, F left ventral striatum.}
\label{mr_pet_voi}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=8cm]{figures/dlpfc.eps}
\includegraphics[width=8cm]{figures/cer.eps}
\caption{Top: left and right dorsolateral prefrontal cortex, bottom: cerebellum.}
\label{dlpfc_cer}
\end{figure}

\bibliography{report}

%\printindex

\end{document}

% Puh, det var jobbigt att skriva allt detta. Men nu är det klart!


