\chapter{Concepts}\index{Concepts}%
\setheader{{\it CHAPTER \thechapter}}{}{}{}{}{{\it CHAPTER \thechapter}}%
\setfooter{\thepage}{}{}{}{}{\thepage}%

\section{Overview}\label{conceptoverview}\index{Concepts,Overview}%
In CTSim, a phantom object, or a geometrical description of the object
of a CT study is constructed and an image can be created.  Then a
scanner geometry can be specified, and the projection data simulated.
Finally that projection data can be reconstructed using various user
controlled algorithms producing an image of the phantom or study object.

In order to use CTSim effectively, some knowledge of how CTSim works
and the approach taken is required. \ctsim deals with a variety of
object, but the two we need to be concerned with are the 'phantom' and
the 'scanner'.

\section{Phantoms}\label{conceptphantom}\index{Concepts,Phantoms}%
\subsection{Overview}\label{phantomoverview}\index{Concepts,Phantoms,Overview}%

CTSim uses geometrical objects to
describe the object being scanned: rectangles, triangles, ellipses,
sectors and segments.  With these the standard phantoms used in the CT
literature (the Herman and the Shepp-Logan) can be constructed.  In fact
CTSim provides a shortcut to construct those phantoms for you.  It also 
allows you to write a file in which the composition of your own phantom is
described.

The types of phantom elements and their definitions are taken from Herman's 1980 
book\cite{HERMAN80}.

\subsection{Phantom File}\label{phantomfile}\index{Concepts,Phantoms,File}
Each line in the text file describes an element of the 
phantom.  Each line contains seven entries, in the following form:
\begin{verbatim}
item cx cy dx dy r a
\end{verbatim}
The first entry defines the type of the element, one
of {\tt rectangle}, {\tt ellipse}, {\tt triangle}, {\tt sector}, or {\tt segment}.  
{\tt cx}, {\tt cy}, {\tt dx} and {\tt dy} have different meanings depending on the element type.

{\tt r} is the rotation applied to the object in degrees counterclockwise, 
and {\tt a} is the X-ray attenuation coefficient of the object.
Where objects overlap, the attenuations of the overlapped objects are summed.



\subsection{Phantom Elements}\label{phantomelements}\index{Concepts,Phantoms,Elements}

\subsubsection{ellipse}
Ellipses use dx and dy to define the semi-major and semi-minor axis lengths,
with the centre of the ellipse at cx and cy.  Of note, the commonly used
phantom described by Shepp and Logan\cite{SHEPP77} uses only ellipses.

\subsubsection{rectangle}
Rectangles use 
cx  and cy to define the position of the centre of the rectangle with respect
to the origin.  dx and dy  are the half-width and half-height of the
rectangle. 

\subsubsection{triangle}
Triangles are drawn with the centre of the base at cx,cy, with a base
width of 2*dx in x direction, and a height of dy.  Rotations are then
applied about the origin. 

\subsubsection{sector}
It appears that dx and dy
define the end points of a radius of the sector, from which the radius and
the angle of the two arms of the sector are calculated.  But then
orientation and centreing of the sector don't make much sense yet.  

\subsubsection{segment}
Segments are the segments of a circle between a chord and the
perimeter of the circle.  This also isn't clear to me, but it appears that 
perhaps the distance from chord to circle perimeter, and circle radius is
defined by dx and dy. Chord is always horizontal through the origin, then
translated and then rotated (???).

\subsection{Phantom Size}
Also note that the overall dimensions of the phantom are increased by 1\%
above the specified sizes to avoid clipping due to round-off errors.  If the phantom is defined as
a rectangle of size 0.1 by 0.1, the actual phantom has extent $\pm$0.101 in
each direction.

\section{Scanner}\label{conceptscanner}\index{Concepts,Scanner}%
\subsection{Geometries}
This is where things get tricky.  There are two possible approaches.  The
simple approach would be to define the size of a phantom which is put at
the centre of the scanner.  The scanner would have it's bore size defined,
or perhaps better, the field of view defined. Here, field of view would be
the radius or diameter of the circular area from which data is collected
and an image reconstructed.  In a real CT scanner, if the object being
scanned is larger than the field of view, you get image artifacts.  And of
course you can't stuff an object into a scanner if the object is larger
than the bore!  In this model, the scanner size or field of view would 
be used as the standard length scale.

