BUGS · REFERENCES · SEE_ALSO

## NAME

mem -- Perform deconvolution of an image by MEM

## USAGE

`mem input psf model output noise adu `

## DESCRIPTION

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.

*psf*

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

*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).
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.

*guess*

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[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 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.

*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 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.

*damping*

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*
.

*tol*

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.

*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.

*|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.

*inner_iter*
:

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[4]*
is displayed in the parentheses.

*|P>> Error*
:

The relative error in the case of Poisson noise only (similar to chi-sq).

*Normalized*
:

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.

*|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 the 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.

*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[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 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).

## PARAMETERS

- input
- Input degraded image to be deconvolved.

- psf
- Input data file of the point spread function.

- model
- Input model (prior) image used initially in the iteration.

- output
- Output restored image.

- noise
- Readout noise in electrons.

- adu
- A/D conversion constant, electrons/DN.

- poisson=no
- Set to yes if the input image has Poisson noise only.

- tp=0.0
- The total power (flux) of the image. Enter 0 for automatic calculation.

- nsub=1
- The factor for subpixelization (subsampling).

- blksum=yes
- To get a normal pixel, block sum nsub**2 subpixels (or simply resample every nsub'th subpixel in each dimension)?

- guess=""
- Input first guess image used initially in the iteration.

- icf=""
- Input data file of the Intrinsic Correlation Function (ICF).

- sigma=0.0,0.0
- 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).

- fwhm=0.0,0.0
- FWHMs of the Gaussian function (ICF).

- hidden=no
- 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.

- aim=1.0
- A factor for setting the actual target chi-square (poisson=no) or error (poisson=yes) value, namely, target value = Ndata * aim.

- maxiter=100
- Prescribed maximum total number of iterations.

- message=1
- 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.

- m_update=1
- Model updating interval, i.e., the model is updated every m_update'th outer iteration.

- method=0
- 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.

- cj_on=2
- The inner iteration number to turn on conjugate method if method=2.

- opt=3
- 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.

- damping=1.0
- The normalized (1.0 - 100.) damping factor used in the approximate Newton method and conjugate method.

- a_sp=10.0
- Speed factor for renewing the Lagrange multiplier alpha for the data fit constraint.

- a_rate=0.5
- Decrease rate for a_sp, used when a_sp is too large (alpha increases too fast).

- b_sp=3.0
- Speed factor for renewing the Lagrange multiplier beta for the total power constraint.

- tol=0.05,0.05,0.10,0.10
- Convergence tolerances for ME solution, data fit, total power, and the solution of linear equations (method=3).

## EXAMPLES

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

## TIME REQUIREMENTS

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).

## BUGS

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.

## REFERENCES

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.