Ask Craig
Combining images: means, medians, and standard deviations
Apr 15 08 08:57
Q: I hear medians are a good way to stack
images as they can remove things like hot pixels,
cosmic rays, or streaks from satellites. Does
Nebulosity support this?
The short answer is no ... but... When combining images we want something that helps reduce the noise. We'd also like something that is tolerant of "outliers". The mean (average) is great at the first part but lousy at the second part. Medians are not so hot at the first part but great on the second part. What we'd like is something that is good at both parts. Nebulosity supports standard-deviation based filtering of your images to let you keep most of the frames and pitch just the bad ones.
OK, so what is it and why is it better? What are these 1.5, 1.75, etc. thresholds I'm being asked about?
If you were to take a perfect image of a target, each pixel would have its "ideal" or "true" value - how much intensity there is from that part of the target. The trouble is, each time we sample the target (aka each image) we get that true value for a pixel but we also get some noise on top of it. We want, of course, the true value. How do we get rid of that noise?
In statistics, we have several typical ways of describing our data. Think right now just about a single pixel (after alignment). So, we have the same spot in the arm of M51 or something. The most common way is the mean (aka average) of all of our samples (aka images, light frames, etc.). It should tell us the central tendency and therefore estimate the truth. The more samples we have, the better the estimate is since we don't have to rely on just one sample (which has truth plus or minus some noise value) or a few samples. With more samples, the noise tends to cancel and we are left with a closer estimate of the truth (the noise, BTW, tends to follow a 1/sqrt(# samples) rule). We can quantify how much noise there is in our samples with a second way of describing our data. The variance (and its square root, the standard deviation) are the typical ways we do this, telling us how much "spread" there is in our samples.
If we assume the data are "normal", about 70% of all samples will lie within one standard deviation (SD) of the mean (that is, 70% are less than one standard deviation above or one standard deviation below the average). About 95% like within 2 SD of the mean. Below, I show the histogram of 5000 normally-distributed random numbers (pretend you had 5000 light frames!). Samples in green lie within 1 SD of the mean. Yellow (and green) lie within 1.5 SD. Orange (and yellow and green) are within 2 SD and red are outside of 2SD. Now, these are all real samples (nothing like an errant satellite) but we could safely pitch those samples in red or orange and still have a good estimate of the mean. We'd not loose too many samples and we'd take out those that are more likely to be outliers. If a sample comes in that is > 2SD, odds are pretty good it's an errant value (like a hot pixel or satellite). Even if it's not, since we don't have 5000 samples - we have far fewer - filtering these out will help keep our estimate of the mean centered where it should be and not skewed by that outlier. Thinking about this diagram will help us a lot in the next step - understanding what happens during stacking. Just remember that with the standard deviation, we know what kind of values we might expect to find and what type of values are really abnormal (e.g., something 5 SD from the mean is very abnormal as there is only a 0.000057% chance this is a real sample and not the result of something else going on).
OK, given that background, here is what happens during the stack. For each (aligned) pixel, we calculate the mean and standard deviation across all of the images in the stack. If your SD threshold is at 1.5, any samples of that pixel that have an intensity beyond 1.5 SD from the mean are removed and a new average, excluding these samples, is calculated. This, BTW, is why hot pixels are often eliminated using SD stacking - those hot pixel values are very abnormal and lie far away from the mean.
With the filter set at 1.75, it takes a more extreme or "outlying" intensity value to be counted as "bad" than at 1.5. At 2.0, it takes even more abnormal a value to be excluded. Thus, more samples go into the final image using a higher threshold (and more things like semi-hot pixels as well). Typically, filtering values at 1.5 or 1.75 will yield the best results.
Standard-deviation based stacking therefore lets in more good samples than a median and takes out more (>0) bad samples than the mean (average). That's what makes it such a nice technique for filtering out bad samples. Note, you're not filtering out whole frames. This math is done on each pixel. So, frame #1 may have a bad value at pixel 109,231 but be great everywhere else. For all other pixels, this frame's data will be used but for pixel 109,231 it won't.
The technique isn't perfect. With a lot of outliers, the estimate of the standard deviation goes way up. So, we have a bigger "spread" that is considered normal and it takes something more aberrant to get filtered out. There are techniques to get around this, of course as well, but that's a topic for another day.
The short answer is no ... but... When combining images we want something that helps reduce the noise. We'd also like something that is tolerant of "outliers". The mean (average) is great at the first part but lousy at the second part. Medians are not so hot at the first part but great on the second part. What we'd like is something that is good at both parts. Nebulosity supports standard-deviation based filtering of your images to let you keep most of the frames and pitch just the bad ones.
OK, so what is it and why is it better? What are these 1.5, 1.75, etc. thresholds I'm being asked about?
If you were to take a perfect image of a target, each pixel would have its "ideal" or "true" value - how much intensity there is from that part of the target. The trouble is, each time we sample the target (aka each image) we get that true value for a pixel but we also get some noise on top of it. We want, of course, the true value. How do we get rid of that noise?
In statistics, we have several typical ways of describing our data. Think right now just about a single pixel (after alignment). So, we have the same spot in the arm of M51 or something. The most common way is the mean (aka average) of all of our samples (aka images, light frames, etc.). It should tell us the central tendency and therefore estimate the truth. The more samples we have, the better the estimate is since we don't have to rely on just one sample (which has truth plus or minus some noise value) or a few samples. With more samples, the noise tends to cancel and we are left with a closer estimate of the truth (the noise, BTW, tends to follow a 1/sqrt(# samples) rule). We can quantify how much noise there is in our samples with a second way of describing our data. The variance (and its square root, the standard deviation) are the typical ways we do this, telling us how much "spread" there is in our samples.
If we assume the data are "normal", about 70% of all samples will lie within one standard deviation (SD) of the mean (that is, 70% are less than one standard deviation above or one standard deviation below the average). About 95% like within 2 SD of the mean. Below, I show the histogram of 5000 normally-distributed random numbers (pretend you had 5000 light frames!). Samples in green lie within 1 SD of the mean. Yellow (and green) lie within 1.5 SD. Orange (and yellow and green) are within 2 SD and red are outside of 2SD. Now, these are all real samples (nothing like an errant satellite) but we could safely pitch those samples in red or orange and still have a good estimate of the mean. We'd not loose too many samples and we'd take out those that are more likely to be outliers. If a sample comes in that is > 2SD, odds are pretty good it's an errant value (like a hot pixel or satellite). Even if it's not, since we don't have 5000 samples - we have far fewer - filtering these out will help keep our estimate of the mean centered where it should be and not skewed by that outlier. Thinking about this diagram will help us a lot in the next step - understanding what happens during stacking. Just remember that with the standard deviation, we know what kind of values we might expect to find and what type of values are really abnormal (e.g., something 5 SD from the mean is very abnormal as there is only a 0.000057% chance this is a real sample and not the result of something else going on).
OK, given that background, here is what happens during the stack. For each (aligned) pixel, we calculate the mean and standard deviation across all of the images in the stack. If your SD threshold is at 1.5, any samples of that pixel that have an intensity beyond 1.5 SD from the mean are removed and a new average, excluding these samples, is calculated. This, BTW, is why hot pixels are often eliminated using SD stacking - those hot pixel values are very abnormal and lie far away from the mean.
With the filter set at 1.75, it takes a more extreme or "outlying" intensity value to be counted as "bad" than at 1.5. At 2.0, it takes even more abnormal a value to be excluded. Thus, more samples go into the final image using a higher threshold (and more things like semi-hot pixels as well). Typically, filtering values at 1.5 or 1.75 will yield the best results.
Standard-deviation based stacking therefore lets in more good samples than a median and takes out more (>0) bad samples than the mean (average). That's what makes it such a nice technique for filtering out bad samples. Note, you're not filtering out whole frames. This math is done on each pixel. So, frame #1 may have a bad value at pixel 109,231 but be great everywhere else. For all other pixels, this frame's data will be used but for pixel 109,231 it won't.
The technique isn't perfect. With a lot of outliers, the estimate of the standard deviation goes way up. So, we have a bigger "spread" that is considered normal and it takes something more aberrant to get filtered out. There are techniques to get around this, of course as well, but that's a topic for another day.
Will DSLR Shutter work with my camera?
Apr 13 08 22:17
Q: I have a WhizBang Foo-Matic Model
XP-Ultra-Mega camera. Will it work with
DSLR
Shutter?
In general, if your camera has a “bulb” port that allows it to be triggered by a simple external device, it should work just fine. DSLR Shutter is really moronically simple. It sends "go" and "don't go" signals to simple parallel port data lines, the RTS/DTR lines of a serial port, or to the ShoeString DSUSB. The former two are very simple, binary signals. The ShoeString is “semi-intelligent" in that I need a software library from Shoestring and need to code up support for it. Doug's libraries are trivial to use and it's still sending very simple commands. It is still just sending "go" and "don't go" signals to this bulb port on the camera.
There is another type of command that can be sent to a camera directly over it's USB port. These "intelligent" commands differ from camera maker to camera maker and within camera makers can differ across models. To send these commands, one needs:
1) SDK (software development kit) from the manufacturer
2) A camera to work on for development
3) A lot more code
The "bulb" port is generic. Heck, you could use DSLR Shutter to turn the lights on and off in your house with a touch of hardware. The simplicity here comes from using this generic style of interface. Were it to bypass this interface and use the USB link to the camera, DSLR Shutter would grow to several times its size just to support the Canon DIGIC II / III cameras (that use the same SDK). So, if you can use this basic kind of trigger signal (i.e., if you have a "bulb" trigger on your camera), you're in luck. If not, you're not in luck and won't be as I have no intention of expanding this right now.
In general, if your camera has a “bulb” port that allows it to be triggered by a simple external device, it should work just fine. DSLR Shutter is really moronically simple. It sends "go" and "don't go" signals to simple parallel port data lines, the RTS/DTR lines of a serial port, or to the ShoeString DSUSB. The former two are very simple, binary signals. The ShoeString is “semi-intelligent" in that I need a software library from Shoestring and need to code up support for it. Doug's libraries are trivial to use and it's still sending very simple commands. It is still just sending "go" and "don't go" signals to this bulb port on the camera.
There is another type of command that can be sent to a camera directly over it's USB port. These "intelligent" commands differ from camera maker to camera maker and within camera makers can differ across models. To send these commands, one needs:
1) SDK (software development kit) from the manufacturer
2) A camera to work on for development
3) A lot more code
The "bulb" port is generic. Heck, you could use DSLR Shutter to turn the lights on and off in your house with a touch of hardware. The simplicity here comes from using this generic style of interface. Were it to bypass this interface and use the USB link to the camera, DSLR Shutter would grow to several times its size just to support the Canon DIGIC II / III cameras (that use the same SDK). So, if you can use this basic kind of trigger signal (i.e., if you have a "bulb" trigger on your camera), you're in luck. If not, you're not in luck and won't be as I have no intention of expanding this right now.
Is my poor guiding the result of flex?
Apr 13 08 20:20
Q: My guiding doesn't seem great and I get
stars that aren't round. How can I tell if this is
flex in my guide setup?
Unless you're using an SBIG camera with two sensors that both get the light from the same OTA, you're bound to have some differential flex between your main imaging camera and your guide camera. Let's say the main OTA is solidly locked onto the mount and the guide scope is attached by rubber bands to the top of the main scope. Since gravity tends to always pull down, as the whole rig rotates, the deflection of the guide scope relative to the main scope changes. With the guide scope atop the main scope, the guide scope will aim a bit too far down. Rotate the whole rig so that the guide scope is now to one side but keep the cameras fixed and the flex in the guide rig makes it aim "left" or "right" on the image (i.e., gravity is now "camera-left" or "camera-right" rather than "camera-down").
All rigs will have some flex. The question is, how much? Is it a big enough factor to hurt your imaging?
Here's a simple way to measure and hopefully rule out flex. Let's assume that PHD (or whatever package you use) isn't loosing the star and that it's guiding somewhere between perfectly and off by, oh, 3-5 arcsec. So, your mount has some error it just can't keep up with, it overshoots, etc (settings being off, mount being crud, etc.). If that is the case, over time, the star should wobble back and forth but on average be in the same place. We overshoot, we undershoot, we oscillate back and forth past the star - whatever. On average, the star is in the right place, but we have line segments instead of points.
Go out some nice, still night (please don't attempt this with 40 MPH gusts...) and shoot, say an hour of anything at something less than say 1 min exposures. We want something short enough that your mount isn't coughing up furballs during the exposure. Be guiding of course during this.
Now, take those images and do two stacks of them:
1) Do an Align and Combine without any aligning (i.e., fixed alignment). Do this for say 1 min worth of shots, 5 min worth of shots, and for the whole shebang. Does the 1 min stack look pretty clean? How much worse is the 5 min? Now, the big question - how much worse is the whole shebang? If the big stack is a lot worse than the 1-5 min stack, you've got flex. Why? The guide scope kept the guide star on target, plus or minus a few arcsec. That error may show up in the 1-5 min (which can show flex too) but if an hour is a lot longer trail, the only explanation is flex (assuming PHD kept a good lock). A 50 arcsec trail there isn't PHD wobbling.
Note, you can use the Measure Distance tool in Nebulosity to see just how many pixels long your trail is. See how wide a star is in a single shot and see how long the trail is, subtract the two (e.g. 117 pixels - 5 pixels = 112 pixels per hour = 1.8 pixels per minute = you'll not be exposing for 10 minutes with clean stars).
2) Do an Align and Combine with Translation (or Translation + Rotation) in Nebulosity. You'll find in your directory an align.txt file with the dx and dy (shifts in x and y) needed to bring that frame into alignment. You can open this up in something like Excel and plot the total distance the stars moved. Ideally, dx and dy would always be 0. If you're having issues, they won't. Use good old Pythagorus to determine the total distance the stars moved: sqrt(dx*dx + dy*dy). If this is a horizontal line on average with some bumps up and down / some noise, you've got no flex. If there is a real drift component, you've got flex.
Now, how bad is it? The easiest way to check is to fit a straight line to your plot. If you're in Excel, you can just have it add a "trend line". Make sure to have it put the equation on the graph if you do this. That equation will tell you just how bad things are. If you're not in Excel and can't automatically add a trend line, print out your plot and just put a ruler on there and draw a line. Your eye is amazingly good at this fit.
The key number that you need is the slope. In Excel, you'll see an equation like "y = 0.4683x - 1.0786" or some such thing. That 0.4683 is what you need. If doing by hand, the slope is the "rise over the run". That is, how much do you move up on the y axis (pixels) over some amount on the x-axis (time). You may find that your line goes up 10 pixels in 15 minutes, making your slope 10/15 or 0.667 pixels per minute.
Here is a sample from two guide rigs I used:
The first slope of 0.4683 pixels per minute means that if I want a max of 1 pixel worth of extra elongation, I can go for about 2.1 minutes in the first setup and 17 minutes in the second. 1 pixel is pretty harsh, so if we double that to 2 pixels (which will still seem pretty round), I'm talking ~4 minutes and over 30 minutes that I can go before differential flex has become an issue for me.
Unless you're using an SBIG camera with two sensors that both get the light from the same OTA, you're bound to have some differential flex between your main imaging camera and your guide camera. Let's say the main OTA is solidly locked onto the mount and the guide scope is attached by rubber bands to the top of the main scope. Since gravity tends to always pull down, as the whole rig rotates, the deflection of the guide scope relative to the main scope changes. With the guide scope atop the main scope, the guide scope will aim a bit too far down. Rotate the whole rig so that the guide scope is now to one side but keep the cameras fixed and the flex in the guide rig makes it aim "left" or "right" on the image (i.e., gravity is now "camera-left" or "camera-right" rather than "camera-down").
All rigs will have some flex. The question is, how much? Is it a big enough factor to hurt your imaging?
Here's a simple way to measure and hopefully rule out flex. Let's assume that PHD (or whatever package you use) isn't loosing the star and that it's guiding somewhere between perfectly and off by, oh, 3-5 arcsec. So, your mount has some error it just can't keep up with, it overshoots, etc (settings being off, mount being crud, etc.). If that is the case, over time, the star should wobble back and forth but on average be in the same place. We overshoot, we undershoot, we oscillate back and forth past the star - whatever. On average, the star is in the right place, but we have line segments instead of points.
Go out some nice, still night (please don't attempt this with 40 MPH gusts...) and shoot, say an hour of anything at something less than say 1 min exposures. We want something short enough that your mount isn't coughing up furballs during the exposure. Be guiding of course during this.
Now, take those images and do two stacks of them:
1) Do an Align and Combine without any aligning (i.e., fixed alignment). Do this for say 1 min worth of shots, 5 min worth of shots, and for the whole shebang. Does the 1 min stack look pretty clean? How much worse is the 5 min? Now, the big question - how much worse is the whole shebang? If the big stack is a lot worse than the 1-5 min stack, you've got flex. Why? The guide scope kept the guide star on target, plus or minus a few arcsec. That error may show up in the 1-5 min (which can show flex too) but if an hour is a lot longer trail, the only explanation is flex (assuming PHD kept a good lock). A 50 arcsec trail there isn't PHD wobbling.
Note, you can use the Measure Distance tool in Nebulosity to see just how many pixels long your trail is. See how wide a star is in a single shot and see how long the trail is, subtract the two (e.g. 117 pixels - 5 pixels = 112 pixels per hour = 1.8 pixels per minute = you'll not be exposing for 10 minutes with clean stars).
2) Do an Align and Combine with Translation (or Translation + Rotation) in Nebulosity. You'll find in your directory an align.txt file with the dx and dy (shifts in x and y) needed to bring that frame into alignment. You can open this up in something like Excel and plot the total distance the stars moved. Ideally, dx and dy would always be 0. If you're having issues, they won't. Use good old Pythagorus to determine the total distance the stars moved: sqrt(dx*dx + dy*dy). If this is a horizontal line on average with some bumps up and down / some noise, you've got no flex. If there is a real drift component, you've got flex.
Now, how bad is it? The easiest way to check is to fit a straight line to your plot. If you're in Excel, you can just have it add a "trend line". Make sure to have it put the equation on the graph if you do this. That equation will tell you just how bad things are. If you're not in Excel and can't automatically add a trend line, print out your plot and just put a ruler on there and draw a line. Your eye is amazingly good at this fit.
The key number that you need is the slope. In Excel, you'll see an equation like "y = 0.4683x - 1.0786" or some such thing. That 0.4683 is what you need. If doing by hand, the slope is the "rise over the run". That is, how much do you move up on the y axis (pixels) over some amount on the x-axis (time). You may find that your line goes up 10 pixels in 15 minutes, making your slope 10/15 or 0.667 pixels per minute.
Here is a sample from two guide rigs I used:
The first slope of 0.4683 pixels per minute means that if I want a max of 1 pixel worth of extra elongation, I can go for about 2.1 minutes in the first setup and 17 minutes in the second. 1 pixel is pretty harsh, so if we double that to 2 pixels (which will still seem pretty round), I'm talking ~4 minutes and over 30 minutes that I can go before differential flex has become an issue for me.