However, CTSim takes another approach.  I believe this approach arose
because the "image" of the phantom produced from the phantom description
was being matched to the reconstruction image of the phantom.  That is, 
the dimensions of the 'before' and 'after' images  were being matched.
The code has a Phantom object and a Scanner object.  The geometry of the
Scanner is defined in part by the properties of the Phantom.  In fact,
all dimensions are determined in terms of the phantom size, which is used
as the standard length scale.     Remember, as mentioned above, the
phantom dimensions are also padded by 1\%.

The maximum of the phantom length and height is used as the phantom
dimension, and one can think of a square bounding box of this size 
which completely contains the phantom.  Let $l_p$ be the width (or height)
of this square. 

\subsubsection{Focal Length \& Field of View}
The two other important variables are the field-of-view-ratio ($f_{vR}$) 
and the focal-length-ratio ($f_{lR}$).  These are used along with $l_p$ to
define the focal length and the field of view (not ratios) according to
\begin{equation}
f_l = \sqrt{2} (l_p/2)(f_{lR})= (l_p/\sqrt{2}) f_{lR}
\end{equation}
\begin{equation}
f_v = \sqrt{2}l_p f_{vR}
\end{equation}
So the field of view ratio is specified in units of the phantom diameter,
whereas the focal length is specified in units of the phantom radius.  The 
factor of $\sqrt(2)$ can be understood if one refers to figure 1, where
we consider the case of a first generation parallel beam CT scanner.

\subsubsection{Parallel Geometry}\label{geometryparallel}\index{Concepts,Scanner,Geometries,Parallel}
\begin{figure}
\includegraphics[width=\textwidth]{ctsimfig1.eps}
\caption{Geometry used for a 1st generation, parallel beam CT scanner.}
\end{figure}

In figure 1A, the excursion of the source and detector need only be $l_p$,
the height (or width) of the phantom's bounding square. However, if the
field of view were only $l_p$, then the projection shown in figure 1B
would clip the corners of the phantom.  By increasing the field of view by
$\sqrt{2}$ the whole phantom is included in every projection.  Of course,
if the field-of-view ratio $f_{vR}$ is larger than 1, there is no problem.
However, if $f_{vR}$ is less than one and thus the scanner is smaller than
the phantom, then distortions will occur without warning from the program.

The code also sets the detector length equal to the field of view in this
case.  The focal length is chosen to be $\sqrt{2}l_p$ so the phantom will 
fit between the source and detector at all rotation angles, when the focal
length ratio is specified as 1.  Again, what happens if the focal length
ratio is chosen less than 1?

The other thing to note is that in this code the detector array is set to
be the same size as the field-of-view $f_v$, equation (2).  So, one has to 
know the size of the phantom to specify a given scanner geometry with a 
given source-detector distance (or $f_l$ here) and a given range of
excursion ($f_v$ here).  

\subsubsection{Divergent Geometries}\label{geometrydivergent}\index{Concepts,Scanner,Geometries,Divergent}
Next consider the case of equilinear (second generation) and equiangular 
(third, fourth, and fifth generation) geometries.  
The parts of the code  relevant to this
discussion are the same for both modes.  In the equilinear mode, a single 
source produces a fan beam which is read by a linear array of detectors.  If
the detectors occupy an arc of a circle, then the geometry is equiangular.
See figure 2. 
\begin{figure}
\includegraphics[width=\textwidth]{ctsimfig2.eps}
\caption{Equilinear and equiangular geometries.}
\end{figure}

