Dark Frame Optimization Algorithm

Juan Conejero

PixInsight Staff
Staff member
We have seen some discussions recently on other forums about dark frame correction and dark frame optimization. Since some of these discussions and analyses involve our image calibration tools and algorithms, and mainly because I've read some uninformed opinions and judgements, I think it is pertinent to stick here a copy of a post I wrote back in 2012, where I describe our current dark optimization algorithm, how it works, and its main limitations.

That said, an improved dark frame optimization routine is currently in the implementation/test phase and will be released in the coming Fall/Winter for PixInsight 1.8.4. The new algorithm is more robust to read noise in the master dark frame and improves on the hot pixel undercorrection problem.

======================================

Let's state the problem first. We may describe an uncalibrated raw light frame with this very simple equation:

I = I0 + D + B​

where I is the uncalibrated raw light frame, D is the dark current, and B is the bias pedestal. All terms in this equation are vectors (or matrices) because we are describing a process for a whole image at the pixel level. I0 is the dark-and-bias-subtracted raw light frame (we know that the next calibration step would be flat fielding, but this is unimportant to this discussion). Our goal here is to obtain I0.

As we have stated it above, this is a trivial problem: just subtract good master bias and dark frames from the uncalibrated light frame, and the result will be a good approximation to I0 (the better the masters, the better the approximation). However, the actual problem is more complex because the dark current D is highly time- and temperature-dependent, so a perfect match between the dark current signals in I and D is only a theoretical goal. A better (but still far from perfect) approximation to the actual problem is more like this:

I = I0 + k*D + B​

where k is a scaling factor that attempts to account for effective exposure time and temperature differences between the light frame and the master dark frame. The process of finding a good approximation to the scaling factor k is what we call dark frame optimization.

Note that I have said that this is far from a perfect solution. It is wrong mainly because the dark current does not vary linearly with exposure time and temperature for the whole numeric range of the data, so a single multiplicative factor cannot be very accurate, especially when time and temperature differences between I and D are large. Dark optimization, however, is always better than nothing. From our experience, dark optimization is beneficial even if exposure times and temperatures are matched between light and dark frames. For this reason the corresponding parameter is enabled by default in our ImageCalibration tool.

Now that we have described the problem, let's describe our solution. Algorithmically this is known as an optimization problem: find the value of a parameter that minimizes (or maximizes) the value of a function. For example, think in terms of the economical cost of a production process. We have to produce an item A that depends on a factor r, which leads to the problem: find the value of r that minimizes the cost of producing A. In real cases r usually represents a complex set of parameters and constraints, such as environmental factors, availability of prime matters, production systems, etc., so that the most difficult part of the solution often consists of identifying a significant set of parameters or model to define a suitable cost function. This leads, for example, to the subject of linear programming and, more generally, to mathematical optimization.

The dark frame optimization problem is a relatively simple case of function optimization with a single parameter. The first step is to define the cost function that we want to minimize. To explain why and how we have selected a particular function we need some visual analysis. Take a look at the following two images:

A.png
B.png


One of them is a crop of a bias-subtracted master dark frame, the other is not. Difficult to say which is which, isn't it? The image to the right is a mix of synthetically generated uniform noise and impulsional noise (salt and pepper noise) with 0.5% probability. It has been generated with the NoiseGenerator tool in PixInsight. With this comparison I am trying to show that a master dark frame looks very similar to random noise: essentially, a master dark frame is composed of pixel-to-pixel intensity variations whose spatial distribution is rather uniform, plus a relatively small amount of hot pixels whose distribution and typical values are similar to impulsional noise. Of course, we know that the thermal noise is a fixed pattern characteristic of a given sensor, so it is not random because it is predictable. However, morphologically a master dark frame is virtually indistinguishable from random noise. It is not random noise, but it behaves like that, and we are going to take advantage of this simple property to implement a purely numerical solution.

So our cost function is just a noise evaluation function. We already have implemented a powerful multiscale noise evaluation algorithm in several PixInsight tools, such as ImageIntegration, to implement a noise-based image weighting algorithm (Jean-Luc Starck and Fionn Murtagh, Automatic Noise Estimation from the Multiresolution Support, Publications of the Royal Astronomical Society of the Pacific, vol. 110, February 1998, pp. 193-199). In this case, however, we have implemented a simpler and very efficient method (Jean-Luc Starck and Fionn Murtagh, Astronomical Image and Data Analysis, Springer-Verlag, 2002, pp. 37-38) known as k-sigma noise estimation:

