Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.stsci.edu/stsci/meetings/irw/proceedings/snyderd.ps.gz Äàòà èçìåíåíèÿ: Mon May 23 18:45:13 1994 Äàòà èíäåêñèðîâàíèÿ: Sun Dec 23 05:19:57 2007 Êîäèðîâêà: IBM-866 Ïîèñêîâûå ñëîâà: ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï

The Restoration of HST Images and Spectra II
Space Telescope Science Institute, 1994
R. J. Hanisch and R. L. White, eds.
Compensation for ReadíOut Noise in HST Image Restoration
Donald L. Snyder
Washington University, Department of Electrical Engineering, Box 1127, St. Louis, MO
63130
Carl W. Helstrom
University of California, San Diego CA 92037
Aaron D. Lanterman and Mohammad Faisal
Washington University, St. Louis, MO 63130
Richard L. White
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218
Abstract. Data acquired with the charge coupled device camera on the HST are modeled
as an additive PoissoníGaussian mixture, with the Poisson component representing cumuí
lative counts of objectídependent photoelectrons, objectíindependent photoelectrons, bias
electrons and thermoelectrons, and the Gaussian component representing readíout noise.
Two methods are examined for compensating for readíout noise. One method is based
upon approximating the Gaussian readíout noise by a Poisson noise and then using the
expectationímaximization (modiîed RichardsoníLucy) algorithm for Poisson distributed
data to effect the compensation. This method has been used for restoring HST images.
The second method directly uses the expectationímaximization algorithm derived for the
PoissoníGaussian mixture data. This requires the determination of the conditionalímean
estimate of the Poisson component of the mixture, which is accomplished by the evaluation
of a nonlinear function of the data. The second method requires more computation than the
îrst, but modest improvements in the quality of the restorations are realized, particularly for
fainter objects.
1. Introduction
Various methods are being used to compensate for the resolution loss that is present in image data
acquired by the Hubble Space Telescope. The common goal of these methods is to invert the excess
point spreading caused by spherical aberration in the primary mirror so as to restore the images to
their original design resolution. Inverse problems of this type are notoriously illíposed. Not only can
noise ampliîcation cause serious degradation as improved resolution is sought for such problems
but, also, the details of the implementation of any restoration method can have a signiîcant impact
on the quality of the restoration, with seemingly insigniîcant design choices having the potential to
cause large changes in the result.
There are a variety of noise sources present in HST image data acquired with a CCD camera.
The photoíconversion process by which object light is converted into photoelectrons introduces
objectídependent noise characterized statistically as a Poisson random process. Nonideal effects
introduce extraneous electrons that are indistinguishable from objectídependent photoelectrons.
Examples of this excess #noise# include objectíindependent photoelectrons, bias or #fat zero#
electrons, and thermoíelectrons. We use the term background counts for the cumulative effect of
139