For these geometries, the following logic is executed:  A variable dHalfSquare
$d_{hs}$ is defined as
\begin{equation}
d_{hs} = (f_v)/(2\sqrt{2}) = (l_p/2) f_{vR}
\end{equation}
This is then subtracted from the focal length $f_l$ as calculated above, and 
assigned to a new variable $\mathrm{dFocalPastPhm} = f_l - d_{hs}$.  Since $f_l$ and 
$d_{hs}$ are derived from the phantom dimension and the input focal length and field of view ratios, one can write, 
\begin{equation}
\mathrm{dFocalPastPhm} = f_l -d_{hs} 
	= \sqrt{2}(l_p/2) f_{lR} - (l_p/2) f_{vR} = l_p(\sqrt{2}f_{lR} - f_{vR})
\end{equation}
If this quantity is less than or equal to zero, then at least for some
projections  the source is inside the phantom.  Perhaps a figure will help at
this point. Consider first the case where $f_{vR} = f_{lR} =1 $, figure 3. The
square in the figure bounds the phantom and has sides $l_p$.  For this case
then, 
\[ 
f_l=\sqrt{2}l_p/2 = l_p/\sqrt{2},
\]
\[
f_v = \sqrt{2}l_p, 
\]
and 
\[
d_{hs} = {l_p}/{2}.
\]
Then 
\[
\mathrm{dFocalPastPhm} = ({l_p}/{2}) (\sqrt{2}-1)
\]
\begin{figure}
\includegraphics[height=0.5\textheight]{ctsimfig3.eps}
\caption{Equilinear and equiangluar geometry when focal length ratio =
field of view ratio = 1.}
\end{figure}
The angle $\alpha$ is now defined as shown in figure 3, and the detector
length is adjusted to subtend the angle $2\alpha$ as shown.  Note that the
size of the detector array may have changed and the field of view is not
used.  
For a circular array of detectors, the detectors are spaced around a
circle covering an angular distance of $2\alpha$.  The dotted circle in
figure 3 indicates the positions of the detectors in this case. Note that 
detectors at the ends of the range would not be illuminated by the source.

Now, consider increasing the focal length ratio to two leaving the
field of view ratio as 1,  as in  Figure 4.  Now the detectors array is
denser, and the real field of view is closer to that specified, but note
again that the field of view is not used. Instead, the focal length is
used to give a distance from the centre of the phantom to the source, and
the detector array is adjusted to give an angular coverage to include the
whole phantom.  
\begin{figure}
\includegraphics[width=\textwidth]{ctsimfig4.eps}
\caption{Equilinear and equiangluar geometry when focal length ratio = 2
and the field of view ratio = 1.}
\end{figure}
Now consider a focal length ratio of 4 (figure 5). As expected, the angle
$\alpha$ is smaller still.  The dotted square is the bounding square of
the phantom rotated by 45 degrees, corresponding to the geometry of a
projection taken at that angle.  Note that the fan beam now clips the top
and bottom corners of the bounding square.  This illustrates that one may
still be clipping the phantom, despite CTSim's best efforts.  You have
been warned.  
\begin{figure}
\includegraphics[width=\textwidth]{ctsimfig5.eps}
\caption{Equilinear and equiangluar geometry when focal length ratio = 4.}

\end{figure}


\section{Reconstruction}\label{conceptreconstruction}\index{Concepts,Reconstruction}%
\subsection{Overview}
\subsection{Filtered Backprojection}
\subsection{Direct Inverse Fourier}
This method is not currently implemented in \ctsim, however it is planned for a
future release. This method does not give as accurate result as filtered
backprojection mostly due to interpolation occuring in the frequency domain rather
than the spatial domain. The technique is comprised of two sequential steps:
filtering projections and then backprojecting the filtered projections. Though
these two steps are sequential, each view position can be processed individually.
This parallelism is exploited in the MPI versions of \ctsim where the data from
all the views are spread about amongst all of the processors. This has been testing
in a 16-CPU cluster with good results.

\subsubsection{Filter projections}
The projections for a single view have their frequency data multipled by
a filter of absolute(w). \ctsim permits four different ways to accomplish this
filtering. Two of the methods use convolution of the projection data with the
inverse fourier transform of absolute(x). The other two methods perform an fourier
transform of the projection data and multiply that by the absolute(x) filter and
then perform an inverse fourier transform.

\subsubsection{Backprojection of filtered projections}