BUGS · REFERENCES · SEE_ALSO
mem -- Perform deconvolution of an image by MEM
mem input psf model output noise adu
This is a task for 1-D or 2-D image restoration by MEM. In search of the maximum point of the objective function, three methods for optimization are available. Data fit and the total power are controlled by the Lagrange multipliers alpha and beta, respectively, to meet the constraints. (See REFERENCES.)
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.
This version D of the program is featured by
1. The preconditioned conjugate method and the accurate Newton method. (See the items \fImethod\fR and \fIcj_on\fR.) 2. The accurate 1-D search in maximization. (See \fIopt\fR.) 3. Poisson noise only case handling. (See \fIpoisson\fR.) 4. Subpixelization technique. (See \fInsub\fR and \fIblksum\fR.) 5. Model updating technique. (See \fIm_update\fR).
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 pixels and sizes
The input degraded image must have normal pixels (without subpixelization). The input PSF, model image, guess image and ICF must have pixels subsampled by a factor of nsub in each dimension (with subpixelization).
Define the input degraded image size as "normal size". Then the sizes of input PSF and other images (if provided by the user) need not be nsub *(normal size) in each dimension. The actual sizes of the arrays holding the images in the program are equal to the maximum of nsub *(normal size) and input PSF size in each dimension. The read-in areas of the model image, guess image and ICF, if provided by the user, will not exceed this maximum.
The output deconvolved image will have the size nsub *(normal size) in each dimension. That is, the output image is subpixelized. Its pixel value datatype will be forced to be real.
As the general guideline, keep the PSF in the smallest size possible. Perform deconvolution only for the degraded image's area of interest. But be aware of the edge effect and the issue of FFT speed.
9 real and 2 (if poisson =no) or 3 (if poisson =yes) complex arrays are needed in the deconvolution procedure. So, for example, to process a 512x512 image when poisson =no, the required core memory is about (1.0MB*13)*nsub **2=13MB*nsub **2. Additional core memory of (1.0MB*2)*nsub **2=2MB*nsub **2 is needed if poisson =yes.
The image sizes need not be a power of two. However, from the point of view of FFT speed, the array size should be a power of two in each dimension if possible. If this is impractical, then it must be avoided to use an array size in any dimension (usually equal to nsub * 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.
The PSF must be space-invariant. It must have been subpixelized by a factor of nsub in each dimension by the user. 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) must have been subpixelized by a factor of nsub in each dimension. It 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 the 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 new alpha and beta 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 method is most accurate (method =1,2), or the accurate solution of the linear equations is easiest (method =3). 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 , adu and poisson
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. If the input image has Poisson noise only, set poisson =yes, in which case noise will be forced to be zero regardless of the user's input.
The defaults noise =13.0, adu =7.5 and poisson =no are for the WF/PC of the Hubble Space Telescope (from the instrument handbook) before the servicing mission. With FOC, use noise =0, adu =1 and poisson =yes. The noise level may be higher if estimated from the degraded image.
Chose the parameters noise and adu carefully for HST images after the servicing mission. They may vary from image to image even for the same camera. Consult the WFPC II Handbook or other documents, and read the image header file.
An image may be resulted from the combination of N images. In the simplest case where these N images are for the same field of view and with the same properties (noise, exposure time etc.), then the following formulas can be used:
Readout noise C/D conversion constant (electrons) (electrons/DN) Individual image noise adu Image from summation sqrt(N)*noise adu Image from average sqrt(N)*noise N*adu
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). If it cannot be found, the program will look for ME_TP in the guess image (if provided by the user). The existence of this keyword indicates that model or guess 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 provided 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 .
nsub and blksum
Set the subpixelization factor nsub to be greater than 1 for using the subpixelization technique. Then each "normal" pixel of the input degraded image will be "split" to nsub*nsub subpixels in the restoration.
If blksum =yes, then the nsub*nsub subpixels are block summed up to get a normal pixel. Accordingly, the intensity of a normal pixel is distributed evenly among nsub*nsub subpixels.
If blksum =no, then the subpixels are simply resampled every nsub 'th one in each dimension to get the normal pixels. Accordingly, the intensity of a normal pixel is assigned to one of the nsub*nsub subpixels with the others having zero values.
Whether using blksum =yes or no depends on the mechanism of image formation. For HST images and PSFs generated by the Tiny Tim package, \blksum =yes should be used because nsub*nsub subpixels are obtained by dividing a normal pixel uniformly.
User may enter a guess image to be used as the first guess image in the iteration, i.e., the initial image with which the iteration starts. It must have been subpixelized by a factor of nsub in each dimension. Without a user input guess image, the model image will be used to start the iteration.
Usually, in the first run of the task, no model or guess image is entered by the user. A program-generated flat image is used as both the model and the first guess image. In the subsequent runs, the ME image from the previous run should be used as model if m_update =1, or as guess if m_update =9999 (meaning no model updating).
icf, sigma, fwhm and hidden
The image formation is modeled like this:
hidden image * ICF * PSF + noise = degraded image,where * denotes 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. In this case the ICF must have been subpixelized by a factor of nsub in each dimension. 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 =0, then fwhm will be read in and sigma will be calculated from fwhm . sigma and fwhm are treated in the same way. By default, sigma- and fwhm- 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.
The criterion for the final convergence of deconvolution is
(current) chi-square or error = Number of data points * \fIaim\fR.The default aim = 1.0 is for using the "standard" critical value. Setting aim < 1.0 to allow the iteration to go further for "over fitting" the data.
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 21 iterations are allowed.
The intermediate ME image (with hidden =yes or sigma and fwhm equal to zero) may be used as the model (if m_update =1) or guess (if m_update =9999) 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.
*** Normally use the defaults for the following parameters ***
method , cj_on and opt
In search of the maximum point of the objective function, a full step(=1) is first taken in the direction determined by one of the following three methods.
\fImethod\fR=1: The approximate Newton method. \fImethod\fR=2: The preconditioned conjugate method, i.e., the conjugate method based on the approximate Newton method. \fImethod\fR=3: The accurate Newton method.
If the conjugate method is chosen, then it will be turned on when the inner iteration number for fixed alpha and beta is equal to or greater than cj_on . The best choice is cj_on =2. This means that every time when alpha and/or beta is renewed, the approximate Newton method is used for one iteration. After that the conjugate method will be used all through until the iteration converges for these particular alpha and beta.
Set method =0 for automatic switching:
if \fIpoisson\fR=no, \fImethod\fR=2 (using the conjugate method), if \fIpoisson\fR=yes, \fImethod\fR=3 (using the accurate Newton method).This is the best default.
Usually the step determined by one of the above methods is not optimal. If opt =2, an optimal step is calculated and taken by (approximate) one-dimensional search using quadratic extrapolation and cubic interpolation. If opt =3 the step calculated using cubic interpolation will be used as the initial guess to find the step leading to the accurate maximum point of the objective function.
In the approximate Newton method, and therefore in the preconditioned conjugate method as well, 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.
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.
At the beginning of each outer iteration, if the normalized magnitude of the initial gradient of the objective function is out of the range [0.20, 0.51], then alpha and a_sp are adjusted automatically with expectation of 0.35 for the normalized magnitude.
If alpha is too large (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 :
\fIa_sp\fR = \fIa_sp\fR * \fIa_rate\fR.
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 .
This array contains the convergence tolerances for ME solution, data fit, the total power and the solution of linear equations in the accurate Newton method, respectively. 0.05, 0.05, 0.10, and 0.10 are reasonable default values for them.
Most output messages are self-explanatory or can be understood on the basis of the above description. Additional information is given in the following.
|P>> , |A>> and |C>> :
Signify that the messages are related particularly to the case of Poisson noise only, the accurate Newton method and the conjugate method, respectively.
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 or (half and minus) log Poisson likelihood 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.
Totally 5 numbers of inner iterations are printed, the current one being outside while the 4 previous ones being inside the parentheses.
|C>> Conjugacy :
Conjugacy factor in the conjugate method.
|A>> Error :
Relative error in the solution of the linear equations for the accurate Newton method. The tolerance tol is displayed in the parentheses.
|P>> Error :
The relative error in the case of Poisson noise only (similar to chi-sq).
The Normalized current value of chi-sq or error, i.e., value / Ndata. Its target value (i.e., aim ) is shown as the second figure in the parentheses.
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 is displayed in the parentheses.
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.
ME solution found :
The inner iteration has converged for these particular alpha and beta.
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 .
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 keywords written by the task are under the header "mem records:" and have a prefix ME_ . The meanings of these keywords are as follows.
ME_NSUB , ME_BSUM (blksum), ME_NOISE , ME_ADU , ME_TP , ME_SIGM1 , ME_SIGM2 , ME_FWHM1 , ME_FWHM2 : parameters used for deconvolution.
ME_ALPHA , ME_A_SP , and ME_BETA : parameters used in computing this output image. They will be used if this image is input as model/guess to run the task again.
ME_HIDDN : a hidden (or visible) image?
ME_MEIMG : an ME image?
ME_CONVG : a converged image?
ME_NITER : the total number of iterations.
ME_NCONV : the total number of convolutions (including correlations).
ME_MAX , ME_MIN : the maximum and minimum values of the hidden image (not necessarily equal to those of the output image).
- Input degraded image to be deconvolved.
- Input data file of the point spread function.
- Input model (prior) image used initially in the iteration.
- Output restored image.
- Readout noise in electrons.
- A/D conversion constant, electrons/DN.
- Set to yes if the input image has Poisson noise only.
- The total power (flux) of the image. Enter 0 for automatic calculation.
- The factor for subpixelization (subsampling).
- To get a normal pixel, block sum nsub**2 subpixels (or simply resample every nsub'th subpixel in each dimension)?
- Input first guess image used initially in the iteration.
- Input data file of the Intrinsic Correlation Function (ICF).
- Sigmas of an elliptic Gaussian function used as ICF if no ICF file is provided. Enter 0 for using another parameter fwhm (full width at half maximum).
- FWHMs of the Gaussian function (ICF).
- The restored hidden image may be convolved with ICF to get a visible image. Set hidden=yes to output the hidden image, or use hidden=no to output the visible image.
- A factor for setting the actual target chi-square (poisson=no) or error (poisson=yes) value, namely, target value = Ndata * aim.
- Prescribed maximum total number of iterations.
- For the verboseness of the output messages when running the task,
1 (least) - 3 (most).
Output an input summary and some parameters in all cases. message=1: detailed for converged inner iterations, brief for non-converged ones. message=2: detailed for all iterations. message=3: very detailed for all iterations.The meanings of messages will be explained later.
- Model updating interval, i.e., the model is updated every m_update'th outer iteration.
- Methods to determine the search direction.
method=0: Automatic, i.e., method=2 if poisson=no, method=3 otherwise. method=1: Approximate Newton method. method=2: Plus conjugate method. (Preconditioned conjugate method.) method=3: Accurate Newton method.
- The inner iteration number to turn on conjugate method if method=2.
- Optimal 1-D search of the maximum point of the objective function.
opt=1: no (step=1.0). opt=2: using quadratic extrapolation and cubic interpolation. opt=3: using quadratic extrapolation and accurate interpolation.
- The normalized (1.0 - 100.) damping factor used in the approximate Newton method and conjugate method.
- Speed factor for renewing the Lagrange multiplier alpha for the data fit constraint.
- Decrease rate for a_sp, used when a_sp is too large (alpha increases too fast).
- Speed factor for renewing the Lagrange multiplier beta for the total power constraint.
- Convergence tolerances for ME solution, data fit, total power, and the solution of linear equations (method=3).
Because quite a few parameters are needed for running the task, some of which are hidden, we advise the user to edit the parameters before running the task. Nevertheless, we give the following examples in which parameters are entered on the command line.
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.
cl>mem r136[1:90,31:120] psf136 "" me136 13 7.52. 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 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).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 b1 psf "" output 1.0e30 1.0e30 aim=1e-4 maxiter=10004. Restore a PC image with a subsampled PSF by a factor of 4 generated by Tiny Tim.
cl>mem pcimage psf_s4 "" output 13 7.5 nsub=45. Restore an FOC image:
cl>mem focimage focpsf "" output 0 1 poisson=yes
The number of convolutions (including correlations) in each iteration is:
for approx. Newton and conjugate methods: 2, for accurate Newton method: indefinite.In any case the total number of convolutions, including some additional ones rather than in the iterations, is displayed as the last message.
There are two FFTs in one convolution. In the case where the image size is a power of two, the other math operations in one iteration is equivalent roughly to 1.0 FFT in CPU time; otherwise the CPU time for the former is negligible. Thus, the CPU time is estimated to be
CPU time for an FFT * (num. of iterations + 2 * num. of convolutions)As an example, for processing a 256x256 image with poisson=no and using the conjugate method, on Sun 4/370 the CPU time for an FFT is about 2 seconds, so 100 iterations (totally about 200 convolutions) 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).
More reasonable criteria for convergence, especially in the case of Poisson noise only, are to be invented and used.
A more efficient method for solving the linear equations in the accurate Newton method is to be used.
The world coordinates should be corrected if nsub>1.
Nailong Wu, SCARS/STSDAS, STScI.
Restoration: Newsletter of STScI's Image Restoration Project, Summer 1993/Number 1.
Wu, Nailong, "Model updating in MEM algorithm", in _The Restoration of HST Images and Spectra II_, pp.58-63, STScI, 1994, R.J. Hanisch and R.L. White (eds.).
Wu, Nailong, "MEM Task for Image Restoration in IRAF", presented at the ADASS IV Conference, Baltimore, Sept. 1994.
Cornwell, T.J. and Evans, K.F., "A simple maximum entropy deconvolution algorithm", Astron. Astrophys., 143, 1985, pp.77-83.
Marheim J. and Rue H., "New Algorithms for Maximum Entropy Image Restoration", CVGIP: Graphical Models and Image Processing, vol.54, No.3, 1992, pp.223-238.