next up previous contents
Next: Summary Up: Tomography Previous: Error sensitivity   Contents

Stopping criteria

It has long been known that iterative solutions to tomographic problems improve initially but deteriorate after some point ( Gordon and Herman, 1974). This highlights the importance of knowing when it is time to stop.

What we want to do is to test whether the reconstruction gives projections that are statistically compatible with the observed images. If we had a large number of samples $ x_n$ from one random process, it is possible to test whether the random process has a hypothetical probability distribution $ P_h(x)$. The method to do this is the $ \chi^2$-test of fit for which the range of the random variable is divided into $ L$ intervals, $ \Delta_l$ with theoretical probabilities $ p_l$ for a random sample $ x_i$ to fall into the interval $ l$ if the hypothesis $ P_x(x) \equiv P_h(x)$ is true.

With a total of $ n$ samples $ x_1,\dots,x_n$ and $ M_l$ samples in interval $ l$, the $ \chi^2$-test function

$\displaystyle \chi^2 = \sum (M_l - n\,p_l)^2/(n\,p_l)$ (3.15)

has an asymptotic $ \chi^2$-distribution with $ L-1$ degrees of freedom. $ \chi^2$ is a measure of the deviation between the true distribution and the hypothetical one. The hypothesis is accepted if the value calculated from equation (3.15) is less than a critical value corresponding to a chosen significance level ( Bronstein and Semendyayev, 1997).

In the case of tomography, the pixels in the projection images of the reconstructed three-dimensional source distribution are the expected values of the probability distributions from which the measured pixel values are the random samples. Here we only have one sample from each probability distribution; it should be noted that the probability distribution for a pixel value is known and assumed to be determined only from the expected value.

For auroral imaging with ALIS, the pixel intensity has a probability distribution that is $ p(I(u,v)) = N(I(u,v),\sqrt{I(u,v)})$ as derived in section 3.1. With the above relation between expected pixel intensity and the random distribution of the measured pixel intensity, it is possible to use the above $ \chi^2$-test after the modification that the random process for each pixel is divided into intervals with equal probabilities $ p_l$ and that the measured pixel intensity is put into its corresponding interval.

The hypothesis that the measured images $ I_s$ are random samples from a set of random processes with expected values $ \tilde{I}_s$ is accepted if the $ \chi^2$ value calculated from equation (3.15) is less than is required at the corresponding significance level. (The values of $ \tilde{I}_s$ are calculated from the source distribution.) For processes that have continuous probability distributions, it is no problem to divide the range into equal intervals for different expected values and widths but for discrete probability distributions this can be a problem; for Poissonian processes a method to overcome this problem is described by Veklerov and Llacer (1987).

An example that can illustrate the working procedure outlined above is to look at the outcome of the roll of four non-standard dice. What separates these dice from ordinary dice is that they have unknown numbers of dots; that is, they might have dots from 1 to 6 or from 3 to 8 or any other sequence of dots $ n,n+1,\dots,n+5$. If we make a measurement of the dot numbers and then get two models which estimate the expected value of each die, as presented in table 3.4, the algorithm for this stopping criteria is as follows.

Table 3.1: The number range $ x---x+5$ denotes that the die have faces with $ x,x+1,x+2,x+3,x+4,x+5$ eyes.
		    \centering \begin{tabularx}{100mm}{XXXXX}%{CRRRR}
		    Die ...
		    ... & 1{---}6 & 2{---}7 & 1{---}6 & 3{---}8 \\
		    \end{tabularx} \end{table}

First we divide the probability distribution into a number of intervals, for this case two intervals $ x < E(X)$ with $ p_1 = 0.5$ and $ x > E(X)$ with $ p_2 = 0.5$. For model 1 all events fall into interval 1 and the $ \chi^2$ value according to equation (3.15) is:

$\displaystyle \chi^2_1 = (4-4\cdot0.5)^2/(4\cdot0.5) + (0-4\cdot0.5)^2/(4\cdot0.5) = 4$ (3.16)

and for model 2 die $ 3$ falls into interval $ 1$ and dice $ 1, 2$   and $ 4$ fall into interval $ 2$ giving a $ \chi^2$ value, according to equation (3.15), of:

$\displaystyle \chi^2_1 = (1-4\cdot0.5)^2/(4\cdot0.5) + (3-4\cdot0.5)^2/(4\cdot0.5) = 1$ (3.17)

For a significance level of $ 0.05$ and one degree of freedom, the critical value of $ \chi^2$ is 3.84, i.e. if $ \chi^2$ is larger than 3.84 we must reject the hypothesis that the model is the cause of the observations. Model 1, with $ \chi^2 = 4$, is not likely given the observations and has to be rejected, and for model 2, with $ \chi^2 = 1$, we cannot reject the hypothesis and should be satisfied and stop the iteration.

If we have come this far the reconstruction gives projections that fit the measured images in a statistical sense. However, there are not yet any guarantees that the spatial distribution of the errors is well behaved, i.e., the large positive errors might be grouped together in one region of the images and small errors might be grouped in another region. One further condition for optimal reconstructions is that the spatial distribution of the errors is also random. This can be tested with the same method as above: the image location of errors in an intensity interval $ l$ is calculated and then the images are divided into 10 by 10 regions and the number of errors $ E_l$ from the interval $ l$ is calculated region by region. If the $ \chi^2$-test function

$\displaystyle \chi^2 = \sum n_{reg}\,E_l^2/n_e - n_e$ (3.18)

where $ n_e$ is the total number of errors in the intensity error interval $ l$ and $ n_{reg}$ is the total number of regions, is also less than the corresponding significance value, there is no reason to continue with the iteration.

By the way, about the cartoon animation on the even pages - it starts of at page 2. The images in the top right corner are ALIS data of HF enhanced airglow from Silkimuotka and the images in the lower left corner is the projected images of the retrieved volume emission. The animation starts at 17:40:15 UT with 10 s between each frame and covers about one HF-pump on-off cycle.

next up previous contents
Next: Summary Up: Tomography Previous: Error sensitivity   Contents

copyright Björn Gustavsson 2000-10-24