next previous
Up: Optimizing ISOCAM data processing redundancy


Subsections

5 New approach for LW-ISOCAM data processing

 

5.1 General description

5.1.1 General inversion method

To address the problems described in Sect. 4 we have developed a method which uses the fact that a position on the sky has been observed by several ISOCAM pixels at different times (see Fig. 1). The redundant information is used to separate the contributions of the sky emission and of instrumental effects to the observed signal. This approach has already been used for the IRAS mission (e.g. Okumura 1991; Wheelock et al. 1997) and may be generalized to every raster type observations, whatever the wavelength of observation. Formally the processing of astronomical data with spatial redundancy could be treated as an inversion problem. We can consider that the data observed O is the result of the convolution of the real sky S with the instrumental function I plus some additive noise N:

\begin{displaymath}O = I \otimes S + N.
\end{displaymath} (1)

When the observation has been conducted with spatial redundancy, it may be possible to address the data processing problem globally by finding I, S and N that will reproduce O and minimize the difference between measurements obtained at the same sky position. Such an approach supposes that the main instrumental effects that affect I and Nare known enough to constrain the inversion method. Presently this is not the case for ISOCAM.

5.1.2 Sequential approach

ISOCAM's response variations are complex and presently we are not able to model it with a reasonable number of parameters which would allow us to address the data processing problem by a global inversion. We have thus adopted a sequential approach where the instrumental effects are treated one at a time. Nevertheless, even if we cannot use such a global method, the idea of minimizing the differences between pixels that have seen the same sky position is the mainstay of our approach. In particular, the ISOCAM long term drift problem is addressed by a least square minimization technique based on the fact that, if the detector is slowly reaching stabilization, the measured flux is also approaching the real observed flux. Here we suppose that every pixels of the array is affected by the same long term drift. The implication of this assumption will be discussed in Sect. 6.


  \begin{figure}
\epsfxsize=9cm
\epsfbox{ds8759f4.eps}
\end{figure} Figure 4: Data reduction chain illustrating the data processing steps from left to right. The time dependent flat field and bad pixel identification operations can be done iteratively

Once the LTT is removed, the other response variations (slow glitches, ghosts and residual transient effects) are corrected by comparing the data cube to the sky image. This is done in two steps. First we compute a variable flat field that takes into account pixel-to-pixel response variations. Second, pixels that have a flux that departs significantly from the sky image are flagged (these are called bad pixels in the following). This last operation removes slow glitches and ghosts. These two operations (variable flat field and bad pixel identification) can be done several times to improve the sky image (see Fig. 4).

   
5.2 Long term transient

5.2.1 Gain or Offset effect?

First of all it is important to determine how the long term transient, the flat field F and the observed flux $I_{\rm obs}$ are related to the incident flux $I_{\rm sky}$. We have considered the three following possibilities for the LTT:

1.
Gain effect: $I_{\rm obs} = G(t) F I_{\rm sky}$;
2.
Offset effect affected by flat: $I_{\rm obs} = F (I_{\rm sky} + \Delta(t))$;
3.
Pure offset effect: $I_{\rm obs} = F I_{\rm sky} + \Delta(t)$.

The two GRB observations are of great help in identifying the nature of the LTT. As the two observations were done in exactly the same way, it is possible to subtract directly the two data cubes. In Fig. 5A we show the flux history of two pixels of the GRB2 observation. As one can see there is no long term drift detectable in this observation. The flux difference between the two selected pixels is due to a $\sim$25% flat field difference. The difference between observation 1 and 2 for these two pixels is shown in Fig. 5B. Both pixels have a very similar drift; the 25% difference between the two pixels flux does not appear in Fig. 5B. We then conclude that the LTT does not depend on the signal $FI_{\rm sky}$. In this context, the only valid description of the LTT is the pure offset effect: $I_{\rm obs} = F I_{\rm sky} + \Delta(t)$. In the following we consider the LTT as a single offset over all the detector; in other words we suppose that all pixels are affected by the same offset.


  \begin{figure}
\epsfxsize=9cm
\epsfbox{ds8759f5.eps}
\end{figure} Figure 5: A) Flux history of two pixels of the GRB2 observation. This observation is not affected by the LTT. The flux difference between the two pixels is due to an intrinsic response difference (flat field). B) Difference between the two GRB observations for the two pixels considered in A). Even if the two pixels have very different responses, the LTT amplitude is the same

5.2.2 Formalism

 

For a given pixel at a position (x,y) on the detector array and at a given time t, the observed flux $I_{\rm obs}(x,y,t)$ is related to the temporally varying flat field F(x,y,t), the incident flux $I_{\rm sky}(x,y,t)$ and the long term drift $\Delta(t)$ by the following equation:

 \begin{displaymath}I_{\rm obs}(x,y,t) = F(x,y,t)I_{\rm sky}(x,y,t) + \Delta(t).