140 Snyder et al.
these extraneous electrons. Background counts are Poisson distributed. Readíout noise further
contributes to the degradation of images acquired with the charge coupled device camera aboard
the HST. This noise is characterized as a Gaussian random process.
We use the following mathematical model due to Snyder, Hammoud, and White (1993) for
describing CCD image data acquired by the HST:
r(j) = n obj (j) + n 0 (j) + g(j); j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1 (1)
where r(j) are the data acquired by reading out pixel j of the CCDícamera array, n obj (j) is the
number of objectídependent photoelectrons, n 0 (j) is the number of background electrons, g(j) is
readout noise, and J is the number of pixels in the CCD camera array. The statistical description
of these quantities is as follows. The random variables n obj (j), n 0 (j), and g(j) are statistically
independent of each other and of n obj (k), n 0 (k), and g(k) for j 6= k. Objectídependent counts
fn obj (j); j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1g form a Poisson process with meanívalue function fï obj (j); j =
0; 1; \Delta \Delta \Delta ; J \Gamma 1g, where
ï obj (j) = fi(j)
I \Gamma1
X
i=0
p(jji)Ö(i) ; (2)
where ffi(j); j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1g accounts for nonuniform #atíîeld response, detector efîciency,
and the spatial extent of the detector array,
fp(jji); i = 0; 1; \Delta \Delta \Delta ; I \Gamma 1; j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1g
is the HST point spread function (PSF), fÖ(i); i = 0; 1; \Delta \Delta \Delta ; I \Gamma 1g is the object's intensity
function. The #atíîeld response function fi(\Delta) is assumed to be known, for example, through
a #atíîeld calibration measurement. The PSF p(\Deltaj\Delta) is also assumed to be known, for example,
from a theoretical model of HST optics or through observations of an unresolved star. Backí
ground counts fn 0 (j); j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1g form a Poisson process with meanívalue function
fï 0 (j); j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1g; this meanívalue function is assumed to be known, for example,
through a darkíîeld calibration measurement or from HST images taken near but not including the
object of interest. Readíout noise fg(j); j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1g forms a sequence of independent,
identically distributed, Gaussian random variables with mean m and standard deviation oe; the
constants m and oe are assumed to be known.
One approach for compensating for blurring while recognizing the statistical properties of noise
encountered in charge coupled device cameras is based on the method of maximumílikelihood
estimation, as discussed by Snyder, Hammoud, and White (1993). This leads to the following
iteration for producing a maximumílikelihood estimate of the object's intensity function:
Ó
Ö new (i) = Ó Ö old (i)
1
ï
fi(i)
J \Gamma1
X
j=0
ß fi(j)p(jji)
Ó
ï old (j)
Ö
F
h
r(j); Ó
ï old (j); oe
i
; (3)
where
ï
fi(i) =
J \Gamma1
X
j=0
p(jji)fi(j) ; (4)
Ó
ï old (j) = fi(j)
I \Gamma1
X
i=0
p(jji) Ó Ö old (i) + ï 0 (j) ; (5)
and where
F
\Theta
r(j); ï; oe
\Lambda
= E
h
n obj (j) + n 0 (j)jr; Ó Ö old
i
= E
h
r(j) \Gamma g(j)jr; Ó Ö old
i
(6)

Compensation for ReadíOut Noise 141
is the conditionalímean estimate of n obj (j) + n 0 (j) in terms of the data r and the intensity estimate
Ó Ö old to be updated. Evaluation of the conditional mean (6) yields
F [r; ï; oe] = ï
p(r \Gamma m \Gamma 1 : ï; oe)
p(r \Gamma m : ï; oe)
; (7)
where p(x : ï; oe) is the Poisson‘Gaussian mixture probability density given by
p(x : ï; oe) =
1
X
n=0
1
n!
ï n e \Gammaï 1
p
2‹oe
exp
ß
\Gamma
1
2oe 2 (x \Gamma n) 2
Ö
: (8)
This iteration is used as follows:
STEP 0. Initialization: select Ó Ö old (\Delta)
STEP 1. Compute Ó Ö new (\Delta) using (3)
STEP 2. If done, display Ó
Ö new (\Delta) else Ó
Ö old (\Delta) Ó Ö new (\Delta) and perform STEP 1
Special cases of the iteration (3) have appeared previously in the literature. If fi(j) = 1 and
ï 0 (j) = 0 for all j, then (3) is the same as that given by Llacer and N#u#nez (1990). As oe tends
to zero, the Gaussian exponentials in Eq. (7) become very concentrated about n = r \Gamma m, and
F [r; ï; oe] tends towards r \Gamma m, in which case the iteration (3) is that of Politte and Snyder (1991)
for image recovery from Poissonídistributed data in the presence of nonuniform detector response
and background counts. If, in addition to oe = 0, there hold fi(j) = 1 and ï 0 (j) = 0 for all j,
then (3) becomes the RichardsoníLucy iteration (Richardson 1972, Lucy 1974) and the SheppíVardi
(1982) expectationímaximization algorithm for Poissonídistributed data.
2. Implementation
The computational burden of implementations of the iteration (3) has motivated the use of approxí
imations for the function (7). However, as we have commented, there is a potential for errors
introduced by approximations to be magniîed when realizing solutions to illíposed stochastic
inverseíproblems. For this reason, we have studied two implementations of (3). The îrst, which is
already in use for restoring HST images, relies on an approximation to F [r; ï; oe] obtained by subí
stituting a Poisson variate for the readíout noise g(j) in Eq. (1). The second relies on an evaluation
of F [r; ï; oe] via a saddlepoint approximation.
We describe these two approaches in the remainder of this section. Results we have obtained
using them are presented and compared in the next section.
2.1. Implementation via Poisson Approximation
For this approximation, we follow the discussion given by Snyder, Hammoud, and White (1993,
p. 1017). The readíout noise g(j) from the jth pixel is a Gaussianídistributed random variable
with mean m and standard deviation oe. If the data r(j) are modiîed to form r(j) \Gamma m + oe 2 ,
thereby converting the readíout noise to g(j) \Gamma m+ oe 2 , which is Gaussian distributed with mean
and variance equal to oe 2 , and if oe 2 is large, then according to Feller (1968), g(j) \Gamma m + oe 2 is
approximately equal in distribution to a Poissonídistributed random variable with mean oe 2 . With
this approximation, the CCD data (1) become
r(j) \Gamma m+ oe 2 = n obj (j) + n 0 (j) + n read\Gammaout (j); j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1 (9)
where n read\Gammaout (j) is Poisson distributed with mean oe 2 . Under this approximation, the data
\Phi
r(j) \Gamma m+ oe 2 ; j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1
\Psi form a Poisson process with meanívalue function
n
ï obj (j) + ï 0 (j) + oe 2 ; j = 0; 1; \Delta \Delta \Delta ; J \Gamma 1
o
:

