As we start supporting non-linearity correction for more cameras I would like to propose some rules, to keep things simple and allow a default version of non-linaerity correction in IsrTask.
We need to support correction based on lookup tables and have existing code that uses polynomials. We may want to support variations of these or other techniques. Thus I propose to make non-linearity correction a subtask of IsrTask that can be retargeted. I further propose that default task support polynomial correction, because it is simple and all the necessary data exists in the AmpInfo table.
For the polynomial-based correction there are a few issues:
- We have four coefficients available. This can easily be increased if needed, but there is no sign of that, rather the opposite...
- Unused coefficients are often set to nan (e.g. in obs_subaru and obs_test)
- The settings for various cameras have some quirks:
- obs_subaru, which has working non-linearity correction code, sets linearityType= "PROPORTIONAL" (or "NONE" if no correction is wanted) in the AmpInfo table, and the associated linearityCoeffs are as follows. This seems a very unnatural order, and I dislike forcing a value to 0 and raising an exception if it is not:
- coeff sqareCoeff: corrIm = uncorrIm + sqareCoeff*uncorrIm^2
- coeff threshold: this must be 0 else an exception is raised
- coeff maxUncorr: if maxUncorr > 0 then flag pixels as SUSPECT where uncorrIm > maxUncorr
- lsstSim sets linearityType="polynomial" and linearityCoeffs=(0, 1, 0, 0), presumably the first four polynomial coefficients. This data is ignored.
- The concept of maxUncorr from obs_subaru seems useful, but it is not a natural fit to a lookup table.
I propose the following:
- Add field suspectLevel to the AmpInfo table and set the SUSPECT flag using specific code for it, rather than relying on non-linearity correction code to do it. This makes the field more obvious than tucking it away as an element of linearityCoeff and also simplifies writing variant tasks to correct non-linearity, since they can focus on doing one thing well.
- Offer constants for linearityType names and use those in code. This catches typos earlier. Suggested names include NONE, POLYNOMIAL and LOOKUP_TABLE. These can be enums in python 3.
- The standard non-linearity correction tasks must accept NONE and do nothing for such detectors. Mosaic cameras may have detectors for which we don't want this correction.
- The standard non-linearity correction tasks will only support one linearityType or NONE. More heterogeneous cameras will require specialized tasks.
- Encode the number of coefficients by adding an extra field, e.g. numLinearityCoeffs. The specified number of fields are required to be finite. Yes we could (and perhaps should) also force the camera to make all unused fields 0, but this is an extra step and is clearly easy to forget. In any case, it is useful information; for instance it can save time to know the order of a polynomial in advance.
- For the standard polynomial non-linear correction task: the coefficient array starts with the x^2 coefficient and goes up from there. Constant and scale correction are explicitly omitted because they are not useful.
- is triggering
DM-5462 Add non-linearity correction to ISR task
I agree about the truncate-floats-and-lookup argument. But it means that you have to use the form value += correction as otherwise you'd be truncating the answers.
Having digested the above, here is my proposal:
AmpInfoRecord gains a new field suspectLevel. ISR task uses this to flag values as SUSPECT by testing before applying non-linearity correction. If the value is 0 then the flagging step is omitted.
We will support two non-linearity types (and others can be added as wanted):
- table lookup, as described above.
- The second-order correction used by obs_subaru, but with the suspect level specified by AmpInfo.suspectLevel instead of by a coefficient. This will require a minor change to the code that creates subaru's amp info tables and a minor change to Subaru's ISR task.
For lookup table based non-linearity correction: use an array that has 32k or 64k-ish values along the quickly varying axis (one value per possible pixel value, truncated to the nearest int) and at most one entry per amplifier along the slowly varying axis. The non-linearity coefficients in Detector are then:
- index info amplifier axis of the non-linearity table. This allows multiple amplifiers to share the same table (especially useful if there is just one table for all amplifiers). Note: if we are sure we always want one table per amplifier then it would be simpler and less confusing to use the index of the amplifier in the amp info table as the index of the amplifier data in the lookup table.
- an offset
- a scale if folks think it desirable. In theory this allows finer or coarser resolution than one value per ADU, but it slows down lookup a bit and makes the lookup table a bit harder to understand.
What to do about values that are out of range? I believe every value must be checked to avoid mysterious segmentation faults or worse. I prefer to raise an exception but would truncate if folks felt that was better. (A bad amplifier index must raise an exception; truncation makes no sense.)
Pseudocode for table-based correction:
amp_index = int(float_amp_index) # float_amp_index is coeff 0
if amp_index out of range:
raise an exception
value_index = int((non_corrected_value + offset)*scale) # offset is coeff 1, scale is coeff 2
if value_index out of range:
raise an exception or truncate; which?
corrected_value = nonlinearity_array[amp_index, value_index] + non_corrected_value
A registry will be used to construct the code that applies the non-linearity correction, as requested by Paul Price. linearityType will contain the registry key.
Where to store the non-linearity tables? The non-linearity table for LSST will be roughly 4MBytes in size (64k values x 16 amps x 4-bytes per float entry). This seems too large to me to sensibly put into the Detector object, where it would always be in memory when we read any calexp or similar image using the butler. It's even more worrisome that the Camera object has a Detector for every CCD in the focal plane! There are several alternatives:
- Treat non-linearity tables as a separate product (like flat fields) that is read as wanted. This is my suggestion. It is simple and well supported and seems appropriate for such large data.
- Save the information in the Detector object, but use lazy loading and flush it after ISR runs. Detector is a C++ object, so this seems unlikely to be easy.
Having talked to Robert Lupton and Kian-Tat Lim today, here is a modification to the proposal described in the previous comment:
- Non-linearity correction is provided by a Python object stored as a butler dataset type (e.g. using pickle). This object contains both the code and the data needed to correct the amplifiers of one CCD.
- This requires each obs_* package to define a new dataset type for any camera that supports non-linearity correction.
In particular, the ISR task will be rewritten to unpersist the non-linearity corrector when it unpersists the flats and other calibration data, and will pass it along to the main data processing method just like the other calibration products.
I think we have converged. The proposal is as described in my previous two comments.
Paul Price asks "aren't the images in floating-point". Yes, by the time we apply non-linearity correction they are. The plan is to truncate each pixel value, use that to determine the index in the lookup table (e.g. by adding an offset if needed) and adding (or subtracting, if you prefer) the looked-up value to the pixel data. The lookup table will presumably contain floats.