PixInsight 1.5: Prerelease Information

[Jack Harvey]

Agree

[Nocturnal]

I should add that as long as you're not loosing any data summing reduces noise just as much as averaging. After all averaging is nothing more than summing and then dividing by the number of elements. The division does not change noise it merely re-scales the numbers back to the range of numbers we started with. Since PI works with a range of 0-1 anyway the whole difference between summing and averaging blurs. Both would result in the same data once rescaled to 0-1.

[georg.viehoever]

Hello,

I think some prefer to use median instead on mean because the median is robust against temporary artifact such as planes, flares, insects, ... . Is that correct?

Georg

[vicent_peris]

Hello,

I think some prefer to use median instead on mean because the median is robust against temporary artifact such as planes, flares, insects, ... . Is that correct?

Georg

Yes, but you lose a lot of signal to noise ratio. I think is better to use a rejection algorithm combined with an average.


Vicent.

[Nocturnal]

More sophisticated algorithms look at the stddev for the pixel and then decide on how to average them. If the stddev is low mean/average can be used. If it's high it uses median. Some folks do dithering and median combine to get rid of hot/cold pixels. I think it's better to get rid of them with a hot pixel map and then combining the resulting clean frames with mean or sigma/kappa.

[Juan Conejero]

The key word here is rescaling. Of course, what you say is true: if you add a set of images and then rescale the sum to fit all the data to the available range ([0,1] in the case of PI), then the result is the same (ignoring changes in the redistribution of the data throughout the available dynamic range) as if you divide the sum by the number of images.

Consider a simplified version of the image formation equation:

I = S + N