\end{displaymath} (2)

Here we suppose that $\Delta(t)$ is not pixel dependent.

The offset function $\Delta(t)$ is found using Eq. (2) and the spatial redundancy inherent to raster mode observations. We determine $\Delta(t)$ by solving a set of linear equations obtained by comparing flat field corrected intensities of the same sky positions but obtained at different times. The flat field F(x,y,t) is computed using the adaptive flat field method (see Appendix A.2). The flat field corrected intensities can be written:

 \begin{displaymath}\frac{I_{\rm obs}(x,y,t)}{F(x,y,t)} = I_{\rm sky}(x,y,t) + \frac{\Delta(t)}{F(x,y,t)}\cdot
\end{displaymath} (3)

In raster mode the camera does not always point at the same position on the sky. To compare two pixels that have seen the same position on the sky, we must put all readouts in a common spatial reference frame. Practically, we project $I_{\rm obs}(x,y,t)$ and F(x,y,t) on the plane of the sky (right ascension ($\alpha$) - declination ($\delta$) reference frame). In this reference frame and following Eq. (3), the difference between two pixels that have seen the same sky $I_{\rm sky}(\alpha, \delta)$at different time (ti and tj) is:

 \begin{displaymath}\frac{I_{\rm obs}(\alpha, \delta, t_i)}{F(\alpha, \delta, t_i...
...elta, t_i)} -
\frac{\Delta(t_j)}{F(\alpha, \delta, t_j)}\cdot
\end{displaymath} (4)

To ensure that we compare fluxes obtained at the same sky position, the field of view distortion must be accurately taken into account in the projection operation.

The function of interest $\Delta(t)$ is estimated using a least-square minimization technique with the following minimization criterion:

 
$\displaystyle \chi^2$ = $\displaystyle \sum_{\alpha, \delta, t_i, t_j} \left[\frac{I_{\rm obs}(\alpha, \...
... t_i)}
- \frac{I_{\rm obs}(\alpha, \delta, t_j)}{F(\alpha, \delta, t_j)}\right.$  
    $\displaystyle - \left.\frac{\Delta(t_i)}{F(\alpha, \delta, t_i)} +
\frac{\Delta(t_j)}{F(\alpha, \delta, t_j)} \right]^{2}\cdot$ (5)

Here the sum is over all the possible pixel pairs that have seen the same region of the sky.

The function $\Delta(t)$ which minimizes the value of $\chi^2$ is found by solving the system determined by:

 \begin{displaymath}
\frac{\partial\chi^2}{\partial\Delta(t_i)} = 0.
\end{displaymath} (6)

Equation (6) represents a standard set of linear equations which can be written in a matrix form: $A \Delta(t) = B$, with:

\begin{displaymath}A(i,j)_{i\neq j} = \sum_{\alpha, \delta}\frac{1}{F(\alpha, \delta,
t_{i})F(\alpha, \delta, t_{j})}
\end{displaymath} (7)


\begin{displaymath}A(i,j)_{i=j} = \sum_{\alpha, \delta, t_{j\cap i}}\frac{-1}{F(\alpha, \delta,
t_{i})F(\alpha, \delta, t_{j})}
\end{displaymath} (8)


\begin{displaymath}B(i)\! =\!\!\! \sum_{\alpha, \delta, t_{j\cap i}}\!\! \frac{1...
...obs}(\alpha, \delta, t_i)}{F(\alpha, \delta, t_i)}\right]\cdot
\end{displaymath} (9)

The A(i,j) and B(i) terms are computed with all pixels that overlap. For each pair (i,j), the sums are always computed on the positions ( $\alpha, \delta$) where $I_{\rm obs}(\alpha, \delta, t_i)$ and $I_{\rm obs}(\alpha, \delta, t_j)$ are both defined.

Finally $\Delta(t)$ is found from: $\Delta(t) = A^{-1}B$. As the second derivative of $\chi^2$ is always positive, the solution found for $\Delta(t)$ necessarily minimizes the $\chi^2$ criterion. We add an offset to all values of $\Delta(t)$ to force the correction to be zero at the end of the observation: $\Delta(t_{\rm end}) = 0$. This is justified by the fact that the long term transient tends to stabilize at the end of an observation. However, there are several observations not stabilized at the end, so that the absolute brightness may be systematically shifted (generally below a few %). For the observations of low contrasted clouds on top of a flat large scale emission (at least the zodiacal emission), this point is not critical as we always remove the large scale emission at the end of our processing.

5.2.3 Practical implementation

 

