A brief note on non-linearity correction

Document Number: VDF-TRE-IOA-00008-0003 Version 1

Mike Irwin, CASU

In which some options for non-linearity correction and its measurement are explored.

Correcting for non-linearity

In default DSC reset-read-read (RRR) mode, downstream of the data aquisition system (DAS) the output that we see is

ΔI′ = I′_{2} − I′_{1} = f(I_{2}) − f(I_{1})

where I′_{1} and I′_{2} denote the non-linear first (ie. the reset-frame) and second readouts respectively and I_{1} and I_{2} the desired linear quantities. The non-linear function f(I) maps the distortion of the desired linear counts to the non-linear system I′. If we define the inverse transform g(I′) that maps measured counts I′ to linearised counts I as the inverse operator g() = f^{−1}() then

I = g(I′) and I_{1} = g(I′_{1}) I_{2} = g(I′_{2})

If I′_{1} and I′_{2} were directly available this is a one-to-one mapping and can be done efficiently and accurately using Look Up Tables (LUTs). This is the conventional way of implementing the correction prior to other image manipulation operations.

However, if I′_{1} and I′_{2} are not separately available and all we have to work from is the difference ΔI′ then a simple LUT transformation is not possible.

For example, taking the simplest case where the illumination level across the detector has not changed during the course of the RRR and no on-board co-addition is happening then, in principle given only ΔI and knowledge of the timing of the RRR operations we can deduce I_{1} and I_{2} by using the effective integration time for each to estimate their scaling to the measured difference ΔI such that,

I_{1} = kΔI and I_{2} = (1 + k)ΔI

Unfortunately, the ratio k will not be constant for the non-linear quantities I′_{1} and I′_{2} forcing us to adopt a scheme along the following lines.

Given ΔI′ and defining the non-linear operator f() as a polynomial with coefficients a_{m} (typically up to quartic) we have

ΔI′ = ∑_{m} a_{m}(I_{2}^{m} − I_{1}^{m}) = ∑_{m} a_{m}[(1 + k)^{m}ΔI^{m} − k^{m}ΔI^{m}]

The quantity we want ΔI is buried in the non-linearity of the RHS and we have to solve an equation like this for every pixel. This is possible, and relatively simple to program using something like a Gauss-Seidel iterative scheme, but is more inefficient than a direct mapping. ^{1}

If we wanted to use a completely general LUT approach we would require a 2D LUT for all possible values of I_{1} and I_{2} ie. 65k × 65k in size, or 4.3×2 Gbytes. Most likely we would need a different correction for each “channel” making a total of 128 (WFCAM) or 256 (VISTA) × 8.6 Gbytes = 1.1 or 2.2 Tbytes of LUT for the cameras! Of course if the range of values of k is limited via exposure time quantisation this decreases the size of the total no. of LUTs required considerably for the constant illumination case, but would be an ugly and possibly impractical solution.

Practical considerations (eg. data volume), suggest two alternative solutions for nonlinearity correction: either correct the individual frames directly in the DAS by measuring and downloading the appropriate LUTs, or polynomial coefficients, to the DAS; or use a non-linear inversion on the reset-corrected frames as outlined here. This methodology is not generally applicable, eg. to multi-NDR/gradient fitting readouts, but is directly applicable to coadded (or coaveraged) frames of the same exposure times, assuming constant illumination over the series.

For reset-corrected data, the non-linear inversion is competitive with complex operations on LUTs and much simpler to implement. It also has the added advantage of removing all aspects of the non-linearity correction from the DAS. The main disadvantages are the method is restricted to DCS RRR mode, and if the illumination level is rapidly varying (eg. twilight) the effective scale factors k_{i} may be hard to compute accurately - although for all realistic practical situations the knock-on effect is likely to be negligible.

Measuring the non-linearity

If all that is available are reset-corrected data from say a time series of dome flats, it is still feasible to directly compute the non-linearity coefficients.

Given a series of measurements {i} of ΔI′_{i} and using the previous notation and polynomial model

ΔI′_{i} = ∑_{m} a_{m}(I_{2}^{m} − I_{1}^{m}) = ∑_{m} a_{m}ΔI_{i}^{m} [(1 + k_{i})^{m} − k_{i}^{m}]

where k_{i} are the exposure ratios under the constant illumination. In general ΔI_{i} = s t_{i} where t_{i} is the exposure time of the ith reset-corrected frame in the series and s is a fixed (for the series) unknown scale factor. The k_{i} are computable from a knowledge of the exposure times and the reset-read overhead, t_{i} and ΔI′_{i} are measured quantities leaving the polynomial coefficients a_{m} and the scaling s to be determined.

1A simulation on a twin processor 3GHz Xeon PC indicates, for typical values, a time of ≈1/8th second per 2k × 2k detector.

Thus the model is defined by

ΔI′_{i} = ∑_{m} a_{m}(I_{2}^{m} − I_{1}^{m}) = ∑_{m} a_{m} s^{m} t_{i}^{m} [(1 + k_{i})^{m} − k_{i}^{m}]

and can be readily solved by standard linear least-squares methods using the following sleight-of-hand. Since the scaling s and hence the polynomial solution a_{m} are coupled, by simply (and logically) requiring in the final solution a_{1} = 1, computation of s can be completely avoided.

Rewriting the previous equation in the following form makes this more apparent

ΔI′_{i} = ∑_{m}(a_{m} s^{m}) t_{i}^{m} [(1 + k_{i})^{m} − k_{i}^{m}] = ∑_{m} b_{m} t_{i}^{m} [(1 + k_{i})^{m} − k_{i}^{m}]

where now b_{m} are the coefficients to be solved for. The final step is to note that

a_{m} = b_{m}/s^{m} = b_{m}/b_{1}^{m}

since by definition a_{1} = 1.

Individual frames available

If both reset and data frames are available then the problem is much simpler since it can be phrased in the following way

I_{1} = ∑_{m} a_{m} I′_{1} ^{m} = s t_{1}

I_{2} = ∑_{m} a_{m} I′_{2} ^{m} = s t_{2}

where t_{1} and t_{2} are the known effective exposure times. Given a series of measurements {i} these can either be solved as a coupled pair of linear least squares minimisations or by forming

Δt = ∑_{m} a_{m}/s (I′_{2} ^{m} − I′_{1} ^{m})

as a standard linear least-squares problem. Again the requirement that a_{1} = 1 automatically defines s. These coefficients can then be directly applied to the observed flux to give the desired linearised units.

Summary

There are two general options for non-linearity correction. The first is to compute and apply the non-linearity correction on individual readouts, preferably in the detector controller system, for data volume considerations, prior to computing the resetcorrected frames. The second method forgoes direct access to the individual frames and instead computes and applies the non-linearity correction at the front-end of the VDFS pipeline.

Each option has certain advantages and disadvantages. The former method is more complex from an operational point-of-view since the non-linearity measures will need to be computed using a different readout mode to that in normal use and then fed back to the controlling system for application. Leaving the non-linearity operation for the pipeline simplifies the operational aspects but makes various (not unreasonable) assumptions about the timing of reset and readout operations and the stability of the illumination.

The only clear drawback of correcting post-facto is that it would probably be very difficult, if not impossible, to correct multiple NDR gradient fitting output that way. It is not obvious that this is a problem for VISTA, or WFCAM, since an adequate alternative for narrow-band imaging is simply to integrate for long enough to get a signal dominated by sky background and stick with the DCS RRR mode.