- Compute a single-layer wavelet transform of the image. This transform consists of the finest wavelet scale w1 plus a residual cJ, which we simply discard. We use the standard B3 spline wavelet scaling function as a separable low-pass filter.

- Iterate a k-sigma process as follows: At each iteration n > 0, denote by d(n) the subset of pixels in w1 such that |d(n-1)ij| < k*sigma(n-1), where sigma(n-1) is the standard deviation of the current subset of pixels (from the previous iteration). In our implementation we have d(0)=w1, k=3, and this process is iterated until no significant difference is achieved between two consecutive iterations. For dark frame optimization, we iterate up to 1% accuracy and normally 3 or 5 iterations are sufficient.

- The final noise estimate is the standard deviation of the resulting subset of pixels. Note that this is an unscaled noise estimate of Gaussian white noise, since it depends on the wavelet scaling function used. However, since we don't want to compare noise estimates among different images, but only to minimize the noise in a particular image after dark subtraction, an unscaled estimate is perfectly suitable.

Note that the use of a wavelet transform is actually irrelevant in this process; we use it just because it is very fast and we have it readily available on the PixInsight/PCL platform. A simple convolution with a low-pass filter followed by a subtraction from the original image would be equivalent.

This noise evaluation algorithm is robust and 'cannot fail' The rest of our algorithm is a routine to find the minimum of a unimodal function. For the sake of robustness, I haven't implemented anything fancy here: just an adaptation of the classical golden section search algorithm. Perhaps not the fastest alternative, but golden section search is absolutely robust, which is what I want here.

Summarizing, the algorithm can be described at a high level as follows:

- Find a range [k0,k1] of dark scaling factors that bracket the minimum of the noise evaluation function. This is the initial bracketing phase.

- Iterate the golden section search algorithm to find the minimum of the noise evaluation function. At each iteration, compute the calibrated raw frame as:

I0 = I - k*D - B​

and evaluate the noise in I0 with the k-sigma iterative method.

- Iterate until the minimum of the noise evaluation function is determined up to a prescribed accuracy (1/1000 fractional accuracy in our implementation).

- The final value of the dark scaling factor is the value of k from the last iteration.

From our tests, this algorithm has proven extremely robust and accurate. However it has two known problems:

- It is very sensitive to bad quality calibration frames. In particular, if the master bias and/or the master dark frames include significant amounts of read noise, the optimization algorithm will find a value of k that overcorrects thermal noise to compensate for the additive random noise component. In extreme cases, this may lead to dark holes as a result of overcorrected hot pixels. However, if the master calibration frames are bad quality ones, the whole image reduction process makes no sense, actually...