Other than detector noise, two additional sources of uncertainties affect the comparison of raster images and thus the LTT correction through the minimization algorithm described in the previous section: slow glitches and flat field variations along the observation. The signal measured on point sources is not identical in the individual readouts, essentially due to the undersampling of the the point spread function for the 6'' pixels. Bright sources are thus an additional source of error but they can easily be discarded. The two other noise sources make the practical implementation of the formalism a non-straightforward procedure.

An extensive use of the LTT inversion method presented here on real data has demonstrated the extreme importance of an accurate flat field. As the LTT correction is based on the comparison of the brightness measured by different parts of the array, its result depends on the accuracy of the flat field. In particular, when $F(\alpha,\delta,t)$ is not well estimated, oscillations in phase with the rastering of the observations may appear in the correction found. To overcome this difficulty we compute an approximate LTT correction (as oppose to the exact correction which is the solution of Eq. (6)). The study of low contrast ISOCAM data affected by long term drifts has shown that, in some cases, the LTT can be approximated by the sum of two exponential functions, one upward and one downward:

\begin{displaymath}\Delta(t) \approx P \times \exp{\left(- Q \times t^R\right)} -
S \times \exp{\left(- T \times t^U\right)}
\end{displaymath} (10)

where P, Q, R, S, T and U are strictly positive. By approximating the LTT in this way, we greatly simplify the problem as $\Delta(t)$ now depends only on six parameters. As for the exact solution, the approximated solution is found by minimizing the criterion of Eq. (5), using an IDL curve fitting program (MPFIT). The obvious advantage of this method is that the approximated LTT correction is smooth, getting rid of the oscillations found in some exact solutions.

On the other hand, we have observed cases where this approximation of the LTT does not hold. Sometimes, the LTT shows oscillations that can dominate the emission in low contrasted fields (e.g. cirrus clouds - see Miville-Deschênes et al. 2000). In these cases, an accurate estimation of the flat field and the use of the exact solution for the LTT are mandatory.

5.2.4 LTT correction for the GRB1 observation

Figure 6 presents the LTT corrections found for the GRB1 observation. A adaptive single flat field (see Appendix A.2) was used and point sources were discarded to compute the exact and approximated corrections. The approximated LTT correction is smoother than the exact one, which oscillates with a $\sim 0.2$ ADU/g/s amplitude and a period of $\sim 1800$ s. In this case we have applied the approximated LTT correction. The sky image computed after that correction has been done is presented in Fig. 8B. One sees that it is mandatory to apply the LTT correction since its amplitude (3 ADU/g/s) is nearly three times the amplitude of the emission of interest (1.1 ADU/g/s). At this stage we have used a single flat field to compute the sky image. It is necessary to use a variable flat field to correct the artefacts (e.g. periodic patterns) seen in Fig. 8B.


  \begin{figure}
\epsfxsize=9cm
\epsfbox{ds8759f6.eps}
\end{figure} Figure 6: Exact and approximated LTT correction for the GRB1 observation. Estimating the error on the LTT correction is difficult as it is very sensitive to the accuracy of the flat field used. The error bar represented here is the statistical error (0.08 ADU/g/s - see Appendix B). A more realistic error determination would be related to the accuracy of the flat field. For the GRB observations, a 1% uncertainty on the flat field corresponds to a 0.4 ADU/g/s error on the LTT correction

   
5.3 Variable Flat Field

After the LTT has been corrected, we then take into account the pixel-to-pixel temporal variations of the detector response. These response variations are observed at various timescales. At short timescales they are due to bad short term transient correction, particularly on point sources. On longer time scales they are mainly caused by slow glitches (see Sect. 4). The use of a single flat field (see Sect. 3) does not take into account these temporal variations (that represent $\sim$1-3% of the average flat field) preventing the benefit of the optimal ISOCAM sensitivity. To go further in the data processing, we try to correct these pixel-to-pixel response variations with a time-dependent flat-field (called "variable flat field'' in the following).


  \begin{figure}
\epsfxsize=9cm
\epsfbox{ds8759f7.eps}
\end{figure} Figure 7: Ghost identification. Flux history of a pixel which has observed a point source. The flux values after the point source are affected by a memory effect and are identified as a ghost (flagged out)


  \begin{figure}
\includegraphics[width=18.5cm]{ds8759f8.eps}
\end{figure} Figure 8: LW10 images of a the first GRB 970402 observation after deglitching, dark subtraction, transient correction A), long term transient correction B), variable flat field C) and bad pixel removal D). Image E) is the final map of the second GRB 970402 observation and image F) is the difference between D) and E). For these observations, 1 ADU/g/s corresponds to 0.242 mJy/pix or 0.286 MJy/sr