142 Snyder et al.
As there is now no Gaussian readíout noise, the iteration (3) becomes
Ó
Ö new (i) = Ó
Ö old (i)
1
ï
fi(i)
J \Gamma1
X
j=0
''
fi(j)p(jji)
fi(j)
P I \Gamma1
i=0 p(jji) Ó Ö old (i) + ï 0 (j) + oe 2
# h
r(j) \Gamma m+ oe 2
i
; (10)
which is the RichardsoníLucy iteration modiîed to accommodate nonuniformity of #atíîeld
detectoríresponse and background noise. This is equivalent to making the following approximation
to F [r; ï; oe] in Eq. (7) when the standard deviation of the readíout noise is large:
F [r; ï; oe] ‹ ï
ï + oe 2
i
r \Gamma m+ oe 2
j
: (11)
We note also from Eq. (7) that as oe tends to zero, F [r; ï; ffi] tends towards r \Gamma m, suggesting that
the approximation (11) for large values of oe will be a good choice for small values of oe as well.
2.2. Implementation via Saddlepoint Approximation
An alternative to the Poisson approximation is to evaluate the function (7) as needed for use in
Eq. (3). However, this function depends on r, ï, and oe in a complicated way, so some numerical
method of evaluation is necessary. We have used the saddlepoint integration and saddlepoint
approximation methods introduced by Helstrom (1978) for calculating tail probabilities.
The Laplace transform of the PoissoníGaussianímixture density (8) is
L[p(\Delta : ï; oe)](s) =
Z 1
\Gamma1
p(x : ï; oe)e \Gammasx dx = exp[ï(e \Gammas \Gamma 1) +
1
2 oe 2 s 2 ] ; (12)
and, by the inverse Laplace transform,
p(x : ï; oe) =
1
2‹j
Z c+j1
c\Gammaj1
exp[\Phi(s)]ds ; (13)
where
\Phi(s) = xs +
1
2 oe 2 s 2 + ï(e \Gammas \Gamma 1) : (14)
The saddlepoint integration and saddlepoint approximation methods can be used with (13) to
produce accurate evaluations of the function F [r; ï; oe] in (7). The saddlepoint ~
s, which is real, is
the root of
\Phi 0 (~s) =
d\Phi(s)
ds
j s=~s = x + oe 2 ~ s \Gamma ïe \Gamma~s = 0 : (15)
This root can be determined numerically using a Newton iteration,
~
s new = ~
s old \Gamma \Phi 0 (~s old )
\Phi 00 (~s old )
; (16)
where \Phi 00 (~s) = oe 2 + ïe \Gamma~s . We iterate Eq. (16), from starting values described below, until
fi fi \Gamma
~ s new \Gamma ~ s old
\Delta
=~s new
fi fi ! ffl, where we have used ffl = 10 \Gamma6 as the tolerance parameter. The saddlepoint
approximation to Eq. (13) is
p(x : ï; oe) ‹
1
p
2‹\Phi 00 (~s)
e \Phi(~s) : (17)

