lucy -- Restore an image using the Lucy-Richardson method algorithm.
lucy input psf output adu noise
This task generates a restored image from an input image and point spread function (PSF) using the algorithm developed independently by Lucy (1974, Astron. J. 79, 745) and Richardson (1972, J. Opt. Soc. Am. 62, 55) and adapted for HST imagery by Snyder (1990, in Restoration of HST Images and Spectra, ST ScI Workshop Proceedings; see also Snyder, Hammoud, & White, JOSA, v. 10, no. 5, May 1993, in press). Additional options developed by Rick White (STScI) are also included.
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 includes effects of 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, 6 electrons for WFPC2), and the conversion factor between electrons and data numbers is specified by the adu parameter (typically about 7.5 for WF/PC, 7 or 14 depending on gain for WFPC2). 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 appropriate 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 WF/PC restorations require at least 40-50 iterations, while WFPC2 restorations may require as few as 15-20 accelerated 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 WF/PC 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 the best 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.
- input [file name]
- Name of the input image of blurred data. Can be either one- or two-dimensional.
- psf [file name]
- Name of the PSF image. Dimensionality must match that of the input image, although the number of pixels per data axis may be larger or smaller than the input image.
- output [file name]
- Name of the output image of restored data. By default, this will be the same size as the input image.
- adu [real]
- Conversion constant from electrons or photons to data numbers (analog to digital unit). Also known as the gain for a CCD. Typical values are 7.5 for WF/PC and 1.0 for FOC.
- noise [real]
- Additive Gaussian noise, or read-out noise, given in electrons. Typical values are 13 for WF/PC and 0 for FOC.
- niter [integer]
- Maximum number of iterations to be performed (optional). The default is 20. Execution will halt before reaching this limit if `convergence' is reached, i.e., (chi-sq <= limchisq).
- limchisq [real]
- Value of chi-squared acceptable for convergence (optional). The default is 1.0. If chi-squared is less than or equal to `limchisq', iterations stop regardless of the value of `niter'.
- model [file name]
- Name of the model image (optional). To resume iterations from the last result, set the model image name to the previous output image name, which will then be used as the first guess. If no model image is specified, the first estimate is an image with all pixel values equal to the image mean value. The model image can also be from another observation, but must be properly registered with the input image.
- background [file name or real value]
- Name of a background image which has been subtracted from the input image, or a number giving the constant offset which has been subtracted from the input image (optional). If a background has been subtracted from the input image it should be specified here in order for the restoration to compute the image statistics properly.
- weight [file name]
- Name of an input weight image which determines the weight of each pixel during the restoration. If the flat_field response varies across the chip, the noise does also. This can be handled in calibrated data where the flat field has been applied by inputting the flat field used in calibration as the weight image. Beware that HST flat field reference files should be inverted before being used in the restoration since weight is a multiplicative image. In cases where long and short exposure images are combined to correct for saturated pixels or bad columns, the relative pixel weight can be adjusted to be proportional to the relative exposure time for each pixel. The weight image is independent of the mask image.
- maskin [file name]
- Name of an input mask image which flags bad pixels in the input image (optional). The mask is assumed to be binary---all 0's or 1's. Make sure that the parameter `goodpixval' is also set appropriately for the mask. Internally the mask is stored as good pixels with value of 1 and bad pixels as 0. When 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.
- goodpixval [integer]
- The flag value corresponding to good pixels in the input mask (optional). Default value is 1 (good pixels flagged as 1, bad pixels flagged as 0).
- nsigma [real]
- The maximum negative deviation of pixel values in the input image, given as a multiple of the expected standard deviation as determined from the `noise' and `adu' parameters (optional). Pixels exceeding this threshhold are flagged as bad data in order to prevent DC offsets from being introduced into the restored image, and the associated holes surrounding strong sources.
- maskout [file name]
- Name of an output mask image which flags any additional bad pixels found in the input image (optional). This mask is a "logical or" of the input mask with any additional pixels masked according to the `nsigma' parameter. Mask files are always output with good pixels flagged as 1 and bad pixels flagged as 0.
- accel_method [none,standard,turbo]
- If set to "standard", (the default), the accelerated iteration algorithm as defined by Hook and Lucy (1992, in Science with the Hubble Space Telescope) is used. If set to "none" the original (unaccelerated) lucy algorithm is used. The value, "turbo" causes the turbo acceleration algorithm as developed by Rick White to be used. The turbo method uses a conjugate gradient approach to modify the search directions from the standard Lucy method in order to speed convergence. The standard accelerated method is always faster than the unaccelerated Lucy method, usually by factors of 3 to 5. The turbo method gains another factor of 2 to 3 for some images, but does not gain much for others. The biggest gains for the turbo method over standard acceleration come for WF/PC images of stellar fields; note that these are also the type of images that require the largest number of iterations.
- xsizeout [integer]
- The size of the output image along the first data axis (optional). The default value of 0 means that the output image will be the same size as the input image. If a non-zero value is specified, it must be larger than the size of the input image.
- ysizeout [integer]
- The size of the output image along the second data axis (optional). The default value of 0 means that the output image will be the same size as the input image. If a non-zero value is specified, it must be larger than the size of the input image.
- center [boolean]
- Center the data in the extended output image? (optional)
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.
- nsave [integer]
- The interval with which to write out intermediate results (optional). The default value of 0 means that no intermediate results are saved. Values of 1 or larger cause the task to write out the current restored image every `nsave' iterations. These output images are named by adding an underscore and the iteration count to the output image root name. For example, if `output = restored.hhh', `nsave = 5', and `niter = 25', the intermediate results would be written as files 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.
- update [integer]
- The interval with which to update the output file (optional). The default value of 0 means that no results are written to the output file until all iterations have been completed, or the chi-squared limit has been reached. Values of 1 or larger cause the task to write out the restored image every `update' iterations, OVERWRITING the previous result. Use this parameter if you want to be sure you do not have to rerun an entire restoration that you decide to abort for some reason. Setting `update = 1' is very cautious, and will cost in execution time. See the parameter `nsave' if you want to preserve intermediate results.
- verbose [boolean]
- If no, (the default), the task will print a statement regarding the total
number of bad pixels detected instead of listing each of the pixels found
(verbose = yes).
Also, in the accelerated mode, additional information about the
determination of the acceleration factor is printed when verbose = yes.
1. Restore the WF/PC image `wf_image' using the PSF `wf_psf' and write the output to an image called `wf_lucy50'. A default model (constant image) is assumed.
re> lucy wf_image.hhh wf_psf.hhh wf_lucy50.hhh 7.5 13 niter=502. 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_lucy100.hhh 7.5 13 niter=50 model=wf_lucy50.hhh
The task does not properly update the world coordinate system information if the input image is a image subset with lesser dimensionality than the parent image (e.g., 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.
The algorithm used in this task was developed independently by Lucy (1974, Astron. J. 79, 745) and Richardson (1972, J. Opt. Soc. Am. 62, 55) and adapted for HST imagery by Snyder (1990, in Restoration of HST Images and Spectra, ST ScI Workshop Proceedings; see also Snyder, Hammoud, & White, JOSA v. 10, no. 5, May 1993, in press). Additional options have been developed by Rick White (ST ScI). Betty Stobie wrote the STSDAS implementation including enhancements. Bob Hanisch authored this help file.