where I is the image, S represents the (degraded) signal and N is the additive noise term. Assuming that all images in a set are of the same subject, follow the same response curve (for example, all of them are linear), are properly registered, and have the same exposure, we agree that S is approximately constant for all images (actually it isn't constant due to the fact that S is a degraded version of the input image, formed by the convolution of the input image with a point spread function due to atmospheric turbulence, instrumental limitations, etc.). We know also that N, the noise, follows a random distribution and is uncorrelated. So for a set of n images we can write, ignoring signal degradation:

I1 = S + N1
I2 = S + N2
I3 = S + N3
.
.
.
In = S + Nn

Then the mean combination M of this set of n images is:

M = n*S/n +  (N1+N2+N3+...+Nn)/n = S + (N1+N2+N3+...+Nn)/n

The formula above tells that the signal remains invariant, irrespective of the number of images. The noise has a random distribution in each image, so its mean across a set of images tends to cancel out. By averaging a set of n images, the signal-to-noise ratio improves with the square root of n. So by averaging two images we get a SNR improvement of sqrt(2) = 1.41, by averaging four images, the SNR improves by a factor of two, and so on.

With PixInsight, it is very easy to design an experiment that demonstrates noise reduction by averaging (or, equivalently, adding/rescaling) a set of images, using the PixelMath and NoiseGenerator tools.

Here is an example. This is the standard Lena image in 32-bit floating point format:

http://forum-images.pixinsight.com/legacy/noise/01.jpg

Let's duplicate it seven times (to get eight identical images), rename them as A, B, ..., H, and add the same amount of Gaussian noise to all images. In this example I have used the NoiseGenerator tool with the help of ImageContainer:

http://forum-images.pixinsight.com/legacy/noise/02.jpg

Here is the result for one of the eight images:

http://forum-images.pixinsight.com/legacy/noise/03.jpg

Now, instead of calculating the mean of the eight images, we'll just add them. To prevent overflow in the summed result, we can divide each image by eight with PixelMath:

http://forum-images.pixinsight.com/legacy/noise/04.jpg

and here is the result:

http://forum-images.pixinsight.com/legacy/noise/05.jpg

Now let's perform the straight addition of the eight images:

http://forum-images.pixinsight.com/legacy/noise/06.jpg

This is a side-to-side comparison of one of the noisy images (before division) and the summed image:

http://forum-images.pixinsight.com/legacy/noise/07.jpg

Of course, the result would be exactly the same if instead of dividing each image by 8 and adding the eight images, we calculated the arithmetic mean directly, with this PixelMath expression:

(A+B+C+D+E+F+G+H)/8

and rescaling disabled. Or, alternatively, we could just add the eight images (without dividing them) and enable rescaling in PixelMath. The result of the three operations is the same:

- Divide by 8 each image and add the 8 images.
- Add the 8 images and divide the sum by 8.
- Add the 8 images and rescale.

The qualitative comparison above tells us that the SNR improvement is quite high. Now let's evaluate the result quantitatively. The following script implements an iterative noise estimation algorithm due to Jean-Luc Starck:

Code: [Select]

/**
 * Evaluation of Gaussian noise by the iterative k-sigma thresholding method.
 *
 * Reference:
 *    J.L. Starck, F. Murtagh, Astronomical Image and Data Analysis, Springer,
 *    1st ed., 2002, pp. 37-38.
 */

// Working array to gather noise at each iteration.
var N = new Array;

/**
 * Image callback routine to extract noise wavelet coefficients.
 *
 * s     is the value of the current sample, which is actually a wavelet
 *       coefficient in the first wavelet layer
 *
 * ks    k*sigma threshold, with k=3 and sigma is the standard deviation of the
 *       noise in the previous iteration.
 */
function GetNoise( x, y, c, s, ks )
{
   if ( Math.abs( s ) < ks )
      N.push( s );
   else
      this.currentSample = 1.0; // tag this coefficient as significant
}

/**
 * Iterative calculation of the standard deviation of Gaussian noise noise.
 * We use k=3 and return an estimate to within 1% accuracy.
 */
function Noise( img )
{
   // Perform a one-layer wavelet transform using a B3 spline scaling function.
   var W = img.aTrousWaveletTransform( [ 1.0/256, 1.0/64, 3.0/128, 1.0/64, 1.0/256,
                                         1.0/64,  1.0/16, 3.0/32,  1.0/16, 1.0/64,
                                         3.0/128, 3.0/32, 9.0/64,  3.0/32, 3.0/128,
                                         1.0/64,  1.0/16, 3.0/32,  1.0/16, 1.0/64,
                                         1.0/256, 1.0/64, 3.0/128, 1.0/64, 1.0/256 ], 1 );
   for ( var i = 0, s0; ; )
   {
      // Standard deviation of the noise for this iteration.
      var s;
      if ( ++i == 1 ) // first iteration?
         s = s0 = W[0].stdDev();
      else
      {
         // Standard deviation of the set of thresholded coefficients.
         s = Math.stddev( N );

         // Return an estimate to within 1% accuracy.
         var converges = (s0 - s)/s0 < 0.01;
         var enough = i == 20;
         if ( converges || enough )
         {
            console.writeln( (converges ? "Convergence reached" :
                                          "*** Warning: No convergence"),
                             " after ", i, " iterations." );
            return 10*s/0.8908; // back to the image space
         }

         s0 = s;
      }

      // Extract the set of noise pixels for the next iteration.
      N.length = 0;
      W[0].forEachSample( GetNoise, 3*s );
      if ( N.length < 2 )
         return 0;
   }
}

function main()
{
   // Get access to the current active image window.
   var window = ImageWindow.activeWindow;
   if ( window.isNull )
      throw new Error( "No active image" );

   console.show();
   console.writeln( "<b>" + window.currentView.fullId + "</b>" );
   console.writeln( "Calculating noise standard deviation..." );
   console.flush();

   // Obtain an estimate of the noise standard deviation for the current view.
   var s = Noise( window.currentView.image );

   // Output results.
   console.writeln( format( "Noise standard deviation = %.8f, N = %u", s, N.length ) );
}

main();

To use the script, create a new JavaScript file with the Script Editor, copy the code and paste it to the new file, then select the image that you want to evaluate, and run the script (press F9). After a few seconds (or minutes if the image is large), the script will write the estimated standard deviation of the noise on the console.

In my experiment with the Lena image, the script has given the following noise estimates:

A: 0.50855925
B: 0.50687131
C: 0.50699363
D: 0.50706419
E: 0.50645103
F: 0.50811908
G: 0.50709520
H: 0.50594039
sum_of_8: 0.20412443

So the achieved signal-to-noise improvement is 2.5 approximately. Our result is slightly worse than the theoretical improvement, which is the square root of eight, 2.83. I think the main reason is the fact that we are working with discrete images. For example, the addition of noise causes a small but irreversible degradation of the signal in each image.

[Juan Conejero]

Hi Sander,

More sophisticated algorithms look at the stddev for the pixel and then decide on how to average them. If the stddev is low mean/average can be used. If it's high it uses median.

Yes, I know those methods. I'll explain why I haven't implemented one of these techniques into ImageIntegration.

Consider a set of images and a stack of pixels at the same geometrical position p. Our goal is to reject outlier pixels as accurately as possible. To that purpose, it is crucial to obtain a good estimate of the true pixel value at p.

Now calculate the histogram of the values in the stack. Since we are integrating registered images of the same subject, and we have carried out an initial normalization of the images (so that their values are statistically compatible), the histogram has a strong central tendency (a unique and well formed peak). Then we have:

- If the dispersion of values is high, then the median is a robust estimate of the true value of the pixel at p.

- If the dispersion of values is low, then the median and the mean of the stack are nearly equal, and are both equally good estimates of the true value of the pixel at p.

Hence, we conclude that the median is always the best estimate Once we have calculated the median, we can apply an iterative rejection scheme such as k-sigma clipping, min/max, percentile-based clipping, or algorithms based on the properties of CCD detectors. This is what ImageIntegration does.

[vicent_peris]

Hi Sander,

More sophisticated algorithms look at the stddev for the pixel and then decide on how to average them. If the stddev is low mean/average can be used. If it's high it uses median.

Yes, I know those methods. I'll explain why I haven't implemented one of these techniques into ImageIntegration.

Consider a set of images and a stack of pixels at the same geometrical position p. Our goal is to reject outlier pixels as accurately as possible. To that purpose, it is crucial to obtain a good estimate of the true pixel value at p.

Now calculate the histogram of the values in the stack. Since we are integrating registered images of the same subject, and we have carried out an initial normalization of the images (so that their values are statistically compatible), the histogram has a strong central tendency (a unique and well formed peak). Then we have:

- If the dispersion of values is high, then the median is a robust estimate of the true value of the pixel at p.

- If the dispersion of values is low, then the median and the mean of the stack are nearly equal, and are both equally good estimates of the true value of the pixel at p.

Hence, we conclude that the median is always the best estimate Once we have calculated the median, we can apply an iterative rejection scheme such as k-sigma clipping, min/max, percentile-based clipping, or algorithms based on the properties of CCD detectors. This is what ImageIntegration does.

I think Sander means at time of combining the pixel sets. At p, if the StdDev is high, the resulting value will be Med(p). If low, we will do Avg(p).


Vicent.

[Nocturnal]

Indeed, that's what I meant.

[Juan Conejero]

I think Sander means at time of combining the pixel sets. At p, if the StdDev is high, the resulting value will be Med(p). If low, we will do Avg(p).

The same applies, IMO: assuming that the images have been normalized, the median of each pixel stack is the best estimate of the true pixel value in both cases. So this method is equivalent, in practice, to a median combination. Of course, if the images are not compatible (no normalization), then neither the median nor the mean are good estimates.

In my opinion, the best option is a k-sigma clipping or percentile clipping rejection (depending on the number of images) followed by the arithmetic mean combination of the surviving pixels. I have implemented three additional pixel rejection algorithms: min/max, averaged sigma clipping, and CCD-clip. The first two work well for small sets of images, while the third one requires exact parameters of the detector used (gain and readout noise) and unmodified data (no interpolation and no flat fielding), so it is mainly useful to integrate calibration images. As you can see, I have implemented essentially the same algorithms available in IRAF.

[mmirot]

The real problem is normalization. You need a good one to get good rejection.
This is not a challange when you skies are clear with stable transparancy.
That is no clouds and not much gradients.
A typical night in my climate will often have good exposures, bad exposures and in between.

On good night you can just normalize the whole image or it's center.

I often have images where 30, 50, 70% of the area is good.  I would think normalizing whole image or even a single background region rejects even the good parts.  




Max

[vicent_peris]

The real problem is normalization. You need a good one to get good rejection.

The ImageIntegration module uses my normalization algorithms. I have been able to put a sigma of 0.02 on my Lulin image to reject stars, and averaging worked very well on the comet.

This is not a challange when you skies are clear with stable transparancy.
That is no clouds and not much gradients.
A typical night in my climate will often have good exposures, bad exposures and in between.

On good night you can just normalize the whole image or it's center.

I often have images where 30, 50, 70% of the area is good.  I would think normalizing whole image or even a single background region rejects even the good parts.

The problem you refer to is the presence of varying gradients in the image sequence? This is not a problem of good normalization; this is a problem of the images themselves. If you have strong gradients, you must remove them prior to doing a good outlier rejection.

Another option can be to do a first average of the images making small groups. Inside the groups, the gradients will be similar, so outlier rejection will work. At time of averaging the resulting images from each group, you can do an average without outlier rejection.


Vicent.

[mmirot]

Vincent,

You are correct it is a problem of the images themselves.
Actually, I generally don't have to much LP gradients.

I am really talking about variations in transparency.
I toss out the really bad ones.

I often I find that images a portion of the area has good S/N but the rest is poor.
It has some good data the rest should be rejected.
Most normalzation and data rejection methods don't do much for these images.

[Nocturnal]

Well I'm glad to hear that sigma-clip with mean is the preferred method as that's what I've been doing with DSS

[vicent_peris]

Vincent,

You are correct it is a problem of the images themselves.
Actually, I generally don't have to much LP gradients.

I am really talking about variations in transparency.
I toss out the really bad ones.

I often I find that images a portion of the area has good S/N but the rest is poor.
It has some good data the rest should be rejected.
Most normalzation and data rejection methods don't do much for these images.

Perhaps our ImageIntegration will do a better work... but I'm not sure about this... If a portion of an image has poor S/N areas, a low sigma will reject noisy pixels. The key is to have a good normalization to be able to lower sufficiently the sigma parameter.

Are you talking about wide field images? Having los S/N inside an image is very strange, IMO.  :?


Vicent.