Compensation for ReadíOut Noise 143
Starting Values Useful starting values for the Newton iteration (16) can be obtained by making
approximations to the exponential in Eq. (15). For example, replacing the exponential by two terms
of its Taylor series,
e \Gamma~s ‹ 1 \Gamma ~ s ; (18)
and then solving x + oe 2 ~ s \Gamma ï(1 \Gamma ~ s) = 0 suggests using
~
s old =
ï \Gamma x
ï + oe 2 (19)
as a starting value. Although this is adequate for some values of x and ï, for low values of x
combined with large values of ï, this resulted in starting values so far from the saddlepoint that
hundreds of Newton iterations were required for convergence (see Table 1). We have developed
two alternative approaches for establishing more efîcient starting values. One of these is obtained
via the Pad#e approximation (Newton, Gould, and Kaiser 1961):
e \Gamma~s ‹
1 \Gamma 1
2 ~
s
1 + 1
2 ~
s
: (20)
Solving this for ~ s in terms of e \Gamma~s and substituting the result into Eq. (15) suggests using
~
s old = \Gamma ln
2
4
\Gamma
x \Gamma 2oe 2 \Gamma ï
\Delta
+
q \Gamma
x \Gamma 2oe 2 \Gamma ï
\Delta 2
+ 4ï
\Gamma
x + 2oe 2
\Delta
2ï
3
5 (21)
as the starting value, where the sign of the radical is selected to make the argument of the logarithm
positive. Although the Pad#e approximation (20) is inaccurate for j~sj ? 2, we have found that the
Newton iterations (16) started with (21) converge, to within a tolerance parameter ffl = 10 \Gamma6 , in no
more than îve iterations for x and ï ranging from 10 \Gamma1 to 10 5 for both oe = 3 and oe = 13 (see
Table 1). For this reason, we have generally used starting values based on this approximation in all
of our studies for a readíout noise level of oe = 13 electrons. Another alternative is to use a small
number of Picard iterations
~
s new = ln
ß ï
x + oe 2 ~
s old
Ö
(22)
to solve Eq. (15). We start these iterations using (19) and only perform a small number of them,
typically three, before converting to Newton iterations. Also, if the argument of the logarithm in
(22) becomes negative during the Picard iterations, we revert to using (19) to start the Newton
iterations. With this approach, no more than îve Newton iterations were required for convergence
to within a tolerance parameter of ffl = 10 \Gamma6 in evaluating either the numerator or denominator of
(7), and often one or two sufîced.
Shown in Table 1 are the starting values and the numbers of Newton iterations required for
convergence to a tolerance parameter of ffl = 10 \Gamma6 for each of the Taylor, Pad#e, and Picard methods
for determining starting values for oe = 13 and a wide range of values of ï and r \Gamma m. While a
variety of starting values and numbers of Newton iterations are seen, all three methods converged
to the same saddlepoints to the degree of numerical precision shown in the table. Both the Pad#e and
Picard methods required no more than 5 Newton iterations, while the Taylor method often required
many more.
Accommodation of Wide DynamicíRange Because of the wide dynamic range encountered in
HST images, with values of ï ranging over several orders of magnitude, it is necessary to modify
Eq. (7) to avoid over#ow and under#ow, writing it as
F [r; ï; oe] = ï exp[ln p(r \Gamma m \Gamma 1 : ï; oe) \Gamma ln p(r \Gamma m : ï; oe)] ; (23)