- It implicitly assumes that the dark optimization function is constant for the whole range of numerical data values, since it consists of a simple multiplicative factor. This is not the case in practice, especially for high intensity values. As a result of this simplistic optimization model, the algorithm tends to undercorrect hot pixels. This is a very minor issue, however, since hot pixels are easily rejected during integration (because you dither your images, don't you?) and can be removed also with cosmetic correction techniques.

In a future version of the dark optimization algorithm we'll compute different optimization factors for separate intensity ranges in a spline fashion. This should fix the hot pixel undercorrection problem.
 
I was just about to ask another question re: dark frames and you answered most of my concerns.  You note that the dark noise components are not linear with time and temperature. So my question is this: is there a reasonable way to pick a time/temperature for your typical images (i.e., due to light pollution, I'm usually limited to under 2 or 3 minutes max binned 1x1 and certainly less binned 2x2). And my working temperature range is -10C (hot summer night) to -25C (winter night).  I have a set of dark frames every 5 degC from -30C to -10C at 3 minutes exposure so I presume I can check the linearity of the dark noise vs temperature, then maybe do the same thing at different exposure times, right?  Would that be of any value?  Or would a simple 3 minutes at -15C (kind of in my middle temperature range) work about as well as anything?  I only image for my own enjoyment so extreme precision is not required.

I understand you can't give me a definitive answer, just point me in the right direction.
 
Dear Juan,

When creating darks and bias files with MaximDL, the default pedestal setting is 100ADU.
If these files are then used in Pixinsight to calibrate lights using the default settings, the program removes the pedestal as a first step.  When the resulting calibrated files are then measured with subframeselector, most of the measured values are zero, i.e. the calibrations are unusable.

I just tried to process a set of images of the Bubble Nebula.  I had 10 red, 10 blue, 10 green and 10 Halpha.  After much head-scratching and hair-pulling, I finally worked out that if I increase the output pedestal as 100 ADU, the reds and greens are correctly calibrated.  To calibrate the blue, I needed to increase the output pedestal to 200, whilst to calibrate Halpha, I needed a pedestal of 300 ADU.

Is this behaviour expected?  I ask because I am unable to find any documentation regarding increasing the pedestal in this way.  The nearest I got was your message from 2010 where you imply that only in very special circumstances it may be necessary to add a small pedestal in the range 50 to 200 ADU "especially in master calibration frames".  Perhaps what I am doing is exactly what you recommended, but I am puzzled as to why the need to do this adjustment is not better documented.  I spent hours and hours trying to work out why the lights were not being properly calibrated.  In fact until recently, I gave up and just skipped the calibration altogether.  Even the new excellent book "Inside Pixinsight" by Warren Keller is silent on this question.

Perhaps this issue is peculiar to my imaging setup for some reason, but I don't think so because I had a very similar problem when using a different camera, a different scope and different imaging software (Astroart).

Sorry if I am raising an issue that everyone knows about except me!

Best regards,

Stephen
 
I have the same issue. I have never been able to use Bias frames for a successful calibration in PI. I get zero background values that render the calibration useless. Dark and flat calibration work fine. If anyone could shed some light on this behavior that would be great.

Nick
 
I use modded Canon XSi and T3i cameras.  I am using the trial PI to look at some data processed by ImagesPlus (IP) and PS3. I did not dither, but I took all the lights, darks and flat, within a four hour window. Flats were taken using an EL plate.  IP stacked 15 in the final image whereas PI stacked only 6.

In PI I received many errors about incompatibility of darks during calibration. The same happened with the flats. Then only six lights were stacked 

My question is what to do?  If the darks and flats were taken one after another as rapidly as possible, what more can a person do to obtain darks and flats which won't cause errors?  In a non-cooled DSLR, the temperature will naturally rise. If this is the problem, one might just forget darks and go to dithering.

Finally, does the lack of good darks and flats have anything to do with the rejection of lights in PI?
 
schulzep said:
In PI I received many errors about incompatibility of darks during calibration. The same happened with the flats.

By errors, do you mean "Warning: No correlation between dark and target frames" messages?  These are common when calibrating flats with darks that have much longer duration.  This is just a warning, not a fatal error.  It is telling you that the result of dark calibration was going to be noisier than the original sub, so it skipped the dark calibration.

schulzep said:
Then only six lights were stacked 

My question is what to do?  If the darks and flats were taken one after another as rapidly as possible, what more can a person do to obtain darks and flats which won't cause errors?  In a non-cooled DSLR, the temperature will naturally rise. If this is the problem, one might just forget darks and go to dithering.

Finally, does the lack of good darks and flats have anything to do with the rejection of lights in PI?

If some of your lights were rejected then it is more likely that image registration failed.  Check the output on the process console.  StarAlignment prints a summary with the number of lights successfully registered and the number of failures.  You may need to tweak the registration parameters if this is the case.  I'd also recommend you Blink check the calibrated lights to see if they all look OK.

I'd also strongly recommend dithering whether you take darks or not.  Personally, I always do both!

Cheers,
Rick.
 
That algorithm is a very good idea because I know from experience with DSLRs that it is almost impossible to take dark frames at the same temperature as the images.  (Or images at the same temperature as each other!) Some cameras record the temperature and tell us how much it is shifting; others leave us guessing.
 
I notice that when I use dark frame optimization, I get calibrated images with higher amplitude than when I don't.  Does this mean the optimization checkbox (in Batch Preprocessing) also does some scaling of the calibrated images?
 
Juan, thank you for this very interesting article. It has answered a question that has been troubling me, namely: why is there no EXPTIME keyword in the FITS header of my master darks? I thought this was needed in order to scale the master dark to the same exposure as the light frames (or flats) being calibrated.

Now I see that PI uses a much more sophisticated algorithm that takes into account temperature and non-linear variations also. So it doesn't need to know the exposure time of the master dark.

Is that correct or am I missing something?

Thanks
Peter.
 
Hi Peter,

Sorry for the late answer, I didn't see your post before.

So it doesn't need to know the exposure time of the master dark.

That's correct. Our dark frame optimization algorithm is purely numeric. It does not use any metadata or acquisition properties stored in the images.
 
In this case, however, we have implemented a simpler and very efficient method (Jean-Luc Starck and Fionn Murtagh, Astronomical Image and Data Analysis, Springer-Verlag, 2002, pp. 37-38) known as k-sigma noise estimation:

@Juan Conejero We happen to have the 2nd edition of that book, which is from 2006. Any idea if what you implemented from your referenced 1st edition is still the same? For this edition I couldn't get a preview of the table of content to see the equivalent subsection that you are referring at p.37-38 of the 1st edition. I believe that would be somewhere in the "filtering" chapter 2 which talks about wavelet. In my 2nd edition, here's what the table of content look like. Do you see if any of the subsections matches what you were referencing in the 1st edition? (hopefuly it's still in here...)

Thanks

2. Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Multiscale Transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.1 The A Trous Isotropic Wavelet Transform . . . . . . . . . . . 31
2.2.2 Multiscale Transforms Compared
to Other Data Transforms . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2.3 Choice of Multiscale Transform . . . . . . . . . . . . . . . . . . . . 36
2.2.4 The Multiresolution Support . . . . . . . . . . . . . . . . . . . . . . . 37
2.3 Significant Wavelet Coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.3.2 Noise Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.3.3 Automatic Estimation of Gaussian Noise . . . . . . . . . . . . 40
2.3.4 Detection Level Using the FDR . . . . . . . . . . . . . . . . . . . . 48
2.4 Filtering and Wavelet Coefficient Thresholding . . . . . . . . . . . . . 50
2.4.1 Thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.4.2 Iterative Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
2.4.3 Other Wavelet Denoising Methods . . . . . . . . . . . . . . . . . . 52
2.4.4 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.4.5 Iterative Filtering with a Smoothness Constraint . . . . . 56
2.5 Filtering from the Curvelet Transform . . . . . . . . . . . . . . . . . . . . 57
2.5.1 Contrast Enhancement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
2.5.2 Curvelet Denoising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
2.5.3 The Combined Filtering Method . . . . . . . . . . . . . . . . . . . 61
2.6 Haar Wavelet Transform and Poisson Noise . . . . . . . . . . . . . . . . 63
2.6.1 Haar Wavelet Transform . . . . . . . . . . . . . . . . . . . . . . . . . . 63
2.6.2 Poisson Noise and Haar Wavelet Coefficients . . . . . . . . . 64
2.6.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
2.7 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
 
@Juan Conejero Concretely, how do you choose your range of [k0, k1] in the absence of any metadata in the initial bracketing phase that you described in your 1st post:

Find a range [k0,k1] of dark scaling factors that bracket the minimum of the noise evaluation function. This is the initial bracketing phase.

E.g. a dark frame that was exposed at 10 times longer than the raw light may require a different initial bracketing range than a dark frame that was exposed at the same exposure time than the raw light. So given that you do not use that information, does the GSS guarantee that the minimization of the cost function does not get stuck in a local minimum regardless of the initial range (or is I0 convex with respect to k?)? And also for computational efficiency, do you have a euristic to find a reasonable range for [k0, k1]?
 
Hi Raphael,

@Juan Conejero Concretely, how do you choose your range of [k0, k1] in the absence of any metadata in the initial bracketing phase
do you have a euristic to find a reasonable range for [k0, k1]?

We start with the initial range k0=1/2, k1=2. In virtually all reasonable cases this interval contains the minimum of the dark frame optimization function. In the cases where it doesn't, our initial bracketing routine finds a valid range by extending the initial one at the lower or upper bound as necessary. Our implementation is based on a similar routine included in Numerical Recipes in C 2nd Ed. section 10.1.

does the GSS guarantee that the minimization of the cost function does not get stuck in a local minimum regardless of the initial range (or is I0 convex with respect to k?)?

The cost function in this case (noise evaluation) is strictly unimodal over the initial bracketing range, so this risk does not exist.
 
Thanks @Juan Conejero for the clarification.

I'm trying a simple Python implementation for the purpose of testing some automation outside of PI. I had trouble tracking down the bits of code in the Gitlab repo to see how you were implementing the single-layer wavelet transform. Would something like this be a correct way to do it? (the purpose is not efficiency, but understanding of the algorithm).

Python:
from scipy.signal import convolve2d
import numpy as np

def evaluate_noise(image):
   
    kernel_b3_spline = np.array([
    [1/256, 1/64, 3/128, 1/64, 1/256],
    [1/64, 1/16, 3/32, 1/16, 1/64],
    [3/128, 3/32, 9/64, 3/32, 3/128],
    [1/64, 1/16, 3/32, 1/16, 1/64],
    [1/256, 1/64, 3/128, 1/64, 1/256]
    ])
   
    c0 = image
    c1 = convolve2d(c0, kernel_b3_spline, mode='same', boundary='symm')
    w1 = c0 - c1
    # Work with the 1D-unwrapped array
    w1f = w1.flatten()
    sigma = w1f.std()
    rejmask = np.abs(w1f) < 3 * sigma
    w1f = w1f[rejmask]
    sigmac = sigma
    sigma = w1f.std()
    ds = abs(sigmac - sigma)/sigmac * 100

    while ds >= 1:
        rejmask = np.abs(w1f) < 3 * sigma
        w1f = w1f[rejmask]
        sigmac = sigma
        sigma = w1f.std()
        ds = abs(sigmac - sigma)/sigmac * 100
    return sigma # divide by {sigma^e}_j = 0.889 for level j=1 ?

Then I apply the GSS algorithm to minimize that noise for I - master_bias - k*master_dark, and my k-values are within 3% the ones given by PI as displayed in the process console during the ImageCalibration process. I wonder if the difference is due to the difference in how the convolution is performed by the python package (convolve2d) or if that's just inescapable due to precision/roundoff errors?
Also, in the k-sigma clipping subsection of the Automatic Estimation of Gaussian Noise in the book referenced above, Starck & Murtagh state:
The noise standard deviation at scale j of the data is equal to the noise standard deviation multiplied by the noise standard deviation at scalej of the simulated data.
For j=1, the simulated one is equal to 0.889 for the bi-dimensional case but I could not see if you were applying it in the code (Gitlab), are you applying this correction throughout the minization or did you leave it out? You mentionned you use unscaled noise for the purpose of minimizing the cost function and I was not clear if you referred to this specifically or something else.
Regarding the noise value displayed in the process console during the ImageCalibration process, is this correction finally applied? In fact, I am using that output to compare them with my noise values that I get from mycode above.

Thanks

Raphael
 
Last edited:
  • Like
Reactions: dld
Hi Raphael,

Although I don't use Python, your code looks basically correct to me.

Then I apply the GSS algorithm to minimize that noise for I - master_bias - k*master_dark, and my k-values are within 3% the ones given by PI as displayed in the process console during the ImageCalibration process. I wonder if the difference is due to the difference in how the convolution is performed by the python package (convolve2d) or if that's just inescapable due to precision/roundoff errors?

A 3% difference could be caused by differences between our implementations of the golden section search minimization algorithm. Also note that the ImageCalibration tool optimizes for a central region of 1024x1024 pixels by default, instead of the whole image. At any rate, this shouldn't be caused by roundoff or truncation errors.

are you applying this correction throughout the minization or did you leave it out?

I am not applying this correction because it is not necessary to minimize the cost function for optimization. To find the value of the independent variable that minimizes (or maximizes) a function, multiplication of function values by a constant scalar > 0 is a no-op.

Regarding the noise value displayed in the process console during the ImageCalibration process, is this correction finally applied? In fact, I am using that output to compare them with my noise values that I get from mycode above.

The noise estimates computed by the ImageCalibration process, which are stored as FITS keywords and XISF properties in the calibrated images, have nothing to do with noise estimates used for dark frame optimization. In the former case the multiresolution support algorithm is used by default and of course the estimates include wavelet Gaussian noise scaling factors; in the latter case the k-sigma method is used for the sake of efficiency, and also because accuracy of noise estimates for dark frame optimization is not very critical.

On a side note, we calculated Gaussian noise scaling factors for several wavelet scaling functions and found values that are slightly different from the ones included in Starck's books. In particular, for the B3 spline function we found 0.8907, 0.2007, 0.0856, 0.0413, 0.0205, 0.0103 and 0.0052, respectively for the first seven wavelet layers (à trous algorithm). We ran simulations with very large samples of synthetic Gaussian noise and high-quality pseudo-random number generators.
 
  • Like
Reactions: dld
Back
Top