Discussion idea - giving StarTools stdev data from stack?

[khyperia]

Some backstory: I looked for quite a while for a good Linux-compatable stacking program, the closest thing I could find was DeepSkyStacker which had severe UI issues running under Wine. So, I made my own - I threw a quick FFT-based phase correlation algorithm at registration, works pretty good for being developed over a few days.

Anyway, something got me thinking when I was implementing it. After the actual stack, I completely throw away standard deviation information and I'm only left with the sigma-kappa mean (or whatever algorithm I decide to use), which then gets fed to StarTools. It seems like a waste to throw away potentially very informative data!

Then again, I don't know that much about how CCDs work, and perhaps standard deviation on a CCD is constant throughout brightness and can be assumed to be constant across the whole image. Still, it might be useful to have that constant number rather than guessing at it.

So that's my idea - possibly allow StarTools to be fed standard deviation data from an external source (normal tiff image or something, I dunno), to allow more accurate error information and processing. Of course, I don't think many stacking programs (like DSS) actually allow you to output this data, so it may be of limited use. Alternatively, one might feed StarTools a bunch of calibrated/registered images (I think DSS does allow you to output these), and it only does the final stack (as opposed to the whole thing - http://startools.org/forum/viewtopic.php?f=8&t=108), keeping all the information it gained.

Discuss?

[khyperia]

Update, I quick slapped out standard deviation information as well as normal on a stack, it's kinda interesting to look at.

At least on my camera, standard deviation is *not* constant with respect to brightness - you can clearly see the structure of the nebula in the standard deviation image!

I've uploaded a zip of the two (16 bit greyscale) to my copy account, https://copy.com/38WmldGvxmpyCN6U under 2014_05_31/M17MeanStdev.zip (this forum has max attachment size of 512 KiB, the images are 5.5 MB, heh)

[admin]

I'm *very* interested in your endeavours here (anything I can do to help - I'd be happy to!). For starters, there is most definitely room for a good stacking program across all platforms!

I think I've said this before (on Reddit), but I really applaud your thinking and getting into the nitty-gritty of things (user Rowland on these forums is similarly inclined).

You're thinking along very similar lines that got me designing the Tracking feature.

How exactly do you calculate the pixel value in the imageStdev.tiff image? What am I looking at here exactly? Accumulative deviation from the sigma/kappa mean for all sub frames?

Some pseudo code would really help!

[admin]

Am I right in concluding that, for the stddev image, we are looking at a plotted measure of 'uncertainty' for the correct pixel value? (which we should be able to establish)
This would be very valuable indeed.

[Amaranthus]

That would be correct Ivo. It may be valuable to also know the standard error of the mean pixel value, SD/sqrt(n), which would indicate how precisely a given pixel's intensity was estimated. For pixels with a large SE, even after sigma rejection (to remove outliers), it could indicate other problems such as a imperfect well, etc. The SD will give an indication of spread from image-to-image, although I'm not sure quite what this would tell us. Maybe faint pixels are likely to have higher SDs, because not every image will capture those rare photons (i.e. frequency of zero detections might be a useful metric). Just musing... (I use statistical modelling a lot in my day job)

[khyperia]

Here's the code that computes the stacking given pixels from each image (list item 1 is 1st image, 2nd is 2nd image, etc.)

Code: Select all

        public static float Stack(List<float> pixels)
        {
            pixels.Sort();

            var median = pixels[pixels.Count / 2];
            float mean, stdev;
            pixels.MeanStdev(out mean, out stdev);

            const float kappa = 2.0f;
            for (var i = 0; i < pixels.Count; i++)
                if (Math.Abs((pixels[i] - median) / stdev) > kappa)
                    pixels[i] = median;

            //commented is imageStdev.tiff, uncommented is imageMean.tiff
            //pixels.MeanStdev(out mean, out stdev);
            //return stdev;
            return pixels.Sum() / pixels.Count;
        }

I *think* that's what median-sigma-kappa is, although I kinda just flailed together a guess and it worked pretty well at removing hot pixels. Ignoring the kappa threshhold, it's basically just taking the mean and standard deviation of each pixel in the stack. The above code just throws out outliers.

