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.jpgLet'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.jpgHere is the result for one of the eight images:

http://forum-images.pixinsight.com/legacy/noise/03.jpgNow, 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.jpgand here is the result:

http://forum-images.pixinsight.com/legacy/noise/05.jpgNow let's perform the straight addition of the eight images:

http://forum-images.pixinsight.com/legacy/noise/06.jpgThis 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.jpgOf 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:

`/**`

* 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.