Package: stsdas.analysis.restore
Help file updated: Apr93
The Lucy-Richardson method can be derived from the maximum likelihood expression for data with a Poisson noise distribution. Thus, it naturally applies to optical imaging data such as HST. The method forces the restored image to be positive, in accord with photon-counting statistics.
The Lucy-Richardson algorithm generates a restored image through an iterative method. The essence of the iteration is as follows: the (n+1)th estimate of the restored image is given by the nth estimate of the restored image multiplied by a correction image. That is,
original data image = image --------------- * reflect(PSF) n+1 n image * PSF n
where the *'s represent convolution operators and reflect(PSF) is the reflection of the PSF, i.e. reflect((PSF)(x,y)) = PSF(-x,-y). When the convolutions are carried out using fast Fourier transforms (FFTs), one can use the fact that FFT(reflect(PSF)) = conj(FFT(PSF)), where conj is the complex conjugate operator.
You can specify a model image to be used if you have prior information about the image. Normally, the model image is simply a constant value image whose value is the mean of the input data.
The restored data will be properly normalized, and the integrated flux in the image is conserved. The total flux for a particular object should also be conserved.
Snyder's modification to the algorithm lets you estimate the additive
noise (e.g., CCD read-out noise) of the detector, which is added to
both the numerator and denominator in the above equation.
Additive noise is specified by the noise parameter (typically about
13 electrons rms for WF/PC), and the conversion factor factor
between electrons and data numbers is specified by the adu parameter
(typically about 7.5 for WF/PC).
For true photon-counting detectors
such as the FOC, a noise value of 0 and gain factor of 1 are
appropriate.
If the data were generated by averaging or summing several images then
the `noise' and `adu' parameters must be modified. For either the sum or
average of N equal exposures the `noise' parameter should be sqrt(N)
times the noise for a single image. For the sum of N equal exposures the
`adu' parameter is the same as for a single image, while for the average
of N equal exposures the `adu' parameter should be set to N times the adu
value for a single exposure.
Note that when the readout noise is non-zero, the flux in the image is no longer conserved exactly. This is reasonable because the best estimate of the total flux in an image with readout noise is not simply the sum of all the data pixels, but rather is a sum with the noisy pixels near zero counts weighted less than the pixels well above the readout noise.
You can specify both the number of iterations to be run and the limiting reduced chi-squared. The program will stop if the chi-squared reaches the specified limit, or if the maximum iteration count is reached. Reaching a chi-squared of 1 is one measure of convergence, but should not be taken too literally. The chi-squared only measures the overall agreement between the restored image and the data, and does not reflect situations in which low-level noise backgrounds or smooth, broad objects are over-fit (e.g., with noise amplification) and high-level sources are not completely restored. Typical HST restorations require at least 40-50 iterations. If the chi-squared is changing significantly with each iteration, convergence has probably not been achieved. Note that the initial iterations will produce a restored image whose resolution is actually worse than the original; in essence, the first iteration yields an image that is the convolution of the raw data and the PSF.
If you specify a high number of iterations, you can cause artifacts to be formed in the restored frame. These artifacts are the result of building noise spikes into apparent sources, or noise amplification. For some images, such as images of star fields, the iteration may be continued for hundreds or even thousands of iterations without introducing objectionable noise amplification in the restored frame. On the other hand, for images with very diffuse, extended emission, such as comets, even 50 iterations can lead to noisy restorations.
The question of when to stop the iteration is especially difficult when the readout noise is included. In that case the convergence of the iteration is somewhat uneven: brighter objects converge more quickly to their ultimate shapes and brightnesses than do faint objects. The convergence slows below a brightness threshold of noise**2/adu counts, or about 25 counts for WFPC images. (The worst case is one in which one has averaged many images together. In that case you can have high signal-to-noise regions of the image that fall well below the threshold and that take many iterations to converge.) You might therefore like to stop the iterations in some parts of the image while continuing to iterate in other regions. This is not currently an option in this task, but you may be able to achieve this result by judiciously choosing the regions of data to restore.
Note that the Lucy-Richardson technique forces all data values in the restored frame to be positive, and performs best if the background data values are close to zero. Achieving good results with this technique requires some preparation of the data frame:
1) Remove any strong DC offset from the frame by fitting a smooth background function or subtracting a constant value. Add this back after restoration if you consider it important for subsequent analysis. Failure to do this can lead to the formation of circular `holes' around strong point sources. The task parameter `background' allows you to specify either a constant background value or a background image that has been subtracted from the input image. These are required in order for the task to compute the Poisson noise statistics properly. In the iteration process, the background is another additive term in both the numerator and denominator of the equation above.
2) Flag or repair any blatantly bad pixels (data drop-outs, cosmic ray hits). Failure to do this can introduce artifacts around the flawed pixels. The input parameter `maskin' specifies a data quality mask, and the associated parameter `goodpixval' identifies whether good pixels are flagged as 1's or 0's. The standard data quality files for WF/PC may be used by setting `goodpixval = 0'.
The task will automatically detect bad negative pixels that exceed the
nominal distribution associated with the given read-out noise.
You can specify a threshold (as multiples of the standard deviation)
for wildly negative pixels using the nsigma parameter.
Flagged pixels are listed as they are found and the resulting mask
file (including pixels already flagged in the input mask) can be
saved by specifying a file name in the maskout parameter.
The restoration can be computed on a data grid which extends beyond the input data. This is important if there are any strong sources near the edge of the image, or if there is a gradient across the image which has not been removed in a background subtraction step. This implementation of the algorithm uses FFTs for efficient computations, but FFTs assume that the data is cyclic. A strong source at the image boundary will introduce artifacts at the opposite side of the image unless the array size is increased somewhat, typically by at least half the size of the PSF. Extending the image also allows for better restorations when a strong source is located just outside the image boundary, but is visible via the edge of its PSF. The parameters `xsizeout' and `ysizeout' allow the user to specify the size of the output image (larger than size of input image). Again, since FFTs are used for computing the convolutions, users are advised to choose image sizes that are powers of two or that factor well. The tasks `factor' and `listprimes' in the `fourier' package can be used to test values for `xsizeout' and `ysizeout'.
If `xsizeout' and `ysizeout' are specified, you can choose how the output image is positioned in the larger array. By default the image will be centered, with an equal buffer on each side. By setting the parameter `center = no', the buffer will all be applied to the right and top, which may be advantageous if you simply plan to trim the output image back to its original size after running the restoration. If `center = no' then the restored image pixels that would appear on the left or bottom of the image wrap around and appear on the right and top. The resulting restored image is the same regardless of the value of `center' except for this circular shift.
If you want to see intermediate results, the parameter `nsave' can be set to retain the results from every `nsave' iterations. The parameter `update' can be set to write out the new solution to the output file every `update' interations. You can then safely abort the task (if, for example, you set the chi-squared limit too low and the program is computing its 1000th iteration!) without losing your output result. It is slightly more time consuming, however, to update the output file.
The PSF will be properly normalized and centered for the computation. If you are doing WF/PC restorations, make sure that your PSF corresponds to the same or a similar position in the field of view. WF/PC PSFs are strongly space-variant, especially near the corners of the CCD chips. A PSF from one edge of the field will not work well for other parts of the field, or on one of the other chips. If you extract the PSF from observed data, make sure it has been corrected for any pixel defects and has any background removed. Any mismatches between the PSF and the actual data will propagate into the restored image as artifacts. If the PSF has any negative data values, these will be truncated at zero prior to use in the restoration iteration.
The iterations may be resumed from where they were left off on a previous run of the program. Use the parameter `model' to specify the name of the previous result. Be sure you specify the name of the original data frame for the `input' image.
goodpixval is set to 0, i.e., a data quality file, the lucy
task will convert all 0's to 1's and all non-zeroes to 0's.
If `yes' (the default), and if `xsizeout' and/or `ysizeout' are larger than the size of the input image, the restored data will be centered in the larger array. This lets you easily see the restoration solution beyond the original image boundaries. If `no', the image is extended to the right and top only. The computations are identical (owing to the cyclical nature of the FFTs---the only difference is in the way the pixels are wrapped in the output image.
restored_5.hhh, restored_10.hhh,
restored_15.hhh, and restored_20.hhh. The final result is only written
to the designated output file. See the parameter `update' if you only want
to update the output image.
2. Perform another 50 iterations on the `wf_lucy50' image, writing output to the image file `wf_lucy100.hhh':re> lucy wf_image.hhh wf_psf.hhh wf_lucy50.hhh 7.5 13 niter=50
re> lucy wf_image.hhh wf_psf.hhh wf_lucy100.hhh 7.5 13 niter=50 model=wf_lucy50.hhh
data[1,*,*]). Image subsets that do not change dimensionality
do correctly update the WCS information, if necessary. This limitation will
be corrected in a subsequent release.