Package: "stsdas.analysis.restore"
Help file updated: Nov93
message=1: basic iteration summary only. message=2: plus input summary. message=3: plus additional iteration summary.The meanings of messages will be explained later.
opt=1: no (step=1.0). opt=2: using quadratic inter/extrapolation. opt=3: using quadratic extrapolation and cubic interpolation.
The program uses a double iteration scheme: the values of `alpha` and `beta` are changed in the outer iteration, while the inner iteration is for finding the ME solution for the particular `alpha` and `beta`. The major revision in this version of program is the use of the model updating technique. (See the item `model` and `m_update`, and REFERENCES.) Choosing appropriate parameters is a rather difficult matter in running an MEM program. Every effort has been made for the user's comfort in running the task. However, some difficulties may still be encountered. Some suggestions and comments, and examples are given in the following.
`Image sizes`
The input degraded image, PSF, model image and ICF may have different sizes. They need not be a power of two. The actual sizes of the arrays holding the images in the program are equal to the maximum of the degraded image and PSF sizes. The read-in areas of the model image and ICF, if provided by the user, will not exceed this maximum. The output deconvolved image will have the same size as the input degraded image. Its pixel value datatype will be forced to be real.
As the general guideline, keep the PSF, model image and ICF in the smallest sizes possible. Perform deconvolution only for the degraded image's area of interest. But be aware of the edge effect. 9 real and 2 complex arrays are needed in the deconvolution procedure. So, for example, to process a 512x512 image, the required core memory is somewhat more than 1.0*13=13MB.
From the point of view of FFT speed, the array size should be a power of two if possible. If this is impractical, then it must be avoided to use an array size (usually equal to the degraded image's size) having a large prime number factor. As a good example, on Sun 4/370 a 128x128 FFT takes 0.42 second, while a 127x127 FFT takes 6.9 seconds. (A 512x512 FFT takes 8.7 seconds.)
To assist the user in choosing the array size, three tasks are available: `factor` in the package stsdas.analysis.fourier and `pfactor` in the mem0 package used to factorize a natural number, and `irfftes` in the mem0 package used to determine the FFT speed.
`psf`
The PSF must be space-invariant. In the input file the peak of PSF need not be centered. The peak value or volume of PSF need not be normalized. The program takes care of this.
`model` and `m_update`
The user entered `model` image (e.g., a filtered degraded image) will be used as the model initially in the iteration. If no `model` image is entered, a flat model will be generated by the program for this purpose. In the first run of the task, the model is usually a flat image and simply enter a space. In the subsequent runs, the previous ME solutions should be used as the models.
Once the task has started to run, the model will be updated automatically every `m_update`'th outer iteration. For default `m_update`=1, the model will be updated in every outer iteration, i.e., the ME solution found in the previous outer iteration for the particular alpha and beta will be used as the model to start a new iteration with a new alpha increased from zero. In this way the values of alpha in iterations are minimum, and therefore the approximation in solving for the image change in the Newton-Raphson method is most accurate. Consequently, convergence is fastest. Set `m_update` to a large number, say 9999, for no model updating. The advice is to always use the default value for the fastest convergence and least nonlinearity in photometry.
`noise` and `adu`
The parameter `noise` is the CCD readout noise in electrons, while `adu` is the A/D (analogue-to-digital unit) conversion constant or gain for calculating the Poisson noise variance in DNs.
The default `noise`=13.0 and `adu`=7.5 are for the WF/PC of the Hubble Space Telescope (from the instrument handbook) before the servicing mission. With FOC, `noise`=0 and `adu`=1. (The program may not restore FOC images correctly due to the improper noise model.) The noise level may be higher if estimated from the degraded image.
`tp`
The total power of the image is conserved after deconvolution. A non-zero `tp` provided by the user will be used for the constraint. If its value is zero, the program will look for the keyword ME_TP in the `model` image (if provided by the user). The existence of this keyword indicates that `model` is an ME image from the previous run of the task, and its value will be taken as `tp`. In this way, a constant `tp` will be used in a step-by-step deconvolution of an image. If the user proveded `tp` is zero and no ME_TP is found, the total power of the input degraded image will be taken as `tp`. In optical image deconvolution, normally the user doesn't have to take care about `tp`.
`icf, sigma, fwhm` and `hidden`
The image formation is modeled like this:
hidden image * ICF * PSF + noise = degraded image,where * denotes 2-D convolution. The hidden image, for which the entropy is defined, has no correlations between its pixels. The correlations between pixels in the so-called visible image, equal to hidden image * ICF, are introduced by convolving the hidden image with ICF. This visible image is what we really see.
ICF may be input from a data file. If no file is provided, an elliptic Gaussian function is generated as ICF. Its sigma, or FWHM in each dimension may be entered, `sigma`=`fwhm`/sqrt(8ln2). In the first dimension, for example, if `sigma[1]`=0, then `fwhm[1]` will be read in and `sigma[1]` will be calculated from `fwhm[1]`. `sigma[2]` and `fwhm[2]` are treated in the same way. By default, `sigma[1]-[2]` and `fwhm[1]-[2]` are zero. Then the hidden image is identical to the visible image.
The result from deconvolution is a hidden image, which may be output if `hidden`=yes, or convolved with ICF to get the visible image before output if `hidden`=no. `aim`
The criterion for the final convergence of deconvolution is
(current) chi-square = Number of data points * `aim`.
The default `aim` = 1.0 is for using the "standard" critical value
of chi-square. Setting `aim` < 1.0 to allow the iteration to go
farther for "over fitting" the data.
`maxiter`
The maximum number of iterations is prescribed so that the task may stop running before convergence to output an intermediate ME image. After `maxiter` iterations, if an ME solution for the final `alpha` and `beta` has not been found, maximally extra 20 iterations are allowed. The intermediate ME image (with `hidden`=yes or `sigma` and `fwhm` equal to zero) may be used as the model to run the task again. This procedure of step-by-step deconvolution may be repeated until the deconvolution converges (the data and total power constraints are satisfied).
As a matter of fact, this model updating mechanism is built in the source code (model updating). Therefore, with `m_update`=1 and a large `maxiter`, you need not perform deconvolution manually in a step-by-step manner unless you want the image files of intermediate results.
`damping`
In the approximate Newton-Raphson method, the non-diagonal elements of the Hessian of the objective function are simply ignored. A compensation may be made by increasing the diagonal elements (damping). However, this compensation is not straightforward. By default, `damping` = 1.0, meaning no compensation. Set `damping` = 100.0 for a full compensation. The advice is not to change this default value unless you fully understand the algorithm.
`opt`
In search of the maximum point of the objective function, a full step(=1) is first taken in the direction determined by the Newton-Raphson method. This step may not be optimal. If `opt`=2, an optimal step is calculated and taken by one-dimensional search using quadratic inter/extrapolation. Set `opt`=3 for using quadratic extrapolation and cubic interpolation. In either case, the step is limited to be not greater than 4. The advice is not to change the default unless you fully understand the algorithm.
`a_sp`, `a_rate` and `b_sp`
The Lagrange multiplier `alpha` increases in the iteration, starting with zero for the input or updated model. Its increase speed is controlled by the parameter `a_sp`, which may be set to 5 through 20 for reducing the data misfit at reasonable speed.
The parameter `a_sp` will be increased if the following two conditions are satisfied: (1) the numbers of inner iterations for convergence are not greater than 2 successively for 5 times, and (2) the last increase of `a_sp` took place at least 5 outer iterations before. The rate of increase is controlled by `a_rate`:
`a_sp` = `a_sp` / `a_rate`.
`a_rate`=0.5 by default.
`alpha` and `a_sp` are adjusted automatically so that at the beginning of each outer iteration, the normalized magnitude of the initial gradient of the objective function falls in the range [0.10, 0.51]
If `alpha` increases too fast (`a_sp` is too large), then for a particular `alpha` a large number of iterations may be needed to find an ME solution. (This number, `inner_iter`, may be output if `message`=3.) If this happens, the current `alpha` and `a_sp` will be decreased when mod((current) `inner_iter`,8)=0. The rate of decrease depends on the parameter `a_rate`:
`a_sp` = `a_sp` * `a_rate`.
The parameter `b_sp` plays a similar role as `a_sp`. However, the constraint on the total power is not as crucial as the one on data fit, so `b_sp` is set to a smaller value, being 3 by default. No parameter like `a_rate` is needed for `beta` and `b_sp`.
`tol`
This array contains the convergence tolerances for ME solution, data fit and the total power, respectively. 0.05, 0.05, and 0.10 are reasonable default values for them.
`message`
Most output messages are self-explanatory or can be understood on the basis of the above description. Additional information is given in the following.
`Vol/max of ACF of combination of PSF and ICF`:
This figure is an indication of how far the Hessian of the (half) chi-square deviates from a diagonal matrix, and therefore how difficult the deconvolution would be. A small value of it implies that the deconvolution can be done with ease.
`inner_iter`:
Totally 5 numbers of inner iterations are printed, the current one being outside while the 4 previous ones being inside the parentheses. The minus sign indicates that the parameter `a_sp` was increased recently. `Normalized`:
The Normalized current value of chi-square, i.e., chi-sq / Ndata.
`|gradJ|/|1|`:
The normalized `|gradJ|`, i.e., the ratio between the magnitudes of the gradient of the objective function and unit vector, used to indicate the degree of approximation to the ME solution. (This value is zero for an exact ME solution.) The tolerance `tol[1]` is displayed in parentheses.
`test`:
The value 1.0 - cos<grad H2 - `beta` * grad F, `alpha` * grad E>. This is an indication of the parallelism between the two vectors shown in the above. This value is zero for an exact ME solution.
`An ME image obtained`:
An ME solution has been obtained for the final `alpha` and `beta`, that is, the objective function has been maximized, `|gradJ|/|1|` <= `tol[1]`.
`Congratulations for convergence !!`:
The iteration has converged, i.e., the entropy has been maximized AND the data fit and total power constraints have been satisfied.
Some parameters and statistics about the restored image are written into the output image header file. All the cards written by the task are under the header "mem records:" and have a prefix `ME_`. The meanings of the keywords are as follows.
`ME_NOISE`, `ME_ADU`, `ME_TP`, `ME_SIGM1`, `ME_SIGM2`, `ME_FWHM1`, `ME_FWHM2`: parameters used for deconvolution.
`ME_HIDDN`: a hidden (or visible) image?
`ME_MEIMG`: an ME image?
`ME_CONVG`: a converged image?
`ME_NITER`: the number of iterations.
`ME_MAX`, `ME_MIN`: the maximum and minimum values of the hidden image (not necessarily equal to those of the output image).
1. Perform deconvolution on a section of the WFC image r136 with the PSF psf136, using a flat model image and no ICF, outputting an image named me136. The maximum iteration number is 100 by default.
2. This example illustrates how to enter the noise parameters to perform deconvolution on an image g1 having only zero-mean Gaussian noise with rms value equal to 2.0 (DNs).cl>mem r136[1:90,31:120] psf136 "" me136 13 7.5
3. For a noiseless blurred image b1, deconvolution can be done perfectly in principle and the iteration can go indefinitely. Set aim to a small value and maxiter to a large value for this purpose.cl>mem g1 psf "" output 2.0e30 1.0e30
Poisson noise variance = (pixel value) / 1.e30 ~= 0.0 (DN), and Gaussian noise rms = 2.e30 (electrons) = 2.e30 / 1.e30 = 2.0 (DNs).
cl>mem b1 psf "" output 1.0e30 1.0e30 aim=1e-4 maxiter=1000
As an example, for processing a 256x256 image on Sun 4/370, the CPU time for an FFT is about 2 seconds, so 100 iterations requires approximately 1000 seconds or 17 minutes. On VAX-8800, the required CPU time is roughly 2200 seconds or 37 minutes ( 4.4 seconds for a 256x256 FFT).CPU time for an FFT * (4 ~ 5) * the number of iterations.
Although the model updating technique is used, this simplified version of the Newton-Raphson method may result in slow convergence when the PSF is very broad. A better approximation to the method is needed.
Restoration: Newsletter of STScI's Image Restoration Project, Summer 1993/Number 1.
Nailong Wu, Model updating in MEM algorithm, presented at "The Restoration of HST Images and Spectra II" Workshop at STScI, Baltimore, Nov. 1993.
Cornwell, T.J. and Evans, K.F., Astron. Astrophys., 143, 1985, pp.77-83.