Examining the stddev of a single pixel in a series of images

[Nocturnal]

Hi,

I've noticed in my ASI1600MM-Cool that dark pixel values vary quite substantially from one pixel to the next. Stddev can be about 25% of the average. This in itself is not a problem as long as dark current for each pixel is constant for a given temperature/gain/offset/duration. The easiest way to test would be to calculate the stddev of a single (or a series, even better) of pixels over a bunch of images. The software that came with HAIP had this feature. Can I also do this with PI without coding javascript? Ideally I'd like to calculate the average and stddev of all stddevs of a series of images. That would really tell me how predictable the dark current is.

Any other suggestions on how to measure this? Create a master dark and subtract it from a single dark and see what's left?

[pfile]

not the perfect answer as as far as i know it does not give you the average, but if you load a bunch of images into Blink and hit the bar graph icon, you'll get the statistics for all the loaded images. i guess if you want to focus on just a few pixels you'd have to pre-crop the images down to just the small area that you want to analyze.

rob

[Nocturnal]

I wrote my own utility to gather stats on a large set of files and dump it into a .csv. This was 'necessary' because the HAIP software had a memory leak that prevented it from working on large sets of data. It is a bit odd that PI doesn't do more multi-image analysis with presentable results. Clearly the statistics module together with javascript would make this possible but doing GUI work that way is so tedious :-)

I think I'll enhance my utility to get this functionality. I'll share it when it's done.

[mschuster]

A bit tedious, but something like this might work:

For each dark, use pixelmath to subract the master from the dark and square the result: (d - m)*(d - m)

Use pixelmath to sum all of these up and divide by the number of darks.

Use pixelmath and take the square root of the result.

This is your per pixel standard deviation.

Use Statistics to view mean and standard deviation across all pixels.

Regards,
Mike

As a sanity check, I think the per pixel standard deviation should be at least read noise.

Use a good master to minimize introduction of residual master read noise.

[Juan Conejero]

This is very easy to do with a bit of JavaScript-fu

In less than 120 lines:

Code: [Select]

/*
 * A script to read average pixel sample values from a rectangular region of a
 * set of images.
 *
 * https://pixinsight.com/forum/index.php?topic=10991.0
 */

// The input directory
#define INPUT_DIR    "/home/juan/tmp/bayer-drizzle-test/4/bpp-test/registered/"
// File name suffix of the image files
#define INPUT_EXT    ".xisf"
// The output plain text file in CSV format
#define OUTPUT_FILE  "/tmp/values.csv"
// Coordinates of the input region (left, top, right, bottom).
// The bottom row and right column are excluded.
#define X0           100
#define Y0           100
#define X1           101
#define Y1           101
// Input format hints
#define HINTS        ""

#include <pjsr/ColorSpace.jsh>
#include <pjsr/SampleType.jsh>

function FindFiles( dirPath, fileSuffix )
{
   if ( !dirPath.endsWith( '/' ) )
      dirPath += '/';
   if ( !fileSuffix.startsWith( '.' ) )
      fileSuffix = '.' + fileSuffix;
   let files = [];
   let find = new FileFind;
   if ( find.begin( dirPath + '*' + fileSuffix )  )
      do
         if ( find.isFile )
         {
            let filePath = dirPath + find.name;
            files.push( filePath );
            console.writeln( "<end><cbr>", filePath );
         }
      while ( find.next() );

   return files;
}

function MeanOfImageRegion( filePath, rect, inputHints )
{
   if ( rect == undefined )
      rect = new Rect( 1, 1 );
   if ( inputHints == undefined )
      inputHints = '';
     
   let suffix = File.extractExtension( filePath ).toLowerCase();
   let F = new FileFormat( suffix, true/*toRead*/, false/*toWrite*/ );
   if ( F.isNull )
      throw new Error( "No installed file format can read \'" + suffix + "\' files." );

   let f = new FileFormatInstance( F );
   if ( f.isNull )
      throw new Error( "Unable to instantiate file format: " + F.name );

   let d = f.open( filePath, inputHints );
   if ( d.length < 1 )
      throw new Error( "Unable to open file: " + filePath );

   let image = new Image( 1, 1, 1, ColorSpace_Gray,
                          d[0].bitsPerSample,
                          d[0].ieeefpSampleFormat ? SampleType_Real : SampleType_Integer );
   if ( !f.readImage( image ) )
      throw new Error( "Unable to read file: " + filePath );

   f.close();

   let values = [];
   for ( let i = 0; i < image.numberOfChannels; ++i )
      values.push( image.mean( rect, i, i ) );

   return values;
}