For LTT corrected data, the observed flux $I_{\rm obs}$ at position (x,y) on the array and at time t is

 \begin{displaymath}I_{\rm obs}(x,y,t) = F(x,y,t) I_{\rm sky}(x,y,t),
\end{displaymath} (11)

where F(x,y,t) is the instantaneous response of pixel (x,y) at time t. Thus flat field and sky structures are mixed together in $I_{\rm obs}(x,y,t)$. In the following we show how the flat field variations can be extracted from the data by estimating $I_{\rm sky}(x,y,t)$ and by taking advantage of the spatial redundancy.

 

In raster observation mode, many pixels of the data cube have seen the same position ( $\alpha, \delta$) on the sky. By averaging all these pixels we reduce the noise due to instrumental effects on the computation of the sky image $I_{\rm sky}(\alpha, \delta)$. The estimate of $I_{\rm sky}(x,y,t)$ is made by an inverse projection of $I_{\rm sky}(\alpha, \delta)$ on each readout of the data cube. Then F(x,y,t) is computed by averaging $I_{\rm obs}(x,y,t)/I_{\rm sky}(x,y,t)$ with a sliding window on the time axis (see Appendix A.1). Here are the guidelines of this method:

1.
Construct a sky image;
2.
Smooth (median smoothing) the sky image with a $10\times10$ window;
3.
Compute an ideal cube $I_{\rm sky}(x,y,t)$ by projecting the smoothed sky image on each readout of the data cube;
4.
Smooth (median smoothing) $I_{\rm obs}(x,y,t)/I_{\rm sky}(x,y,t)$ on the time axis. The size of the smoothing window should be of the order of the time spent on 5 different sky positions.

The sky image of the GRB1 observation, obtained with the variable flat field, is shown in Fig. 8C. The size of the filtering window on the time axis is 100 (corresponding to $\sim$500 s). As seen in Fig. 8C, the variable flat field removes almost all periodic patterns due to high-frequency variations of the detector response. Compared to other methods we have tested (see Appendices A.1 and A.2), this variable flat field gives by far the best results.

5.4 Bad pixel identification

The deglitching process used at the beginning of the reduction chain (see Sect. 3) removes extremely deviant pixels that have been hit by cosmic rays. To go further in the noise minimization process, we again take advantage of the spatial redundancy.

5.4.1 Ghost removal

Memory effects often appear after a strong flux step (e.g. after a point source) due to improper short term transient correction; this is what we call ghosts. To identify pixels affected by such effects, we look at the flux history of every pixel for a memory effect after a flux step (see Fig. 7). Again, we use the redundancy information to improve the identification of ghosts. Here is how we proceed:

1.
Smooth (median smoothing) the sky image with a 3 $\times$ 3 window;
2.
Compute an ideal cube by projecting the smoothed sky image on each readout of the data cube;
3.
Compute the residual cube = |data cube - ideal cube|;
4.
Identify strong flux steps as residual flux above 10 $\times$ the noise level;
5.
On the time history of each pixel, identify a ghost as a residual flux above or under the noise level after a strong flux step;
6.
Flag ghosts in the original data cube.


  \begin{figure}
\epsfxsize=8.3cm
\epsfbox{ds8759f9.eps}
\end{figure} Figure 9: Structure function of the difference between the two GRB final maps (bottom line) and of the difference between the first map of GRB1 and the final map of GRB2 (top line)

5.4.2 Noise reducer

We can go further in the reduction of the noise level by working sky position by sky position instead of working on the time history of every pixel. The idea is to look at pixels in the data cube that have seen the same sky position and discard deviant flux values.

First the sky image is smoothed (median smoothing) with a 10 $\times$ 10 window. Then, for a given sky position ($\alpha$, $\delta$), we compare the N pixels in the data cube that have seen that position with the flux of the smoothed sky image. For most of the sky position, the N pixels are distributed around the smoothed sky estimate. On the other hand, at point source positions, the N pixels generally will be above the smoothed sky estimate. Furthermore, it is also possible to find sky positions where most of the N pixels have fluxes under the smoothed sky level. Here is how we deal with each case:

1.
Most of the N pixels are around the smoothed sky level (SML);
The new sky flux is the average of ( $SML - 3 \times noise < I < SML + 3 \times noise$);
2.
Most of the N pixels are above the smoothed sky level;
The new sky flux is the average of I > SML;
3.
Most of the N pixels are under the smoothed sky level;
The new sky flux is the average of ( $SML - 3 \times noise < I < SML$).
This method allows us to globally reduce the noise level and to keep a good photometry of point sources. This is the final data processing step. The sky images obtained after bad pixel removal of the first and second GRB observations are presented in Figs. 8D and 8E.


next previous
Up: Optimizing ISOCAM data processing redundancy

Copyright The European Southern Observatory (ESO)