Hi Wade,
Thank you for uploading the images. Here is my solution to this nice image processing problem.
The original data set consists of five images of the 2009/R1 McNaught comet. Each image has been guided on the comet's head and there are time gaps between exposures which allow us to remove the stars by pixel rejection. However, the task is by no means as easy as it may seem from its statement, as we have seen in previous posts on this thread.
Our main goal here is to remove all the stars and keep just the comet's image. This task is not difficult, as we're going to see soon. However, we have an important constraint that poses the true challenge of this problem: we want to preserve, as far as possible, the SNR improvement that results from averaging five images (sqrt(5) = 2.24). Normal integration/rejection methods don't provide a satisfactory result, as we have seen from the original post, so we're going to implement a completely different strategy.
The first problem that we have is that the five images are not
statistically compatible. They have varying background and signal levels, as can be seen on the following screenshot.
Click on the image for a full size versionIn the above screenshot, the five images are being represented with the same screen transfer function (STF) applied. This clearly shows the differences that we have to fix before starting to work on the solution of the problem. Before fixing these differences, any comparison between the five pixels at any given coordinates would be meaningless. This is because each pixel is referred to a different distribution in terms of dispersion and central values. This is what we call
statistically incompatible images. Fixing these incompatibilities is what we call
image normalization.
The ImageIntegration tool implements efficient normalization procedures. In this case, however, the images have strong sky gradients that require a more accurate normalization. The StarAlignment tool implements an accurate mosaic frame adaptation routine that we can use to normalize these images very easily: just register the images with the frame adaptation option enabled, and apply the computed adaptation coefficients with PixelMath. Naturally, the newly registered images will be discarded --we already have the five images registered on the comet's head.
Below is the StarAlignment instance used to compute the adaptation functions.
Click on the image for a full size versionThe frame adaptation transformation is a linear function of the form y = a + b*x, where a is the intercepting Y-axis coordinate and b is the slope dy/dx. These are the computed coefficients for the images in this example:
| Image | a | b |
| 002L | +0.000000 | 1.000000 |
| 003L | +0.001065 | 0.934481 |
| 004L | +0.000950 | 0.915482 |
| 005L | +0.000822 | 0.860648 |
| 006L | +0.002124 | 0.706960 |
These adaptation functions are very easy to apply with PixelMath. The result follows. As before, the STF applied is the same for the five images.
Click on the image for a full size versionThe next step is removing the stars. This can be done very efficiently with a min/max rejection. In this case, we have rejected the four brightest pixels of each stack. Note that this has been equivalent to computing the minimum of each pixel stack.
Click on the image for a full size versionThis achieves our first goal: remove the stars. The result is not perfect, as the brightest stars in the image remain visible, but this is about as good a star removal as we can get. We are not done, however: the signal-to-noise ratio of the resulting image is very poor; approximately as poor as the ratio of one of the individual frames. This happens simply because we have rejected four pixels from each five-pixel stack.
We can improve this result significantly, as we're going to demonstrate. Our next step is integrating the five frames without rejection. This generates the straight average of the five comet images, which includes the comet and
all of the stars.
Click on the image for a full size versionSince we have isolated the comet, we can subtract it from the total integration with PixelMath. This produces a new image with all the stars
without the comet.
Click on the image for a full size versionWe are getting close to our final result. Now we have three images: the comet alone, the stars alone, and everything together. That's all we need to isolate the comet with the best possible SNR improvement. Consider the following equation:
C = I - S if S >= x
C = I - S*(1 - (x - S)/x) if S < x
where C is our final image with the isolated comet, I is the total integration image, S is the image with isolated stars, and x is given by:
x = k*med( S )
where med(f) represents the median value of f and k is a constant.
The above equations subtract the stars from the integrated image. However, they don't subtract the stars equally for all pixels of the integrated image. For relatively bright star pixels (S >= x), the stars are completely removed. For relatively dark star pixels (S < x), the stars are
partially removed as a function of their deviation from a reference value. The reference value is the median of the stars image, and the k constant is an empiric adaptation multiplier that must be found by trial and error. By avoiding a total subtraction, we are preserving most of the signal-to-noise improvement that we have achieved by averaging five frames.
In PixelMath language, the above equations can be written as:
x = k*med( stars );
$T - iif( stars > x, stars, stars*(1 - (x - stars)/x) )with the following symbols defined:
x, k=2.0I've found that k=2 is a good value for the images in this example. Fortunately, the k constant isn't very critical; a good value can be found after three or four iterations. The result after applying this PixelMath transformation can be seen on the following screenshot.
Click on the image for a full size versionFinally, the image has some residual noise that is basically impulsional noise composed of pixel "holes" (pepper noise) as a result of rejected bright pixels in previous stages of our processing. This noise is very easy to fix with a median filter. This is the image before noise reduction:
Click on the image for a full size versionand after noise reduction, my final result:
Click on the image for a full size versionHope this helps!
