3.1.5.1 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, p_i, 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, p₀ 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 p₀ 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 p₀ 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 p_i, so they concentrate on the modal value p₀. 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 p_i values and derive b_i, an estimate of p₀. After the first exposure, the best estimate of b₁ is clearly p₁. After the i'th exposure, two possible algorithms yield improved values of b_i, viz.  

$$b_{i} = \min(b_{i-1}, p_{i})$$ (5)

which guarantees that, after a number of exposures, none of the b_i will retain anomalously high values from X-ray or ionizing particle events. The resulting b_i will generally, of course, be less than p₀ by a few $\sigma$, and it will be seriously compromised if even a few of the p_i 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})$$ (6)

 

$$b_{i} = Cp_{i} + (1 - C)b_{i-1} \mbox{ if } p_{i} < (b_{i-1} + T_{2})$$ (7)

i.e. if the new value p_i is much less than the previously stored value b_i, 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, b_i. Optimum choice of the thresholds T₁ and T₂, 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 b_N that are guaranteed to lie within a few $\sigma$'s of the modal values p₀. 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 p_i are subjected to the following two-step algorithm:

$$p_{i} = 4095 \mbox{ if } p_{i} \gt (b_{i-1} + T_{3})$$ (8)

where T₃ 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,

 

b_i = Cp_i + (1 - C)b_i-1 (9)

Although this algorithm guarantees that b_i will converge to the neighborhood of p₀ for all positive C < 1, b_i 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 b_N+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:

• They converge on p₀, the modal pixel values. This would be satisfactory if the widths $\sigma$ were identical for each pixel over the lifetime of the CCD since a bias threshold value of, say, $p_{0} + 3 \sigma$ would be a robust discriminator of energy deposition. However, this is not even true of undamaged CCDs - the width of flickering pixel values can be many times the average $\sigma$, and the situation is expected to deteriorate as the CCDs absorb high doses of radiation. The choice of bias threshold, p₀ + T_b, must be made carefully - if T_b is too low, flickering pixels will appear as events; if too high, low-energy events will go unreported.
• Since all parts of the final Bias Map are calculated after the last calibration exposure, there is no opportunity to copy the bias values to the downlink telemetry system while the bias is being calculated. The time required to compress and transmit the six Bias Maps, ca. 12 minutes assuming 75% compressibility, must therefore be added to the bias calibration time.
• These algorithms rely on particular properties of the pixel distributions. For instance, if equation 3.3 or 3.4 + 3.5 are used to condition the bias values, the presence of anomalously low pixel values, $p_{i} \ll p_{0}$,will start the second phase of the algorithm in a poorly conditioned state and the number of subsequent exposures, M, may be insufficient for equation 3.7 to change b_i so that $\vert b_{i} - p_{0}\vert \leq C \sigma$, as desired. This situation can be alleviated by applying an additional `grading' process at the end of the conditioning phase in which the anomalously low values of b_N are identified and replaced by better estimates, e.g. by the median of the b_N of the surrounding pixels.

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}$$ (10)

$$\sigma^{2} = \frac{1}{N - 1} \left( \sum_{j = 1}^{N} x^{2} - N\overline{x}^{2} \right)$$ (11)

$$\mbox{Rejection rule: } (x_{j} - \overline{x})^{2} \geq k^{2}\sigma^{2}$$ (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:

• All pixel values are within an acceptance range.
• The maximum number of iterations was reached.
• The number of remaining pixels is less than or equal to the minimum number acceptable.

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 {x_n} are sorted, the formula for the median is:

$$x_{med} = \left\{ \begin{array} {ll} x_{\frac{(N + 1)}{2}} ... ...\frac{N}{2} + 1 \right]}) & \mbox{N even} \end{array} \right.$$ (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}$$ (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:

• All pixel values are within an acceptance range.
• The maximum number of iterations was reached.

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 p_i, 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 p_i 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}$$ (15)

The trial set of p_i are compared to the trial $\vert p \vert$,and any values which do not fall within $m \sigma$ are removed, i.e., any p_i which do not satisfy the condition:  

$$\left\vert p_i - \vert p \vert \right\vert \leq ( \sigma \cdot m )$$ (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 p_i, by the normal equation for the mean:

$$\bar{p} = \frac{1}{N^*-L-M} \sum_{i=0}^{N^*-L-M-1} p_i$$ (17)

where N^* is the number of p_i values remaining after the exclusion for eq. 3.14.