144 Snyder et al.
Table 1. Comparison of starting values for Newton iterations.
log 10 ï log 10 (r \Gamma m) Taylor: Pad#e: Picard: sadípoint
Starting Newtons Starting Newtons Starting Newtons
0 0 0.000e+00 1 í0.000e+00 1 0.000e+00 1 0.000e+00
0 1 í5.294eí02 2 í5.295eí02 2 í5.294eí02 2 í5.293eí02
0 2 í5.824eí01 2 í5.982eí01 3 í5.824eí01 2 í5.811eí01
0 3 í5.876e+00 5 í6.497e+00 5 í5.876e+00 5 í5.021e+00
0 4 í5.882e+01 54 í9.176e+00 4 í9.043e+00 2 í9.044e+00
0 5 í5.882e+02 581 í1.151e+01 3 í1.149e+01 2 í1.149e+01
1 0 5.028eí02 2 5.036eí02 2 5.396eí01 4 5.035eí02
1 1 0.000e+00 1 í0.000e+00 1 0.000e+00 1 0.000e+00
1 2 í5.028eí01 3 í5.042eí01 3 í5.028eí01 3 í4.947eí01
1 3 í5.531e+00 6 í4.208e+00 5 í3.354e+00 4 í3.647e+00
1 4 í5.581e+01 54 í6.873e+00 3 í6.785e+00 2 í6.786e+00
1 5 í5.586e+02 554 í9.207e+00 3 í9.195e+00 1 í9.195e+00
2 0 3.680eí01 3 3.970eí01 3 8.725eí01 4 3.934eí01
2 1 3.346eí01 3 3.582eí01 3 6.895eí01 4 3.555eí01
2 2 0.000e+00 1 í0.000e+00 1 0.000e+00 1 0.000e+00
2 3 í3.346e+00 6 í2.005e+00 4 í1.886e+00 3 í1.912e+00
2 4 í3.680e+01 37 í4.572e+00 3 í4.525e+00 2 í4.526e+00
2 5 í3.714e+02 369 í6.904e+00 3 í6.896e+00 1 í6.896e+00
3 0 8.546eí01 5 1.523e+00 4 1.661e+00 4 1.422e+00
3 1 8.469eí01 5 1.495e+00 4 1.615e+00 4 1.400e+00
3 2 7.699eí01 4 1.248e+00 3 1.279e+00 4 1.197e+00
3 3 0.000e+00 1 í0.000e+00 1 0.000e+00 1 0.000e+00
3 4 í7.699e+00 10 í2.275e+00 3 í2.264e+00 2 í2.264e+00
3 5 í8.469e+01 85 í4.602e+00 3 í4.597e+00 1 í4.597e+00
4 0 9.833eí01 6 3.448e+00 4 3.096e+00 4 2.985e+00
4 1 9.824eí01 6 3.422e+00 4 3.077e+00 4 2.972e+00
4 2 9.735eí01 6 3.191e+00 4 2.908e+00 3 2.846e+00
4 3 8.850eí01 6 2.070e+00 3 2.013e+00 3 2.010e+00
4 4 0.000e+00 1 í0.000e+00 1 0.000e+00 1 0.000e+00
4 5 í8.850e+00 11 í2.300e+00 2 í2.299e+00 1 í2.299e+00
5 0 9.983eí01 8 5.694e+00 5 4.871e+00 3 4.811e+00
5 1 9.982eí01 8 5.667e+00 5 4.860e+00 3 4.802e+00
5 2 9.973eí01 8 5.437e+00 5 4.754e+00 3 4.714e+00
5 3 9.883eí01 8 4.321e+00 4 4.084e+00 3 4.081e+00
5 4 8.985eí01 6 2.275e+00 3 2.265e+00 2 2.265e+00
5 5 0.000e+00 1 í0.000e+00 1 0.000e+00 1 0.000e+00
where in the saddlepoint approximation
ln p(x : ï; oe) ‹ \Phi(~s) \Gamma 1
2
ln[2‹\Phi 00 (~s)] ; (24)
Since the factors of 2‹ cancel in determining (7), we omit these in our computations.
Accuracy of the Approximations The saddlepoint integration of Eq. (13) is performed by integratí
ing along the vertical line c = ~ s + jy, \Gamma1 ! y ! 1, through the saddlepoint to obtain
p(x : ï; oe) =
1
‹
e \Gammaï
Z 1
0
Re
ae
exp
ß
x(~s + jy) +
1
2 (~s + jy) 2 + ïe \Gamma(~s+jy)
Öoe
dy : (25)
The trapezoidal rule was used, taking the step size as \Deltay =
\Theta
\Phi 00 (~s)
\Lambda \Gamma1=2 and stopping the summation
when the absolute value of the integrand fell below ffi \Deltay times the accumulated sum, with ffi = 10 \Gamma5 .
Smaller stepísizes did not noticeably alter the results, which were also checked by evaluating Eq. (7)
with the summation in Eq. (8); the results agreed to at least six signiîcant îgures. The discrepancies
may be attributable to roundíoff error in the summations in Eq. (8).
Shown in Table 2 are values of F [r : ï; oe] obtained for the methods of saddlepoint integration
(Eq. 25), saddlepoint approximation (Eq. 17), and Poisson approximation (Eq. 11) for oe 2 = 200
and several values of r \Gamma m and ï. One purpose of this table is to compare #exact" values obtained
with saddlepoint integration to values obtained with the saddleípoint and Poisson approximations.

