next up previous contents
Next: 3.2 Commanding Up: 3.1.5 Bias Determination Previous: 3.1.5 Bias Determination Bias Computation Algorithms


Algorithms Requiring No Additional Storage

Ford and Somigliana at MIT have investigated a class of algorithms that possesses the considerable advantage that they use no additional FEP memory beyond the Image and Bias Map Buffers and a modest amount of Data Cache. At any moment during the course of the algorithm, each pixel is represented by only two quantities - the 12-bit value in the Image map which will be overwritten by each fresh exposure, and the 12-bit value in the Bias Map whose value must converge to the desired bias threshold level.

The distribution of values, pi, of a pixel observed N times ($1
\leq i \leq N$) may be approximated by a narrow Gaussian, $\exp[-(p_{i}
- p_{0})^{2} / \sigma^{2}]$, with the addition of a few outlying values, especially those with $p_{i} \gg p_{0}$ produced by X-rays or ionizing particles. By definition, p0 is the modal value and $\sigma (\ll p_{0})$ is the width, typically a few data units. For new CCDs, not yet subjected to radiation damage, both p0 and $\sigma$are nearly identical from pixel to pixel, except for a few pathological `hot' or `flickering' pixels.

When a CCD is damaged by radiation, it is anticipated that p0 and $\sigma$ will increase with time, and at rates that vary randomly from pixel to pixel. Such behavior has been reported by Andy Rasmussen from a study of the CCDs aboard the ASCA observatory. Algorithms that use no additional storage are unable to estimate $\sigma$ from pi, so they concentrate on the modal value p0. Once this has been computed for all pixels, its variance $(\langle p_{0}^{2}\rangle/N) -
\langle p_{0} \rangle^{2}$ serves as a good approximation to $\langle
\sigma^{2} \rangle$.

The algorithms start by examining a series of exposures to `condition' the pi values and derive bi, an estimate of p0. After the first exposure, the best estimate of b1 is clearly p1. After the i'th exposure, two possible algorithms yield improved values of bi, viz.  
b_{i} = \min(b_{i-1}, p_{i})\end{displaymath} (5)
which guarantees that, after a number of exposures, none of the bi will retain anomalously high values from X-ray or ionizing particle events. The resulting bi will generally, of course, be less than p0 by a few $\sigma$, and it will be seriously compromised if even a few of the pi possess anomalously low values. A rather more accurate conditioning can be achieved by the following two-step algorithm:  
b_{i} = p_{i} \mbox{ if } p_{i} < (b_{i-1} - T_{1})\end{displaymath} (6)
b_{i} = Cp_{i} + (1 - C)b_{i-1} \mbox{ if } p_{i} < (b_{i-1} + T_{2})\end{displaymath} (7)
i.e. if the new value pi is much less than the previously stored value bi, let it replace the stored value. Otherwise, if the new value isn't much larger than the stored value, use a running average to better approximate the `mean' pixel value, bi. Optimum choice of the thresholds T1 and T2, and of the partitioning factor C, will depend on the expected pixel distribution functions themselves.

After N conditioning exposures, the Bias Map Buffer contains values bN that are guaranteed to lie within a few $\sigma$'s of the modal values p0. These can now be used as rejection thresholds to identify X-ray and ionizing particle events in subsequent calibration exposures, and thereby improve their own values in the process. This is achieved by making a further set of M exposures, i.e. 0 < N < i < (N+M). The new pixel values pi are subjected to the following two-step algorithm:
p_{i} = 4095 \mbox{ if } p_{i} \gt (b_{i-1} + T_{3}) \end{displaymath} (8)
where T3 is a threshold value. In addition, the 8 neighbors of any pixel that meets this criterion may also receive a contribution from the event that caused the central pixel to lie above the threshold. These neighboring pixels are therefore also reset to a value of 4095.

In the second step, the Image Map pixels are re-examined and those with values less than 4095 are used to refine the Bias Map values,

bi = Cpi + (1 - C)bi-1 (9)

Although this algorithm guarantees that bi will converge to the neighborhood of p0 for all positive C < 1, bi will eventually execute a random walk with width $\sim C \sigma$. In practice, the choice of C must be balanced against the exposure count M to optimize the width of the bN+M distribution versus total calibration time $t \propto M + N$.

