Indeed this is very strange. I've made a similar test to yours, and my results are correct.
I've used this script:
/**
* 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 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();
which is a slightly improved version using the B3 spline scaling function and a more accurate value of the noise sigma in the first wavelet layer. However, the original script from the NGC7000 tutorial yields nearly the same result (actually, this algorithm cannot provide more than about 1% accuracy, so only two decimals are really significant).
The 0.8908 constant is the standard deviation of the first wavelet layer for a random Gaussian distribution of zero mean and unit sigma. I calculated this value using high-quality synthetic data and a large number of experiments (it is an improved value with respect to the original one published by Starck in his book, 0.889).