So it's not accumulative deviation, it's `sqrt(sum((pixels - mean) ^ 2))`, where `pixels` is the same corresponding location in all sub frames.

However, I kinda scaled the resulting image to 0-65536 exactly when writing to disk, as in outputPixel=pixel/maxPixelInImage*65536, so the two images aren't at the same scale, unfortunately. If you'd like me to re-run without scaling, I can do that, but for now, the sun just set and my autoguider and filters just came in the mail two days ago!

[admin]

That would be correct Ivo. It may be valuable to also know the standard error of the mean pixel value, SD/sqrt(n), which would indicate how precisely a given pixel's intensity was estimated. For pixels with a large SE, even after sigma rejection (to remove outliers), it could indicate other problems such as a imperfect well, etc. The SD will give an indication of spread from image-to-image, although I'm not sure quite what this would tell us. Maybe faint pixels are likely to have higher SDs, because not every image will capture those rare photons (i.e. frequency of zero detections might be a useful metric). Just musing... (I use statistical modelling a lot in my day job)

Exactly right! This is the sort of measure I'd be after for incorporation into the Tracking data. I'm trying to figure out what I'm looking at in the image khyperia uploaded.
I've just now been playing with it, plugging it into some of Tracking module's algorithms.

The trouble is (and I brought this up with khyperia before) I still don't think his/her(?) data is linear.

What this would be very useful for is for detecting and weeding out read noise as a component of the total signal - I've been experimenting with khyperia's data in that light with some interesting results...

[admin]

If you'd like me to re-run without scaling, I can do that, but for now, the sun just set and my autoguider and filters just came in the mail two days ago!

Keep me posted - I'd love to play with this because you're most definitely on to something. The scaling is not a problem. The (suspected) non-linear nature of your data is however. As with all your data to data, doing a gamma correct (1/2.2) made the tiffs somewhat usable for experimentation, but it would've thrown off the real standard deviation results.

(to perform said gamma correct, normalize your dynamic range to 0..1 , then newpixel = pow(normalizedpixel, 2.2), then scale up again by multiplying by 65535 (or whatever your unity value is). If the source is true sRGB, then the conversion is a little more complicated (but not much) but 2.2 should still approximate well.

[khyperia]

Keep me posted - I'd love to play with this because you're most definitely on to something. The scaling is not a problem. The (suspected) non-linear nature of your data is however. As with all your data to data, doing a gamma correct (1/2.2) made the tiffs somewhat usable for experimentation, but it would've thrown off the real standard deviation results.

I've been busy the past few days and haven't gotten the chance to until now. I just uploaded a gamma stretched (pow 2.2) process to the same folder (https://copy.com/38WmldGvxmpyCN6U, 2014_05_31/M17MeanStdev_gammaStretch.zip). I've developed the stacker a bit more (actually re-wrote it entirely from C# to Python, to both make it more portable and the fact that scipy is wonderful), so if you do a direct comparison for some reason or another and find unexpected differences, that's probably why. Also, the new ones are unscaled, I'm directly converting the float64 to 16 bit tiff.

There's actually four images in that zip, image.tiff, imageStdev.tiff, image_clone.tiff, and imageStdev_clone.tiff. See, trouble is, for some reason or another, the image library I'm using for python (Pillow), when it writes out a tiff, StarTools can't open it. It's not compressed (Windows' Properties dialog tells me so, as well as my code telling it not to compress it), and I looked at the hex and the first two chars are II (so it's little-endian, IBM whatever), so I'm a bit confused why it's not working. Running it through ImageMagick (convert image.tiff image_clone.tiff) fixes the problem, but I've included the original issue images to see if you could shed some light on it.

If you'd like, I could throw up the code on github or someplace so you can run whatever dataset you want through it without going through me. It's in a very hackish state right now, though, so there might be issues. (Python + Pillow + numpy + scipy required)

[admin]

That's awesome - thanks for this! WIll have a play as soon as I can free up some time.