function main()
{
   console.show();

   console.noteln( "<end><cbr><br>=> Searching for " + INPUT_EXT + " files:" )
   console.noteln( INPUT_DIR );
   let files = FindFiles( INPUT_DIR, INPUT_EXT );
   console.noteln( format( "%u file(s) found.", files.length ) );

   console.noteln( "<end><cbr><br>=> Reading mean sample values for region " + format( "{%d,%d,%d,%d}", X0, Y0, X1, Y1 ) );
   let rect = new Rect( X0, Y0, X1, Y1 );
   let data = [];
   for ( let i = 0; i < files.length; ++i )
      data.push( MeanOfImageRegion( files[i], rect, HINTS ) );

   console.noteln( "<end><cbr><br>=> Writing output file:" )
   console.writeln( OUTPUT_FILE );
   let f = File.createFileForWriting( OUTPUT_FILE );
   for ( let i = 0; i < data.length; ++i )
   {
      for ( let j = 0; ; )
      {
         f.outText( format( "%.15lg", data[i][j] ) );
         if ( ++j == data[i].length )
            break;
         f.outText( ',' );
      }
      f.outText( '\n' );
   }
   f.close();

   let values = new Vector( data.length );
   for ( let i = 0; i < data.length; ++i )
      values.at( i, data[i][0] );
   console.writeln( "<end><cbr><br>First channel statistics:" );
   console.writeln( format( "Count ..... %.u",    values.length ) );
   console.writeln( format( "Minimum ... %.6le", values.minComponent() ) );
   console.writeln( format( "Maximum ... %.6le", values.maxComponent() ) );
   console.writeln( format( "Mean ...... %.6le", values.mean() ) );
   console.writeln( format( "Median .... %.6le", values.median() ) );
   console.writeln( format( "StdDev .... %.6le", values.stdDev() ) );
   console.writeln( format( "AvgDev .... %.6le", 1.2533*values.avgDev() ) );
   console.writeln( format( "MAD ....... %.6le", 1.4826*values.MAD() ) );
   console.writeln( format( "Sn ........ %.6le", 1.1926*values.Sn() ) );
}

main();

This script reads mean pixel sample values from a prescribed rectangular region (which can be a single pixel) of a set of images. The images are found automatically on a given directory. The output data are written to a plain text file in good old CSV format: each line corresponds to an image, and each comma-separated value in a line is the mean for an image channel.

You have to modify the INPUT_DIR, INPUT_EXT, OUTPUT_FILE, X0, Y0, X1 and Y1 macros as described at the top of the script to fit your needs.

The script can be easily extended to generate graphical representations, but I'll leave this as an exercise

BTW, beware of the standard deviation because it is a non-robust estimator. You may prefer MAD or something more sophisticated such as Sn. For the same reason, in case you decide to work with rectangular regions larger than a single pixel, you can change the script to compute median values instead of mean values for each image (line 77), or maybe complicate it a bit to compute a trimmed mean, which can be both a robust and efficient estimator (let me know if you are interested in this).

Let me know if this helps.

[Nocturnal]

Thanks Juan for this script! I'll have to play with it a bit because it suffers the same problem as the HAIP software. It runs out of memory :-)

I'll compare stddev and MAD and see what I get.

[Nocturnal]

I tried to force the image object to be deleted (should not be needed but worth a try) but still GC can not keep up and I run out of memory before all 160 files are read. I can see GC knock down memory consumption at one point but it's too slow. I artificially reduced the number of files to scan and have my answers. Bias values vary quite substantially between exposures for individual pixels. I haven't yet done this type of analysis on my old CCD camera but I will.

