\documentstyle[agupp,aguplus]{article} % PREPRINT, 
%\documentstyle[jgrga,aguplus]{article}          % CAMERA-READY 
%\documentstyle[agums,12pt]{article}         % MANUSCRIPT
\def\beps{\mbox{\boldmath{$\varepsilon$}}}
\def\bSigma{\mathbf{\Sigma}}
\columnsep=0.4cm
\input epsf 
\setcounter{secnumdepth}{4}
%\afour

\authoraddr{H.~U. Frey, S. Frey, D. Larson, Space Sciences Laboratory, 
University of California, Berkeley, CA 94720-7450, email: 
hfrey@ssl.berkeley.edu.}

\authoraddr{T. Nygr\'en, Department of Physical Sciences, University of Oulu, 
FIN-90570 Oulu, Finland}

\authoraddr{J. Semeter, Max-Planck-Institute for Extraterrestrial Physics, 
D-85740 Garching, Germany}

\begin{document}
\title{Tomographic methods for magnetospheric applications}
\author{H.~U. Frey, S. Frey, D. Larson}
\affil{Space Sciences Laboratory, University of California, Berkeley, CA 
94720-7450, U.S.A.}
\author{T. Nygr\'en}
\affil{Department of Physical Sciences, University of Oulu, 
FIN-90570 Oulu, Finland}
\author{J. Semeter}
\affil{Max-Planck-Institute for Extraterrestrial Physics, D-85740 Garching, 
Germany}

\begin{abstract} The mathematics of tomography was developed for medical
applications when the patient may be examined under controlled circumstances.
Geophysical applications can generally not meet this situation and result in
erratic or impossible reconstructions if minimum requirements are not fulfilled.

This paper summarizes several tomographic reconstruction methods and shows the
limitations and minimum requirements for magnetospheric applications.  The
reconstruction of a simple phantom is used as a point of comparison for the
various methods.  Under experimental conditions with a very limited angular
range and a small number of observers, algebraic and statistical reconstruction
methods generally yield better results than transform techniques.
\end{abstract}

\section{Introduction}\label{sec:intro} 
The aim of tomography is to reconstruct an unknown n-dimensional source
distribution from a set of (n-1)-dimensional remote measurements.  In the
mathematical sense it is an inverse problem [{\it Natterer,} 1986].  The best
known application is in medical diagnosis whereby images of the interior of the
body are reconstructed from a set of X-ray projections.

Tomography in geophysics has long been used for the reconstruction of geological
structures [{\it Menke,} 1984 and references therein] but just recently it has
found applications in ionospheric research.  Ionospheric tomography reconstructs
the electron density distribution in an altitude-latitude plane from satellite
measurements of the total ionospheric electron content [{\it Raymund et al.,}
1995].  These measurements assume spatial homogeneity in longitude and temporal
stationarity.  Thus, temporal fluctuations can not be studied with this method,
and satellite passes need not necessarily coincide with interesting events.
Nevertheless, this method achieved highly desired gross-scale results in large
polar regions and radars provided fundamental independent data used for the
verification of ionospheric reconstructions [{\it Walker et al.,} 1996; {\it
Nygr\'en et al.,} 1996].

Scanning photometers and cameras can be used in a similar way to obtain
two-dimensional [{\it Aso et al.,} 1990; {\it Vallance Jones et al.,} 1991; {\it
Doe et al.,} 1997] or three-dimensional auroral volume emission reconstructions
[{\it Frey et al.,} 1996; 1998].  As in general geophysical tomography, auroral
tomography with ground-based optical instruments is restricted to a small number
of cameras within a limited angular range.  Proper reconstruction methods have 
to account for these limitations.

The aim of this paper is to summarize different reconstruction techniques which
have been used in ionospheric research and compare the results for
magnetospheric reconstructions from radio sounding measurements.  This
comparison is made with a phantom plane of electron density distributions.  Even
a carefully-designed multi-satellite mission cannot guarantee tomographic
measurements from a complete set of angles around the object.  An important
point of comparison, therefore, is the sensitivity of the various techniques to
incomplete data.

\section{Tomographic methods}\label{sec:method}
\subsection{Remote measurements} \label{sec:remote}

Many investigations determine the integral of some property along
lines from the observer position at ${\bf r}=0$ through the object as the
projection $p_i(\theta, \phi, t)$

\begin{equation} \label{eq:allg_tomo1} 
p_i(\theta,\phi,t)=C \int\limits_{0}^{\infty} 
f({\bf r}(l)) \exp \left\{-\int\limits_{s=0}^{s=l} \mu(s) ds \right\} {\rm d} l 
\end{equation}

\noindent
with the source distribution $f({\bf r})$ and a general absorption
term $\mu(s)$.  The unknown quantity is either the source distribution, or the
distribution of the absorbing matter, or both.  Computed tomography is then used
to determine the distribution of the unknown quantity from a set of line
integral measurements.

In the computerized discrete case, the plane is divided into small pixel, with
the value in each pixel $(x,y)$ representing the mean of the source distribution
$f(x,y)$.  For the purpose of better visualization the following presentation
shows the reconstruction of a 2D-distribution, but in general it can also be
used for 3D problems.

\begin{figure} \vspace*{7.0cm} 
\special{psfile=/disks/sprite/disk1/hfrey/paper/tomo_agu/parallel.eps angle=0 
hscale=40 vscale=40}
\caption{\label{fig:parallel} Parallel beam geometry.}  \end{figure}
\begin{figure} \vspace*{7.0cm} 
\special{psfile=/disks/sprite/disk1/hfrey/paper/tomo_agu/fan.eps angle=0 
hscale=40 vscale=40}
\caption{\label{fig:fan} Fan beam geometry.}  \end{figure}

The parameter $C$ contains measurement parameters.  In this investigation we
want to approximate the absorption as constant and (\ref{eq:allg_tomo1})
simplifies to

\begin{equation} \label{eq:allg_tomo2} p_i(\theta ,t)
\approx C'\sum\limits_{x,y} h^{i \theta t}_{xy} f(x,y)\end{equation}

\noindent
with $h^{i \theta t}_{xy}$ describing the contribution of the pixel (x,y) to 
the region subtended by the ray for the observer $i$.

The measurement geometries are based on parallel (Figure \ref{fig:parallel}) or
diverging fan (Figure \ref{fig:fan}) beams.  Diffracting media need full angular
coverage around the object.  However, for non-diffracting media parallel
projections which are 180$^{\circ}$ apart are just redundant mirror images of
each other [{\it Kak and Slaney,} 1988], and it is sufficient to measure
projections for angles from $180^\circ$. For measurements with a fan-beam angle 
of $\gamma$ projections in an angular region of 

\begin{equation} \label{eq:angle2}
\theta_0 \le \theta < \theta_0 + 180^\circ + 2 \gamma
\end{equation} 

\noindent
are necessary.  Even this angular range can not be provided in most
geophysical applications which leads to the so-called missing cone problem [{\it
Hayner and Jenkins,} 1984].  One further characteristics of geophysical
applications is that the objects are generally smooth with features of similar
contrast and there are no sharp borders like bones in medical tomography.

The goal of tomography is to solve the set of linear equations
(\ref{eq:allg_tomo2}).  Various methods have been developed to find an estimate
$\hat{f}$ of $f$ [for reviews, see e.g., {\it Budinger and Gullberg,} 1974; {\it 
Natterer,} 1986; {\it Verhoeven,} 1993], and they can be classified into three 
main categories:
\begin{itemize}
\item{Fourier transform methods}
\item{Algebraic methods }
\item{Statistical methods}
\end{itemize}

\subsection{Quality parameters for the techniques} \label{sec:parameters}

Any tomographic technique has to be tested with models in order to establish a
quantitative criterion for the quality of reconstruction.  The correlation
coefficient is only a qualitative measure how well the model and the
reconstruction vary jointly, but is does not represent the quantitative
agreement.

Adapted from medical investigations [{\it Shepp and Logan,} 1974], a phantom is
used which might resemble density structures in a magnetospheric plane with
regions of lower and higher electron density than the surrounding area (Figure 
\ref{fig:filtback}).

\begin{figure} \vspace*{5.0cm} 
\special{psfile=/disks/sprite/disk1/hfrey/idl/tests/paper_fbp/filtback.eps 
angle=0 hscale=100 vscale=100}
\caption{\label{fig:filtback} Improvement of the filtered backprojection 
reconstruction with increasing number of observations. The upper left panel 
contains the phantom which is used with all methods and shows the circle of 
observer locations.  The relative error $\Lambda$ in \% is indicated.}  
\end{figure}

In any practical case only the observations (line integrals) are known
quantities and any theoretical treatment of the observation process applied to
the reconstructed data has to provide good agreement with the measured data
after the influence of noise has been taken into account.  As the most
important goal of tomography is the agreement between the reconstructed
$\hat{f} (x,y)$ and the original distribution $f(x,y)$, a quantitative reconstruction
parameter $\Lambda_k$ is defined: 

\begin{equation} \label{eq:lambda} 
\sum_{x} \sum_{y} \sqrt{[f(x,y)
- \hat{f}(x,y)]^2} =  \Lambda_{k} \cdot \sum_{x} \sum_{y} f(x,y).
\end{equation}

\subsection{Fourier transform methods} \label{sec:fourier}

Fourier transform (FT) techniques use the Fourier Slice Theorem which states
that the FT of a projection is equal to the slice through the FT of the original
object [{\it Kak and Slaney,} 1987].  Thus, given FTs ${\cal F}(p_i)$ of
projections $p_i$ from the complete angular range around the object will provide
the original distribution by simple inversion.  However, in reality, only a
finite number of observations is available.  Hence, interpolation between the
FTs of the projections in the frequency domain is necessary to obtain the
complete FT of the object.  Since the FT of the individual projection can be
regarded as a filtered FT of the object, this interpolation is mostly done after
a weighting of the projection FT in the frequency domain by filtering with a
filter function $\tilde c$.  This leads to the Filtered Back-Projection (FBP)
method which is mostly used in medicine [{\it Kak and Slaney,} 1988]

\begin{equation} \label{eq:fbak}
\hat{f} = \sum\limits_{i} b^{xy}_{i \theta t}~{\cal F}^{-1} \left [ {\tilde c} {\cal 
F}(p_i) \right ]
\end{equation}

The weighting factor $b^{xy}_{i \theta t}$ is the inverse of $h^{i \theta
t}_{xy}$ in (\ref{eq:allg_tomo2}) and establishes the geometrical 
correspondence between the object and the projection.

Figure \ref{fig:filtback} shows how the quality of the FBP reconstruction
improves with an increasing number of equally separated observations.  In this
investigation FBP is the only method used with parallel beams.

\subsection{Algebraic methods} \label{sec:algebraic}

Direct inversion methods have only limited use in geophysical inverse problems.  
The matrices involved are often ill-conditioned and so direct matrix inversion 
is not possible. Inversion by singular value decomposition handles many problems 
associated with ill-conditioning but is far too computationally intensive for 
the large dimensions encountered in many real geophysical data sets. 
 
Algebraic reconstruction techniques (ART) solve the problem in an iterative
scheme.  Several types of ART differ in the way the corrective term is
determined while the iteration is always performed on a ray by ray basis until
the calculated projection $p^k_{i,l(\theta, t)}$ of the iteration $k$ is as
close as possible to the measured data $p^e_{i,l(\theta, t)}$ [{\it Gilbert,}
1972].

Mainly, we distinguish three types:\\
additive ART (ART) with the non-negative constraint
\begin{equation} \label{eq:back_art1} 
\hat{f}^{k+1}(x,y) = \hat{f}^{k}(x,y) +  b^{xy}_{i 
\theta t}~(p^e_{i,l(\theta, t)} - p^k_{i,l(\theta, t)}) ~~ \ge 0
\end{equation}
multiplicative ART (MART)
\begin{equation} \label{eq:back_art2} 
\hat{f}^{k+1}(x,y) = \hat{f}^{k}(x,y) \cdot  b^{xy}_{i \theta t}~
\frac{p^e_{i,l(\theta, t)}}{p^k_{i,l(\theta, t)}} 
\end{equation}
simultaneous iterative ART (SIRT) 
\begin{equation} \label{eq:back_art3} 
\hat{f}^{k+1}(x,y) = \hat{f}^{k}(x,y) \cdot \frac{1}{I} \sum\limits_{i} b^{xy}_{i \theta 
t}~
\frac{p^e_{i,l(\theta, t)}}{p^k_{i,l(\theta, t)}} \end{equation}

All methods use back-projection (BP) [{\it Budinger and Gullberg,} 1974] which
projects and sums each projection to the object space filling each line with the
data of the projection [{\it Peyrin,} 1985]
\begin{equation} \label{eq:back_eq1} \hat{f}(x,y) = \frac{T}{T'} \cdot
\sum\limits_{i,\theta ,t} b^{xy}_{i \theta t}~p_{i}(\theta ,t) \end{equation}
with establishing the quantitative correspondence with
the density $T$ of the element and the density $T'$ of the complete array.
This method is very fast but the many simplifications generally return too high
cell contents.

\begin{figure} \vspace*{11.1cm} 
\special{psfile=/disks/sprite/disk1/hfrey/idl/tests/paper_mart/algebraic.eps 
angle=0 hscale=98 vscale=98}
\caption{\label{fig:algebraic} Reconstructions of the phantom by algebraic 
methods from 8, 32, and 256 views, respectively with fields of view of 
90$^\circ$ and 120$^\circ$.}  
\end{figure}

The iterations generally start with a homogeneous first estimate.  Its selection
is of main importance for the applicability of a method and for the number of
iterations.  In a modification of the MART {\it Frey et al.,} [1996] reconstruct
the first estimate from the BP of projections
\begin{equation} \label{eq:minimum}
\hat{f}(x,y) = \min\limits_{i} \left\{ b^{xy}_{i \theta t}~p_{i}(\theta ,t) \right\}
\end{equation}
where each element is assigned the minimum of the back-projected data.
This modification provides the advantage that a cell content once set to 0 
cannot increase. 

Figure \ref{fig:algebraic} shows reconstructions with ART, SIRT, and modified
MART with 90$^\circ$ and 120$^\circ$ field of view.  For the given geometry the
object is bigger than the region seen with the 90$^\circ$ field of view.  This
makes the inversion less accurate than in the case of 120$^\circ$ field of view,
when the outer border of the object is recognized.  

MART has successfully been used for the 3D reconstruction of the volume emission
in an auroral arc [{\it Frey et al.,}  1996a].  The good correspondence between
the reconstructed volume emission and the electron density measured by the
EISCAT radar was later used to cross-calibrate the volume emission and to
calculate the ionospheric conductivity [{\it Frey et al.,}  1998].  This
large-scale conductivity distribution together with the radar measured electric
field can then be used to determine ionospheric properties in and outside an
auroral arc.

If the original image is known to be reasonably smooth then a least squares
(LS) method can be used.  This method assumes the original image can be
represented by a smooth function of two independent variables and some large
number of parameters.  Typically these parameters could be the value of the
function at fixed locations and the function would perform a 2-D
interpolation between these points.  Here we use the function:
\begin{eqnarray} \label{eq:least_sq}
\hat{f} (x,y) = C_0 + C_1 x + C_2 y + \sum\limits_{n=3} C_n D_n ln D_n \\
D_n = \sqrt{(x-X_n)^2 + (y-Y_n)^2} \nonumber 
\end{eqnarray}
where $(X_n,Y_n)$ are the coordinates of the fixed points and $C_n$ are
adjustable parameters.  These points need not to be on a regular grid.  As the
line integrals $p_i$ then become linear functions of the parameters $C_n$,
standard linear least squares analysis can be used to determine the values of
the $C_n$ which best fit the measured values of $p_i$.  Figure
\ref{fig:least_sq} shows results from 16 and 32 observations, and the resolution
in this case can not be better than the distance between the grid points.

Alternatively, the density n(x,y) can be equated to the exponential of
$\hat{f}(x,y)$, i.e.  $n=exp(\hat{f})$.  Although non-linear least squares
analysis is then required, this has the benefit of forcing the density to be
positive definite.

\begin{figure} \vspace*{2.9cm} 
\special{psfile=/disks/sprite/disk1/hfrey/idl/dreidim/davin/least_sq.eps 
angle=0 hscale=47 vscale=47}
\caption{\label{fig:least_sq} Reconstructions from a least squares method. 120 
lines connect 16 observation points and each line represents one line integral 
(left).  The reconstructed image shows the location of fixed points $(X_n ,Y_n)$ 
(center). The reconstruction from 32 observations is shown to the right.}  
\end{figure}

The natural pixel (NP) object representation was developed for tomography from
incomplete data and irregularly viewing geometries [{\it Buonocore et al.,}
1981], and has recently been applied to the problem of ionospheric radio
tomography [{\it Semeter and Kamalabadi,} 1998].  In NP tomography, the object
is represented by a basis expansion.  Consider a set of tomographic measurements
through a region of interest $\Omega$ containing $f(x,y)$.  A natural pixel is a
characteristic function with the following definition:
\begin{equation} \label{eq:char}
\phi_i(x,y)=\left\{ \begin{array}{ll}
          1\hspace{.2cm}\hbox{\footnotesize{(x,y) lies in  
                         field-of-view of measurement $i$}}\\
          0\hspace{.2cm}\hbox{\small{otherwise}}
                    \end{array}
            \right.
\end{equation}

The measurement $p_i$ is specified by
\begin{equation} \label{eq:moment}
p_i=\int_\Omega f(x,y)\phi_i(x,y)dxdy
\end{equation}

Now assume that $\hat{f}(x,y)$ is adequately represented by a linear combination in 
this basis, i.e.,
\begin{equation} \label{eq:fhat}
\hat{f}(x,y)=\sum_{j} q_j\cdot\phi_j(x,y)
\end{equation}
where the $q_j$ are the coefficients of the NP expansion. Substituting 
(\ref{eq:fhat}) into (\ref{eq:moment}) we have
\begin{equation} \label{eq:phiphi}
p_i=\sum_{j} q_j\underbrace{\int_\Omega\phi_j(x,y)\phi_i(x,y)
dxdy}_{C_{ij}}
\end{equation}

\begin{figure} \vspace*{3.0cm} 
\special{psfile=/disks/sprite/disk1/hfrey/idl/tests/paper_josh/josh.eps 
angle=0 hscale=100 vscale=100}
\caption{\label{fig:natpix} Reconstructions of the phantom from natural pixel 
presentation. Each equally-spaced projection uses 32 equally-spaced rays.}  
\end{figure}

The term $c_{ij}$ is the inner product of natural pixel $i$ with natural pixel
$j$.  These form the elements of the NP projection matrix $C$ whose elements are
the areas of polygons defined by overlapping ray boundaries.  Tomographic
reconstruction involves solving the matrix equation $p=Cq$ [{\it Buonocore et
al.,} 1981].  The coefficients $q$ are then backprojected into a displayable
grid.  Fig.  \ref{fig:natpix} shows an example of NP-reconstruction.  Note that
useful features are reconstructed with only 4 projections, but that the results
are generally worse than MART as the number of views increases.

In principle, NP offers many advantages over other techniques for image
reconstruction from incomplete data.  Useful results are achieved at low
computational overhead and are devoid of many of the incomplete data artifacts
present in other techniques.  The NP-representation leads to a framework by
which scale-dependent features in the data can be related to scale-dependent
features in $f(x,y)$ [{\it Bhatia et al.,} 1997].

However, there are also significant drawbacks to the NP representation.  The
matrix $C$ is full and individual projections cannot be processed separately.
Thus tremendous storage is required for large problems.  There is also no direct
way to impose solution positivity, although {\it Semeter and Kamalabadi} [1998]
have recently proposed a technique for refining NP reconstructions under a
positivity constraint.

\subsection{Statistical methods} \label{sec:statistic}

Statistical methods consider the measurements as a stochastic process and 
look for the most probable solution. They mostly use Maximum Likelihood (ML) 
methods which reconstruct the original distribution in an iteration process as 
[{\it Shepp and Vardi,} 1982]
\begin{equation} \label{eq:maxlike}
\hat{f}^{k+1} (x,y) = \hat{f}^k (x,y) \sum\limits_{i} \frac{p_i(\theta, t) c_{xi}}
  {\sum\limits_{n=1}^{N} c_{ni} \hat{f}^k (x,y)}
\end{equation}
with the probability $c_{xi}$ of detection of an emission event in pixel x in 
the projection i.

The ML has successfully been used in positron emission tomography and the
main advantage is the solid theoretical basis, the good convergence, and the
robustness against statistical noise [{\it Shepp and Vardi,} 1982].

\begin{figure} \vspace*{3.2cm} 
\special{psfile=/disks/sprite/disk1/hfrey/idl/dreidim/tuomo/tuomos/tuomo.eps 
angle=0 hscale=100 vscale=100}
\caption{\label{fig:stochastic} Results of stochastic inversion.}  \end{figure}

Stochastic inversion (SI) has been used in reconstructing the F region electron
density from difference Doppler measurements of total electron content, obtained
from navigation satellites and a chain of ground-based receivers [{\it
Nygr\'en et al.,} 1996].  The results have been compared to simultaneous
electron density measurements carried out by the EISCAT incoherent scatter
radar.

In this method all measurements are collected into a single column vector {\bf
p}, their errors into a column vector $\beps_p$, and the unknowns into a column
vector {\bf f} with\ ${\bf p} = {\bf A_p} \cdot {\bf f} + \beps_{\bf p}$, where
${\bf A_p}$ is a matrix of linear coefficients.  Similarly, writing an equation
$0 = f_n - f_m + \varepsilon_r^{(nm)}$ for the unknowns $f_n$ and $f_m$ at each
pair $(n,m)$ of neighbouring grid points gives a matrix equation ${\bf 0} = {\bf
A_r} \cdot {\bf f} + \beps_{\bf r}$.  Here the components of $\beps_{\bf r}$ are
the errors made when the differences $f_n - f_m$ are put to zero.  Both {\bf p},
$\beps_p$, {\bf f}, and $\beps_r$ are treated as random variables.  Assuming 
that the distributions of $\beps_p$ and $\beps_r$ are Gaussian with
covariance matrices $\bSigma_{\bf p} = \langle \beps_{\bf p} \cdot \beps_{\bf
p}^T \rangle$ and $\bSigma_{\bf r} = \langle \beps_{\bf r} \cdot \beps_{\bf r}^T
\rangle$, respectively, and the mean value $\langle \beps_{\bf r} \rangle = {\bf
0}$, the conditional density of the unknowns is [{\it Nygr\'en et al.,} 1997]
\begin{displaymath}
D({\bf f}|{\bf p},{\bf 0}) \propto D_{pr}({\bf f}) \cdot
\exp \left[ -{1 \over 2}({\bf A_r} \cdot {\bf f})^T
\cdot \bSigma_{\bf r}^{-1} \cdot 
({\bf A_r} \cdot {\bf f}) \right] 			
\end{displaymath}
\begin{equation} \label{eq:si_sol}
\cdot \exp \left[ -{1 \over 2} ({\bf A_p} \cdot {\bf f} - {\bf p})^T
\cdot \bSigma_{\bf p}^{-1} \cdot ( {\bf A_p} \cdot {\bf f} - {\bf p}) \right] 
\end{equation}

Here $T$ indicates the transpose.  The inversion solution is obtained as a
maximum point of the conditional density $D({\bf f}|{\bf p},{\bf 0})$ or as the
mean value of {\bf f}.  If no a priori information is used, these solutions are
identical.  Then the inversion problem is linear, and the maximum likelihood
point is simply given by 
\begin{equation} \label{eq:ml_point}
{\bf \hat{f}} = ({\bf A_p}^T \cdot \bSigma_{\bf p}^{-1} \cdot {\bf A_p}
 + {\bf A_r}^T \cdot \bSigma_{\bf r}^{-1} \cdot {\bf A_r})^{-1}
 {\bf A_p}^T \cdot \bSigma_{\bf p}^{-1} \cdot {\bf p}
\end{equation}

Here ${\bf A_r}^T \cdot \bSigma_{\bf r}^{-1} \cdot {\bf A_r}$ is a non-diagonal
regularization matrix which stabilizes the solution.  A smoothly behaving
solution is obtained by means of a proper choice of the covariance matrix
$\bSigma_{\bf r}$.

If the unknown function cannot have negative values, the non-negativity
constraint can be implemented by means of an a priori density $D_{pr}({\bf f})$,
which makes the inversion problem nonlinear.  Then, instead of
(\ref{eq:ml_point}), it is possible to search for a set of random vectors ${\bf
\hat{f}}$ obeying the posteriori distribution in (\ref{eq:si_sol}).  This can be
done by means of Gibbs sampling.  When the set is large enough, the average of
the sample vectors gives a good estimate to the statistical average of ${\bf
\hat{f}}$.  In this case the average is given as the inversion result instead of
the maximum likelihood point (Figure \ref{fig:stochastic}).

\section{Results of the methods}\label{sec:results}

For tests of the reconstruction methods observations of the phantom from 4 to
256 directions were simulated.  The accuracy of tomography may be improved by
the introduction of a-priori knowledge about the object like smoothness, the
extent, or the non-negativity.  For algebraic and statistical techniques it is
generally easy to incorporate this knowledge to limit the number of possible
solutions, or speed-up the iteration [{\it Trussel,} 1984].  FT techniques can
generally not include this information.

If an object is to be uniformly resolved to a resolution d, the number of
discrete views must be at least 

\begin{equation} \label{eq:resolution}
n \approx \pi D / d
\end{equation}

\noindent
with the dimension of the object D [{\it Budinger and Gullberg,} 1974].  With 
our phantom at least 25 equally separated views are necessary to resolve the 
smaller structures.

\begin{figure} \vspace*{5.0cm} 
\special
{psfile=/disks/sprite/disk1/hfrey/idl/tests/paper_mart/fehler_methods.eps 
angle=0 hscale=50 vscale=40}
\caption{\label{fig:fehler_ma} Mean error of methods for an increasing number of 
observers around the phantom.} 
\end{figure}

For this example the MART and SI yield the best results (Figure
\ref{fig:fehler_ma}).  While the SI reached errors $<$ 10\% with as low as 8
observations, the MART got better results than SI with more than 64
observations.

The difference between the methods is that all used 128 rays from each observer
except the LS with only (n-1) rays and NP with 32 rays.  This limited the total
number of useful rays for LS to n(n-1)/2 and from the corresponding small number
of grid points, the results of this method for 16-32 observers are quite good.

The problem with all methods is that there is no confirmed mathematical or
physical reason to set particular pixel elements to zero, except for the
modified MART where this is done at least for all the pixels along a ray from
one of the observers which showed zero content in the image.  Because the total
intensity within the images has to correspond to the total source distribution,
some pixels at the outer surface of the object contain intensity which is then
missing in the center pixels.  This is also most probably the reason why radio
tomography tends to overestimate the electron density on the topside of the
reconstruction [{\it Walker et al.,} 1996] and to produce thicker layers [{\it
Nygr\'en et al.,} 1996].

Comparing the results in Figure \ref{fig:fehler_ma} it has to be kept in mind,
that all methods used iterations (up to 64) in order to obtain the small
resulting errors, except the FBP and SI, which only used one reconstruction.
The FBP is a good method if many observations from the full angular range are
available and if no iterations are used.  Since this is the exception rather
than the rule in geophysical applications, use of this method is limited.

The SIRT is a good method if many observations are available, because in this
case the mean values in the update scheme produce sufficient effects.  While ART
and SIRT encounter problems in determining the areas with lower density than the
surrounding, the modified MART can at least determine the circular area of our
phantom.  The iterations provide an opportunity to recover sufficient
reconstructions already from 16 observations ($ 23^\circ$ separation).  However,
as predicted by Equation (\ref{eq:resolution}) a minimum number of 32
observations ($\approx 12^\circ$ separation) was necessary for most methods to
achieve a satisfactory quantitative reconstruction of the phantom (Figure
\ref{fig:fehler_ma}).

\begin{figure} \vspace*{5.0cm} 
\special{psfile=/disks/sprite/disk1/hfrey/idl/tests/paper_mart/miss_cone.eps 
angle=0 hscale=50 vscale=40}
\caption{\label{fig:filt_miss} Quality parameter $\Lambda$ for
a reconstruction with a missing cone in the observation directions starting with 
32 equally spaced observers around the phantom. FBP has been adjusted to the 
parallel beam geometry.} 
\end{figure}

The algebraic and statistical methods are more robust against missing cones than
the FBP.  This is shown in Figure \ref{fig:filt_miss}, where the mean error is
given for 32 equally separated observations and a missing cone region from
$\theta = 0^\circ$ to $\theta = \theta_1$.  The agreement between the original
phantom and the reconstruction is best for the SI.  It provides the best
reconstruction of the overall shape of the object.  The agreement is still
sufficient for the modified MART as long as the missing cone does not reach more
than 200$^\circ$.

Noise always deteriorates the quality of real images.  As a quantitative test
noise was added to the projections with increasing noise and decreasing
signal-to-noise ratio S/N [{\it Verhoeven,} 1993]
%\begin{equation} \label{eq:sdn}
%S/N = 10 \log \left(\frac{\sum\limits_{i} p_i^2}{\sum\limits_{i} \Delta 
%p_i^2} \right) \end{equation}
\begin{equation} \label{eq:sdn}
S/N = 10 \log \left(\sum\limits_{i} p_i^2 / \sum\limits_{i} \Delta 
p_i^2 \right) \end{equation}
where $\Delta p_i$ is the ray sum error caused by the noise (Figure
\ref{fig:noise}).  The SI is very robust against noise and can handle S/N values
down to 10 dB.  Without noise reduction the results of the MART get already much
worse at signal to noise ratios of less than 50.  The non-negative constraint
improves the situation, but the median filter best improves to acceptable signal
to noise ratios greater than 20.  In this way some knowledge about the nature of
noise sources may help to improve MART results by proper filtering.

\begin{figure} \vspace*{5.0cm} 
\special{psfile=/disks/sprite/disk1/hfrey/idl/tests/paper_mart/noise.eps 
angle=0 hscale=50 vscale=40}
\caption{\label{fig:noise} Effect of noise on the reconstruction parameter.}  
\end{figure}

\section{Conclusions}\label{sec:conclus}

The goal of this paper was to compare different tomographic reconstruction
methods to find those which might work best for a magnetospheric
radio-tomography mission with the limitations of a small number and probably a
limited angular range of observations.  The phantom which was used for this
investigation contained not only regions with higher density than the
surrounding, but also regions with lower density.

Overall the best results were obtained with the SI, not only for a small number 
of observers, but also in view of missing cones and the robustness against 
noise. The method is computationally quite extensive, however still very fast 
with about 10 minutes calculation time because it does not use iterations.

The second best methods are the modified MART, NP and LS.  MART provides very
good reconstructions with more than 32 observers, tolerates large missing cones,
and can adjust noise after preprocessing.  The LS and NP yield good results for
a small number of observers.  

The FBP, ML, and SIRT methods have bigger problems with the reconstruction
quality and might not be very well suited for magnetospheric applications.

This investigation was performed for magnetospheric satellite
radio tomography. For the purpose of ionospheric satellite radio tomography, the
value of the comparisons made in this paper is limited by the fact that, in that
case there is always a large sector around the horizontal direction with no
observations.  These missing observations would contain information on the
vertical ionospheric profile, and the absence of this information poses the main
problem in satellite radio tomography.
\acknowledgement{We are grateful to M. Markkanen for carrying out the SI 
calculations. Support from the Academy of Finland is acknowledged by T.N.}
\begin{references}
\reference
Aso, T., T.  Hashimoto, M.  Abe, T.  Ono, and M.  Ejiri, On the analysis of
aurora stereo observations, {\it J.  Geomagn.  Geoelectr., 42,} 579--595, 1990.

\reference
Bhatia, M., W. C. Karl, and A. S. Willsky, Tomographic reconstruction and
estimation based on multiscale natural-pixel bases, {\it IEEE Trans. Image 
Proc., 6,} 463, 1997.

\reference
Budinger, T.~F., and G.~T.  Gullberg, Three-dimensional reconstruction in
nuclear medicine emission imaging, {\it IEEE Trans.  Nucl.  Sci., NS-21},
2--20, 1974.

\reference
Buonocore, M. H., W. R. Brody, and A. Macovski, A natural pixel decomposition
for two dimensional image reconstruction, {\it IEEE Trans. Biomed. Eng., 
BME-28,} pp. 69-78, 1981.

\reference
Doe, R.~A., J.~D. Kelly, J.~L. Semeter, and D.~P. Steele, Tomographic 
reconstruction of 630.0 nm emission structure for a polar cap arc, {\it Geophys. 
Res. Lett., 24,} 1119-1122, 1997.

\reference
Frey, S., H.  U.  Frey, D.  J.  Carr, O.  H.  Bauer, and G.  Haerendel, Auroral
emission profiles extracted from three-dimensional reconstructed arcs, {\it J.
Geophys.  Res., 101,} 21731--21741, 1996a.

\reference
Frey, H.~U., S. Frey, B.~S. Lanchester, and M. Kosch, Optical tomography of the 
aurora and EISCAT, in press in {\it Ann. Geophys.}, 1998.

\reference
Gilbert, P., Iterative methods for the three dimensional reconstruction of an 
object from projections, {\it J. Theor. Biol., 36,} 105--117, 1972.

\reference
Hayner, D.~A. and W.~K. Jenkins, The missing cone problem in computer 
tomography, in {\it Advances in Computer Vision and Image Processing,} JAI 
Press, Tokyo, 1984.

\reference
Kak, A. C. and M. Slaney, {\it Principles of Computerized Tomographic Imaging}, 
IEEE Press, New York, 1988.

\reference
Menke, W., {\it Geophysical data analysis: Discrete inverse theory,} Academic 
Press, Inc., New York, 1984.

\reference
Natterer, F., {\it The Mathematics of Computerized Tomography}, John Wiley, New
York, 1986.

\reference
Nygr\'en, T., M. Markkanen, M. Lehtinen, E. D. Teresh- chenko, B. Z. Khudukon, 
O. V. Evstafiev, and P. Pollari, Comparison of F-region electron density 
observations by satellite radio tomography and incoherent scatter methods, {\it 
Ann. Geophys., 14,} 1422--1428, 1996.

\reference
Nygr\'en, T., M. Markkanen, M. Lehtinen, E. D. Teresh- chenko and B. Z. 
Khudukon, Stochastic inversion in ionospheric radiotomography, {\it Radio Sci., 
32,} 2359-2372, 1997.

\reference
Peyrin, F.~C., The generalized back projection theorem for cone beam
reconstruction, {\it IEEE Trans.  Nucl.  Sci., 32}, 1512--1519, 1985.

\reference
Raymund, T. D. S.~J. Franke, and K.~C. Yeh, Ionospheric tomography: its 
limitations and reconstruction methods, {\it J. Atmos. Terr. Phys., 56,}
637--657, 1994.

\reference
Semeter, J. and F. Kamalabadi, A natural pixel decomposition for tomographic
imaging of the ionosphere, {\it Proc. Int. Conf. Acoustics, Speech and Signal 
Processing (ICASSP)}, Seattle, Washington, 1998. 

\reference
Shepp, L.~A. and B.~F. Logan, The Fourier reconstruction of a head section, {\it 
IEEE Trans. Nucl. Sci., NS-21,} 21--43, 1974.

\reference
Shepp, L.~A. and Y. Vardi, Maximum likelihood reconstruction for emission 
tomography, {\it IEEE Trans. Med. Imaging, 1,} 113--122, 1982.

\reference
Trussel, H.~J., A priori knowledge in algebraic reconstruction methods, in {\it 
Advances in computer vision and image processing}, JAI Press Inc., 1984.

\reference
Vallance Jones, A., et al., Characterization and modeling of an evening auroral 
arc observed from a rocket and a ground-based line of meridian scanners, {\it 
Planet.  Space Sci., 39,} 1677--1705, 1991.

\reference
Verhoeven, D., Limited-data computed tomography algorithms for the physical 
sciences, {\it Appl. Optics, 32,} 3736--3754, 1993.

\reference
Walker, I. K., J. A. T. Heaton, L. Kersley, C. N. Mitchell, S. E. Pryse, and M. 
J. Williams, EISCAT verification in the development of ionospheric tomography, 
{\it Ann. Geophys., 14,} 1413--1421, 1996.

\end{references}
\end{document}