Compensation for ReadíOut Noise 145
Starting values for the Newton iteration were determined by the Picard method. The saddlepoint
approximation is seen to agree closely with the results of saddlepoint integration. If we take the
values of F [r : ï; oe] obtained in Table 2 by saddlepoint integration as #exact," then the sum of
absolute errors for the saddlepoint and Poisson approximations differ by several orders of magnitude,
being 0.0004 and 19.81, respectively.
Table 2. Evaluations of F (\Delta).
oe 2 ï r \Gamma m spíint spíapprox Píapprox
200 20 í55.0 14.122700 14.122702 13.181818
200 20 í25.0 16.236131 16.236133 15.909091
200 20 5.0 18.638897 18.638900 18.636364
200 20 35.0 21.362958 21.362960 21.363636
200 20 65.0 24.441810 24.441813 24.090909
200 20 95.0 27.910075 27.910078 26.818182
200 200 100.0 153.126879 153.126854 150.000000
200 200 140.0 171.028212 171.028184 170.000000
200 200 180.0 190.001082 190.001050 190.000000
200 200 220.0 209.998996 210.998963 210.000000
200 200 260.0 230.972066 230.972031 230.000000
200 200 300.0 252.868626 252.868590 250.000000
200 500 350.0 396.313041 396.313010 392.857143
200 500 400.0 430.003993 430.003964 428.571429
200 500 450.0 464.555108 464.555081 464.285714
200 500 500.0 499.897950 499.897926 500.000000
200 500 550.0 535.970081 535.970060 535.714286
200 500 600.0 572.714680 572.714662 571.428571
200 500 650.0 610.080121 610.080105 607.142857
Shown in Table 3 are values of F [r : ï; oe] obtained for the methods of saddlepoint approxií
mation, using Eqs. (23) and (24), and Poisson approximation (Eq. 11) with oe = 13, corresponding
to the level of readíout noise present with the HSTs CCD camera, and oe = 3, corresponding to
possible values of readíout noise present with improved CCD cameras, for a wide range of values of
r and ï. The purpose of this table is to compare values obtained with the saddlepoint and Poisson
approximations over a dynamic range encountered in HST image data. In this table, j\DeltaF j is the
magnitude of the value of the difference in values of F determined by the saddlepoint and Poisson
approximations.
The saddlepoint and Poisson approximations are seen from these tables to produce practically
identical results for ï ‹ r \Gamma m. Elsewhere, the Poisson result is consistently less than the saddlepoint
result, with the difference worsening as r \Gamma m strays further and further from ï. For r \Gamma m ! ï, the
error tends to level off to some value, whereas for r \Gamma m ? ï, the error rises rapidly and somewhat
proportionately to r \Gamma m. The percentage error can be quite large, especially when r \Gamma m is much
larger than ï. For instance, for ï = 1, oe = 13 and r \Gamma m = 10 5 , j\DeltaF j/(spíapprox)=99.4%.
In performing Newton iterations to evaluate F (r; ï; oe) in Eq. (7), the saddlepoint for the
denominator was evaluated îrst, and the saddlepoint then determined was used as the starting value
for Newton iterations in evaluating the numerator. Except when x was close to ï, the number of
Newton iterations for the numerator was then no larger than for the denominator and was usually
smaller.
2.3. Computer Implementation
Serial Implementation AC subroutine to compute F [r; ï; oe] using the saddlepoint approximation
with a Newton iteration to determine the saddlepoint is given in the Appendix for starting values
obtained with either the Pad#e or Picard methods. With ï = 10 2 and r \Gamma m = 10 3 , which are typical
values, 2 16 calls of these subroutines take approximately 16 seconds on a SUN Sparcstation IPC
and 8 seconds on a SUN Sparcstation IPX. These times can probably be lowered substantially while
maintaining nearly the same accuracy by using a larger value for the tolerance parameter ffl.