The principal advantage of this class of algorithms is that they are very simple to implement, requiring just a few lines of computer code. Preliminary tests show that the resulting bias maps compare favorably with those generated from `strip-by-strip' algorithms. However, they do have some weaknesses, as follows:

MEAN-I2L2 This algorithm calculates the mean and $\sigma^{2}$of N values of a pixel. It then rejects values that diverge from the mean by more than $2 \sigma^{2}$, and then recalculates the mean and $\sigma^{2}$.

\mbox{Definition: } \overline{x} = \frac{1}{N} \sum_{j = 1}^{N} x_{j}\end{displaymath} (10)
\sigma^{2} = \frac{1}{N - 1} \left( \sum_{j = 1}^{N} x^{2} - N\overline{x}^{2}
 \right)\end{displaymath} (11)
\mbox{Rejection rule: } (x_{j} - \overline{x})^{2} \geq k^{2}\sigma^{2}\end{displaymath} (12)
The algorithm is executed by successive iterations. The mean and $\sigma^{2}$ are calculated, then the outliers are rejected and the mean and $\sigma^{2}$ are recalculated. The algorithm is repeated until the termination rule is satisfied. The value k can be thought of as an estimated value derived by Gaussian distributions; e.g. k = 2 implies a probability of 5% that good pixel values are being rejected. The termination rules are defined as any of the following:

The algorithm is expected to work best with Gaussian or symmetric heavy-tailed distributions in the presence of outlier points for which the assumption of normality does not hold. These outlier points are rejected and the remaining data are treated as Gaussian.

The calculation is accelerated by saving the intermediate values $\sum
x_{j}$ and $\sum x_{j}^{2}$. During the rejection phase, rejected values are subtracted from $\sum
x_{j}$ and $\sum x_{j}^{2}$ and the sample number is updated.

MEDIAN The algorithm calculates the median of N pixel values. It sorts them in ascending order, and identifies the central value, thereby automatically rejecting the outliers which will be sorted to one end of the list or the other. Once the {xn} are sorted, the formula for the median is:

x_{med} = \left\{ \begin{array}
 x_{\frac{(N + 1)}{2}} ...
 ...\frac{N}{2} + 1 \right]}) & 
 \mbox{N even} \end{array} \right.\end{displaymath} (13)
A suitable rejection rule is based on the interquartile range (IQR), which indicates the distance between the upper and lower quartile:
x \leq x_{med} - k \mbox{ IQR and } x \geq x_{med} + k \mbox{ IQR}\end{displaymath} (14)

A value of k = 3 can be considered to be a good approximation to reject severe outliers. The termination rule for the median is defined as follows:

The algorithm is expected to work best with Gaussian or symmetric heavy-tailed distributions with few outlier points which have large $\sigma$, but zero mean. It fails if the area in the tails is large (Press, William H., et al., Numerical Recipes in C), and is somewhat influenced by points with very large values, but the rejection of contaminated sources is particularly simple.

MEDMEAN This algorithm calculates the median value after excluding some of the highest and lowest values, and excluding values based on the variance about a trial mean. Specifically, it takes the N values for a specific pixel pi, i=0,N-1, where N is the number of conditioning exposures, and removes the L largest values and the M smallest values. Using the remaining N-L-M values, it computes the median value, $\vert p \vert$, and the variance, $\sigma^{2}$, of those pi less than $\vert p \vert$:
\sigma ^2 = \frac{2}{N - L - M - 1} \left. \sum_{i=0}^{N-L-M...
 ... p_i - \vert p \vert) ^2 \right\vert _{p_i \leq \vert p \vert} \end{displaymath} (15)

The trial set of pi are compared to the trial $\vert p \vert$,and any values which do not fall within $m \sigma$ are removed, i.e., any pi which do not satisfy the condition:  
\left\vert p_i - \vert p \vert \right\vert \leq ( \sigma \cdot m )\end{displaymath} (16)
where m, L, and M are all adjustable parameters, which will be set by AXAF operations and the ASC, and will be infrequently changed.

Finally, the bias value is computed based on the remaining pi, by the normal equation for the mean:
\bar{p} = \frac{1}{N^*-L-M} \sum_{i=0}^{N^*-L-M-1} p_i\end{displaymath} (17)
where N* is the number of pi values remaining after the exclusion for eq. 3.14.

next up previous contents
Next: 3.2 Commanding Up: 3.1.5 Bias Determination Previous: 3.1.5 Bias Determination

John Nousek