Thanks for providing this code!

[aworonow]

Just a "remember this": the error of the mean (eom) is what measures the uncertainty of an average value.

Alex

[Nocturnal]

Hi,

I quick look at the stats shows that the pixel values vary quite a bit from image to image:

Count ..... 100
Minimum ... 4.150454e-003
Maximum ... 6.591897e-003
Mean ...... 5.100175e-003
Median .... 5.127031e-003
StdDev .... 5.626460e-004
AvgDev .... 5.538349e-004
MAD ....... 5.429526e-004
Sn ........ 5.823331e-004

MAD and StdDev are close enough that both give the same indication. The absolute number doesn't matter much in this case. I was hoping for an stddev 2 orders of magnitude smaller than the average. Right now we're at less than 1.

[Juan Conejero]

Hi Sander,

Glad to know it's useful. Try inserting this line:

   image.free();

between lines 77 and 79, just before returning from MeanOfImageRegion(). That will improve memory consumption during script execution.

[Nocturnal]

Thanks Juan, will do that. I also changed the rectangle to 1x1 as that was the goal of the exercise. I now get more reasonable numbers although I'm not sure right now why that would be since you were averaging the pixels in the rectangle.

[Nocturnal]

Hi Juan,

Since I'm dealing with raw files from the camera I'd like to keep my pixel values in integers just like the source format. I multiplied the dumped values * 65535 but this does not result in whole numbers. Because the floats are truncated and I'm not sure how to change that I printed the raw mean values where they get put in the data array:

=> Reading mean sample values for region {200,200,200,200}
Reading FITS image: 16-bit integers, 1 channel(s), 4656x3520 pixels: done
0.005203575316703955
Reading FITS image: 16-bit integers, 1 channel(s), 4656x3520 pixels: done
0.005200042422321341
Reading FITS image: 16-bit integers, 1 channel(s), 4656x3520 pixels: done
0.005201757825855853
Reading FITS image: 16-bit integers, 1 channel(s), 4656x3520 pixels: done
0.005202179120371493
Reading FITS image: 16-bit integers, 1 channel(s), 4656x3520 pixels: done
0.005203593148093765

Multiplying those numbers times 65535 also does not result in whole integers, not even close. The mean of a single integer should be an integer, right? I know you are rescaling 16b int FITS into 0-1 float but that should be reversible.

Can you help me understand what is going on?

[Juan Conejero]

Hi Sander,

=> Reading mean sample values for region {200,200,200,200}

This is the problem. You are defining an empty rectangle because the bottom-right corner is always excluded by convention in PCL (also in the JavaScript runtime, which is based on PCL internally). By passing an empty rectangle to Image.mean(), you are computing the mean of the whole image (actually, of the current selection, which is the whole image by default).

To define a single pixel region at coordinates x=200 y=200, define the rectangle X0=200, Y0=200, X1=201, Y1=201. Let me know if it works this way.

[Nocturnal]

Hi,

That is not very intuitive. A box with coordinates 100, 100, 101, 101 has 4 pixels if the edges are included or 0 if it is not. Certainly the new code works much nicer (pixel values closer to what it expected) so if it is wrong I wonder what it is doing.

Anyway, I changed the box size because I was getting fractional pixel values when multiplied with 65535. That is reasonable when you have the mean of 4 pixels.

So instead of using a rectangle, let's just use a single pixel and avoid any confusion. I'm confused about how to use getPixelValue() to get the value of a pixel at a particular coordinates. If you can show me how to get the value at x, y that would be great.

[Juan Conejero]

It is a usual convention in many (most?) graphics engines. See for example the definition of the RECT structure in the Windows GDI API:

https://msdn.microsoft.com/es-es/library/windows/desktop/dd162897(v=vs.85).aspx

The reason for excluding the bottom right corner from rectangles defined by integer coordinates is that in this way we can compute the width and height by simple subtractions:

width = right - left
height = bottom - top

Thanks to this convention, the same operations can be applied to real and integer rectangles. For real coordinates, the bottom right corner is not excluded.

Use this method to read the value of a single pixel:

Number Image.sample( x, y, channel )