146 Snyder et al.
Table 3. More Evaluations of F (\Delta).
log 10 ï log 10 (r \Gamma m) spíapprox Píapprox j\DeltaF j spíapprox Píapprox j\DeltaF j
(oe = 13) (oe = 13) (oe = 13) (oe = 3) (oe = 3) (oe = 3)
0 0 9.97eí01 1.00e+00 2.92eí03 9.56eí01 1.00e+00 4.43eí02
0 1 1.05e+00 1.05e+00 1.66eí03 2.26e+00 1.90e+00 3.60eí01
0 2 1.78e+00 1.58e+00 2.01eí01 6.27e+01 1.09e+01 5.18e+01
0 3 1.51e+02 6.88e+00 1.45e+02 9.38e+02 1.01e+02 8.37e+02
0 4 8.47e+03 5.98e+01 8.41e+03 9.92e+03 1.00e+03 8.92e+03
0 5 9.81e+04 5.89e+02 9.75e+04 9.99e+04 1.00e+04 8.99e+04
1 0 9.48e+00 9.50e+00 1.34eí02 5.72e+00 5.26e+00 4.57eí01
1 1 9.97e+00 1.00e+01 2.63eí02 9.87e+00 1.00e+01 1.26eí01
1 2 1.64e+01 1.50e+01 1.33e+00 8.11e+01 5.74e+01 2.37e+01
1 3 3.84e+02 6.53e+01 3.18e+02 9.59e+02 5.31e+02 4.28e+02
1 4 8.85e+03 5.68e+02 8.29e+03 9.94e+03 5.27e+03 4.67e+03
1 5 9.84e+04 5.60e+03 9.29e+04 9.99e+04 5.26e+04 4.73e+04
2 0 6.74e+01 6.32e+01 4.18e+00 1.69e+01 9.17e+00 7.68e+00
2 1 7.00e+01 6.65e+01 3.44e+00 2.31e+01 1.74e+01 5.63e+00
2 2 9.99e+01 1.00e+02 1.17eí01 1.00e+02 1.00e+02 3.75eí02
2 3 6.77e+02 4.35e+02 2.42e+02 9.79e+02 9.26e+02 5.38e+01
2 4 9.24e+03 3.78e+03 5.45e+03 9.96e+03 9.18e+03 7.76e+02
2 5 9.88e+04 3.72e+04 6.16e+04 9.99e+04 9.18e+04 8.19e+03
3 0 2.41e+02 1.45e+02 9.57e+01 3.19e+01 9.91e+00 2.20e+01
3 1 2.46e+02 1.53e+02 9.34e+01 3.91e+01 1.88e+01 2.03e+01
3 2 3.02e+02 2.30e+02 7.20e+01 1.19e+02 1.08e+02 1.11e+01
3 3 1.00e+03 1.00e+03 6.18eí02 1.00e+03 1.00e+03 4.33eí03
3 4 9.62e+03 8.70e+03 9.19e+02 9.98e+03 9.92e+03 5.96e+01
3 5 9.92e+04 8.57e+04 1.35e+04 1.00e+05 9.91e+04 8.42e+02
4 0 5.05e+02 1.67e+02 3.38e+02 4.88e+01 9.99e+00 3.88e+01
4 1 5.12e+02 1.76e+02 3.36e+02 5.65e+01 1.90e+01 3.75e+01
4 2 5.81e+02 2.65e+02 3.16e+02 1.38e+02 1.09e+02 2.96e+01
4 3 1.34e+03 1.15e+03 1.90e+02 1.02e+03 1.01e+03 1.24e+01
4 4 1.00e+04 1.00e+04 7.81eí03 1.00e+04 1.00e+04 0.00e+00
4 5 9.96e+04 9.85e+04 1.11e+03 1.00e+05 9.99e+04 6.02e+01
5 0 8.14e+02 1.70e+02 6.44e+02 6.67e+01 1.00e+01 5.67e+01
5 1 8.21e+02 1.79e+02 6.43e+02 7.47e+01 1.90e+01 5.57e+01
5 2 8.97e+02 2.69e+02 6.28e+02 1.58e+02 1.09e+02 4.90e+01
5 3 1.69e+03 1.17e+03 5.23e+02 1.04e+03 1.01e+03 3.22e+01
5 4 1.04e+04 1.02e+04 2.31e+02 1.00e+04 1.00e+04 1.26e+01
5 5 1.00e+05 1.00e+05 0.00e+00 1.00e+05 1.00e+05 0.00e+00
Parallel Implementation We have implemented the EM iteration (Eq. 3) on a DECmpp 12000Sx /
Model 200 singleíinstructionmultipleídata (SIMD) computer with a 64\Theta64 processor array. Ideally,
images would be mapped directly with one pixel per processor. Since, in these simulations, we
perform computation on images of dimension 256\Theta256 pixels, zero padded to 512\Theta512, 64 pixels
are stored per processor. Convolutions are computed with FFT algorithms tailored to take advantage
of the SIMD architecture. The additions, multiplications, and division in (3) are performed for 64 2
pixels in parallel. Bo