Interesting thread. I'd like to add just a few resources to facilitate some useful testing work.
To demonstrate the accuracy of the noise estimation algorithms we have implemented in PixInsight, let's consider the following script:
#include <pjsr/ImageOp.jsh>
/*
* A script to check accuracy of noise estimates computed with the
* multiresolution support method of Jean-Luc Starck et al. (known as MRS noise
* evaluation in PixInsight).
*
* Thanks to Mike Schuster for pointing out important errors in the initial
* implementation.
*/
/*
* This is the standard deviation of the additive Gaussian noise that will be
* added to the target image, in the normalized [0,1] range. Change it when
* performing several tests. Something around 0.001 can be reasonable for
* simulation of raw deep-sky images.
*/
#define NOISE_SIGMA 0.0005
/*
* Simulated sensor gain in e-/DN
*/
#define GAIN 1.0
/*
* Set to true to apply Poisson distributed noise to the image.
*/
#define APPLY_POISSON_NOISE true
/*
* Set to true to add Gaussian distributed noise to the image.
*/
#define ADD_GAUSSIAN_NOISE true
/*
* Random Gaussian deviates.
*/
function GaussianRandomDeviates()
{
this.V1 = 0;
this.V2 = 0;
this.S = 0;
this.phase = 0;
/*
* Returns a random deviate from a Gaussian distribution with zero mean and
* unit standard deviation.
*/
this.next = function()
{
/*
* Marsaglia polar method.
*/
let X;
if ( this.phase == 0 )
{
do
{
let U1 = Math.random();
let U2 = Math.random();
this.V1 = 2*U1 - 1;
this.V2 = 2*U2 - 1;
this.S = this.V1*this.V1 + this.V2*this.V2;
}
while ( this.S >= 1 || this.S == 0 );
X = this.V1 * Math.sqrt( -2*Math.ln( this.S )/this.S );
}
else
X = this.V2 * Math.sqrt( -2*Math.ln( this.S )/this.S );
this.phase = 1 - this.phase;
return X;
};
}
/*
* Random Poisson deviates.
*/
function PoissonRandomDeviates()
{
/*
* Returns a random Poisson deviate for a given pixel value.
*/
this.next = function( value )
{
if ( value < 30 )
{
/*
* Implementation of the algorithm by Donald E. Knuth, 1969.
*
* Slow (unusable) for large values.
*/
let L = Math.exp( -value );
let k = 0;
let p = 1;
do
{
++k;
p *= Math.random();
}
while ( p >= L );
return k-1;
}
else
{
/*
* Code adapted from 'Random number generation in C++', by John D. Cook:
* https://www.johndcook.com/blog/cpp_random_number_generation/
*
* The algorithm is from "The Computer Generation of Poisson Random
* Variables" by A. C. Atkinson, Journal of the Royal Statistical Society
* Series C (Applied Statistics) Vol. 28, No. 1. (1979)
*
* Slow (unusable) for small values.
*/
let c = 0.767 - 3.36/value;
let beta = Math.PI/Math.sqrt( 3*value );
let alpha = beta*value;
let k = Math.ln( c ) - value - Math.ln( beta );
for ( ;; )
{
let u = Math.random();
let x = (alpha - Math.ln( (1 - u)/u ))/beta;
let n = Math.trunc( Math.floor( x + 0.5 ) );
if ( n >= 0 )
{
let v = Math.random();
let y = alpha - beta*x;
let tmp = 1 + Math.exp( y );
let lhs = y + Math.ln( v/tmp/tmp );
let rhs = k + n*Math.ln( value ) - this.logFactorial( n );
if ( lhs <= rhs )
return n;
}
}
}
};
this.logFactorial = function( n )
{
let x = n + 1;
return (x - 0.5)*Math.ln( x ) - x + 0.5*Math.ln( 2*Math.PI ) + 1/(12*x);
};
}
/**
* Estimation of the standard deviation of the noise, assuming a Gaussian
* noise distribution.
*
* - Use MRS noise evaluation when the algorithm converges for 4 >= J >= 2
*
* - Use k-sigma noise evaluation when either MRS doesn't converge or the
* length of the noise pixels set is below a 1% of the image area.
*
* - Automatically iterate to find the highest layer where noise can be
* successfully evaluated, in the [1,3] range.
*/
function UnscaledNoiseEvaluation( image )
{
let a, n = 4, m = 0.01*image.selectedRect.area;
for ( ;; )
{
a = image.noiseMRS( n );
if ( a[1] >= m )
break;
if ( --n == 1 )
{
console.warningln( "<end><cbr>** Warning: No convergence in MRS noise evaluation routine - using k-sigma noise estimate." );
a = image.noiseKSigma();
break;
}
}
this.sigma = a[0]; // estimated stddev of Gaussian noise
this.count = a[1]; // number of pixels in the noisy pixels set
this.layers = n; // number of layers used for noise evaluation
}
function main()
{
let window = ImageWindow.activeWindow;
if ( window.isNull )
throw new Error( "No active image" );
console.show();
console.writeln( "<end><cbr><br><b>" + window.currentView.fullId + "</b>" );
console.writeln( "Generating synthetic Gaussian+Poisson noise..." );
console.flush();
let G = new GaussianRandomDeviates;
let P = new PoissonRandomDeviates;
let view = window.currentView;
let image = view.image;
view.beginProcess();
for ( let c = 0; c < image.numberOfChannels; ++c )
{
image.selectedChannel = c;
image.initSampleIterator();
do
{
let v = image.sampleValue();
// Apply Poisson noise to the current pixel value.
if ( APPLY_POISSON_NOISE )
v = P.next( 65535*GAIN*v )/65535/GAIN;
// Add Gaussian white noise.
if ( ADD_GAUSSIAN_NOISE )
v += NOISE_SIGMA*G.next();
image.setSampleValue( Math.range( v, 0, 1 ) );
}
while ( image.nextSample() );
}
view.endProcess();
console.writeln( "<end><cbr><br>Unscaled Gaussian noise standard deviation estimates:" );
console.writeln( "<br>Ch | noise | count(%) | layers |" );
console.writeln( "---+-----------+-----------+--------+" );
for ( let c = 0; c < image.numberOfChannels; ++c )
{
console.flush();
image.selectedChannel = c;
let E = new UnscaledNoiseEvaluation( image );
console.writeln( format( "%2d | <b>%.3e</b> | %6.2f | %d |", c, E.sigma, 100*E.count/image.selectedRect.area, E.layers ) );
console.flush();
}
console.writeln( "---+-----------+-----------+--------+" );
}
main();
The script adds a mix of synthetic Gaussian and Poisson distributed noise to the active image. The NOISE_SIGMA macro defines the standard deviation of the generated Gaussian noise. After adding noise, the script computes noise estimates measured on the image.
If the original image has no significant noise, the computed estimates should be close to NOISE_SIGMA, ideally with errors smaller than a 1% approximately. The Poisson component degrades noise estimation.
The easiest way to use this script is to apply it to a noise-free image. By noise-free I mean an image where the standard deviation of the noise is in the range of a few DN units, say about 10
-5 in the normalized [0,1] range. For your convenience,
here is one such image that simulates a typical raw deep-sky image.
You can use the script also on real raw frames. In such case you should get final noise estimates close to a quadrature sum:
Sqrt( s12 + s22 )
where s
1 is the standard deviation of the noise before applying the script, and s
2 is the standard deviation of the noise added by the script (the value of NOISE_SIGMA, which is 0.0005 by default).
As for the scaled noise estimation scripts mentioned above, I'll copy them here for your convenience:
For normal (non-mosaiced) frames:/**
* Estimation of the standard deviation of the noise, assuming a Gaussian
* noise distribution.
*
* - Use MRS noise evaluation when the algorithm converges for 4 >= J >= 2
*
* - Use k-sigma noise evaluation when either MRS doesn't converge or the
* length of the noise pixels set is below a 1% of the image area.
*
* - Automatically iterate to find the highest layer where noise can be
* successfully evaluated, in the [1,3] range.
*
* Returned noise estimates are scaled by the Sn robust scale estimator of
* Rousseeuw and Croux.
*/
function ScaledNoiseEvaluation( image )
{
let scale = image.Sn();
if ( 1 + scale == 1 )
throw Error( "Zero or insignificant data." );
let a, n = 4, m = 0.01*image.selectedRect.area;
for ( ;; )
{
a = image.noiseMRS( n );
if ( a[1] >= m )
break;
if ( --n == 1 )
{
console.writeln( "<end><cbr>** Warning: No convergence in MRS noise evaluation routine - using k-sigma noise estimate." );
a = image.noiseKSigma();
break;
}
}
this.sigma = a[0]/scale; // estimated scaled stddev of Gaussian noise
this.count = a[1]; // number of pixels in the noisy pixels set
this.layers = n; // number of layers used for noise evaluation
}
function main()
{
let window = ImageWindow.activeWindow;
if ( window.isNull )
throw new Error( "No active image" );
console.show();
console.writeln( "<end><cbr><br><b>" + window.currentView.fullId + "</b>" );
console.writeln( "Calculating scaled noise standard deviation..." );
console.flush();
console.abortEnabled = true;
let image = window.currentView.image;
console.writeln( "<end><cbr><br>Ch | noise | count(%) | layers |" );
console.writeln( "---+-----------+-----------+--------+" );
for ( let c = 0; c < image.numberOfChannels; ++c )
{
console.flush();
image.selectedChannel = c;
let E = new ScaledNoiseEvaluation( image );
console.writeln( format( "%2d | <b>%.3e</b> | %6.2f | %d |", c, E.sigma, 100*E.count/image.selectedRect.area, E.layers ) );
console.flush();
}
console.writeln( "---+-----------+-----------+--------+" );
}
main();
For Bayer CFA frames:/**
* Estimation of the standard deviation of the noise, assuming a Gaussian
* noise distribution.
*
* - Use MRS noise evaluation when the algorithm converges for 4 >= J >= 2
*
* - Use k-sigma noise evaluation when either MRS doesn't converge or the
* length of the noise pixels set is below a 1% of the image area.
*
* - Automatically iterate to find the highest layer where noise can be
* successfully evaluated, in the [1,3] range.
*
* Returned noise estimates are scaled by the Sn robust scale estimator of
* Rousseeuw and Croux.
*/
function ScaledNoiseEvaluation( image )
{
let scale = image.Sn();
if ( 1 + scale == 1 )
throw Error( "Zero or insignificant data." );
let a, n = 4, m = 0.01*image.selectedRect.area;
for ( ;; )
{
a = image.noiseMRS( n );
if ( a[1] >= m )
break;
if ( --n == 1 )
{
console.writeln( "<end><cbr>** Warning: No convergence in MRS noise evaluation routine - using k-sigma noise estimate." );
a = image.noiseKSigma();
break;
}
}
this.sigma = a[0]/scale; // estimated scaled stddev of Gaussian noise
this.count = a[1]; // number of pixels in the noisy pixels set
this.layers = n; // number of layers used for noise evaluation
}
/*
* Returns a Bayer CFA frame converted to a 4-channel image. Individual CFA
* components are written to output channels 0, ..., 3 as follows:
*
* 0 | 1
* --+--
* 2 | 3
*
* Where CFA element #0 corresponds to the top left corner of the input frame.
* The output image will have half the dimensions of the input frame.
*
* If specified, the k0, ..., k3 scalars will multiply their respective output
* channels 0, ..., 3.
*/
function BayerCFAToFourChannel( image, k0, k1, k2, k3 )
{
if ( k0 == undefined )
k0 = 1;
if ( k1 == undefined )
k1 = 1;
if ( k2 == undefined )
k2 = 1;
if ( k3 == undefined )
k3 = 1;
let w = image.width;
let h = image.height;
let w2 = w >> 1;
let h2 = h >> 1;
let rgb = new Image( w2, h2, 4 );
for ( let y = 0, j = 0; y < h; y += 2, ++j )
for ( let x = 0, i = 0; x < w; x += 2, ++i )
{
rgb.setSample( k0*image.sample( x, y ), i, j, 0 );
rgb.setSample( k1*image.sample( x+1, y ), i, j, 1 );
rgb.setSample( k2*image.sample( x, y+1 ), i, j, 2 );
rgb.setSample( k3*image.sample( x+1, y+1 ), i, j, 3 );
}
return rgb;
}
#define K0 1.0
#define K1 1.0
#define K2 1.0
#define K3 1.0
function main()
{
let window = ImageWindow.activeWindow;
if ( window.isNull )
throw new Error( "No active image" );
console.show();
console.writeln( "<end><cbr><br><b>" + window.currentView.fullId + "</b>" );
console.writeln( "Scaled Noise Evaluation Script - Bayer CFA Version." );
console.writeln( "Calculating scaled noise standard deviation..." );
console.flush();
console.abortEnabled = true;
let image = BayerCFAToFourChannel( window.currentView.image, K0, K1, K2, K3 );
console.writeln( "<end><cbr><br>Ch | noise | count(%) | layers |" );
console.writeln( "---+-----------+-----------+--------+" );
for ( let c = 0; c < image.numberOfChannels; ++c )
{
console.flush();
image.selectedChannel = c;
let E = new ScaledNoiseEvaluation( image );
console.writeln( format( "%2d | <b>%.3e</b> | %6.2f | %d |", c, E.sigma, 100*E.count/image.selectedRect.area, E.layers ) );
console.flush();
}
console.writeln( "---+-----------+-----------+--------+" );
}
main();
As noted many times on this forum, scaled noise estimates
must be used for noise-based comparisons performed among different images.
