Linear estimator optimum for data on finite interval for restoration of image

Rino, C. L.

English

NASA Technical Reports Server

Bandlimited Image Restoration by 
Mean-Square Estimation 
Charles L. Rino 
Department of Applied E le  c t r  ophys i c s  
University of California, San Diego 
La Jolla, California, 92037 
A l inear  mean-square es t imator  optimum fo r  data available only on 
a finite interval  is derived for  the restoration of images degraded by a sys t em 
with a bandlimited spread  function. The analysis is ca r r i ed  out in one 
dimension using a prolate spheroidal wavefunction expansion of the image 
process .  
function, the expansion represents  the image process  with ze ro  mean square 
e r r o r  on the ent i re  interval,  and the estimate of the geometrical  image 
If the noise is bandlimited t o  the same bandwidth as the spread 
achieves the same mean-square e r r o r  as the optimum est imate  for data on 
the infinite interval. The rate  at which the est imate  converges is discussed 
and an example presented. 
1 
https://ntrs.nasa.gov/search.jsp?R=19690006402 2020-03-12T08:21:11+00:00Z
A solution to  the image restoration problem using l inear  mean- 
1 
p a r e  estimation has been presented by Hels t rom, who demonstrated that 
data are available over an infinite interval,  Four i e r  methods could be 
sed to  obtain the res tor ing  filter. 
quare method to  obtain a restor ing f i l t e r  optimum for  data f r o m  a sys tem 
ith a stochastic spread  function. His method uses  data only over a finite 
sgion, but he does not assume that the sys tem spread  function is band- 
mited. 
roblem for  sys t ems  with bandlimited spread functions. 
Slepian has used the l inear  mean- 
2 
Other authors have recently presented solutions t o  the restoration 
3,4,5 
In this paper a restor ing f i l ter  is derived that is optimum for  data 
vailable only over a finite interval.  The analysis is performed in one 
imension for simplicity; this imposes no fundamental limitation on the 
:chnique. 
andwidth as the spread  function, the optimum f i l te r  for  data on the finite 
i terval  achieves the same  minimum mean-square e r r o r  as the optimum 
Iter fo r  data on the infinite interval. 
In the special  case  where the noise is  bandlimited t o  the same 
The analysis makes use of the fac t  that a bandlimited stationary 
tochastic process  known only over  a finite interval can be estimated with 
e r o  mean-square e r r o r  at any point outside the interval.  
y expanding the process  in a series of prolate spheriodal wave functions 
PSWF). 
This can be done 
A simple proof is given in the Appendix. 
2 
1. The Formulation of the Problem 
The optical sys tem will be modeled mathematically by the convolu- 
onal integral  equation 
CO 
J(x) = S(x-x') J (x ' )dx '  t N(x) .  (1.1) 
0 
-CO 
(x), the illuminance of the geometr ical  image, will be taken t o  be a zero-  
3 
lean stationary 
Le point spread  
stochastic process  with spec t ra l  density cp ( w )  . S(x) is 
function of the system. The Four i e r  t r ans fo rm s (w)  of 
0 
[x) is bandlimited to  the spatial  frequency ( w  ) interval  (-R / 2 ,  s2 / 2  ). The 
ystem noise N (x) is a l so  a zero-mean stationary stochastic process  in this 
iodel. It has spec t ra l  density cp ( w ) ,  and for  simplicity it is  assumed t o  
e uncorrelated with the geometr ical  image. 
e l s t rom.  
N 
This model is discussed by 
1 
We assume that an experimenter  has observed the image-plane 
.luminance J(x) over the finite interval  (-L/2, L b ) .  His task is to  use 
lese data  t o  f o r m  an est imate  of the geometrical  image J ( x ) at some 
oint f . He does s o  by performing a l inear  operation on the available- 
0 
1 
ata that minimizes the mean-square e r r o r .  Hence, he must  find a weighting 
inction M ( 5 ,  x) such  that 
ninimizes 
3 
3ere E 5 denotes mathematical  expectation. 
At this point we could proceed by formal  methods t o  obtain an 
ntegral  equation fo r  M( 5 ,  x). There is, however, an alternative procedure 
hat makes use of the optimum filter fo r  data on the infinite interval derived 
)y Hels t rom , 
If the P S W F ' s .  
1 
First it will be necessary  to  establish'  some propert ies  
There are taken directly f r o m  Slepian and Pollak. 6 
7 The P S W F ' s  are solutions of the integral  equation 
1 
2 s in  - -R (x - x') 
Y ( c , x l )  dx' = hn(c)  Yn(c,x) (1.4) n s L'2 -L/2 'rr (x-x') 
vhere c = R  L / 4  . W e  will make use of the following two properties of the 
P R O P E R T Y  I. The P S W F ' s  Yn (c ,x)  f o r m  a complete orthonormal 
,asis for  the class of bandlimited functions. 
P R O P E R T Y  11. The P S W F ' s  Y (c,x) are orthogonal over the n 
n t e rva l  ( -L/2,  L /2 )  with 
I'hey form a complete orthogonal basis for the class of functions square 
integrable on the x-interval (,-L/2, L/2). 
4 
8 
he following equation will a lso be used, 
2. The Solution for Noise Bandlimited to  ( -  R / 2 ,  R / 2 )  
If the noise N ( x )  is  bandlimited to  the interval  (-a /2, R /2) ,  J(x ) is 
r o m  Eq. (1. I)  a bandlimited process .  
he finite interval  (-L/2, L/2), it can be estimated with ze ro  mean-square 
x r o r  at any point by its P S W F  expansion. 
Jence 
Now if  J(x) is available only over 
This is proved in the appendix. 
m - 
J (x) = 7 Jn yn ( c , x )  , 
L 
n = O  
vhe re 
L / 2  
-L/2 J n = 1 J(x') Y n (c ,x ' )  d x ' / h n ( c )  , (2.2) 
:epresents J(x) with ze ro  mean-square e r r o r  for  all x. 
inear  mean-square estimation ?(x) is equivalent t o  J(x). 
' ilter for  data on the infinite interval  can be applied t o  (2. 1). 
F o r  the purpose of 
Hence the restor ing 
W e  have 
a, m m 
L 
j o ( s ) =  M m ( s - x ' ) J ( x * ) d x l =  1 J, Mm ( 5  - x') Y, ( c , x '  )dx' 
,m 
( 2 . 3 )  n = o  
5 
:t 
r using Eq. (1. 5) and the convolution relation f o r  Four ie r  transforms, 
1. (2.3) can be simplified to 
The minimum mean square e r r o r  achieved by this estimate is 
his can be verified by substituting Eq. (2 .6 )  into Ey.  (1. 3). Note that 
i the analysis it was not necessary to assume that 5 (x) itself was band- 
mited. 
0 
Also, Eq. (2.8)  is (as i t  mus t  be from the analysis) the minimum 
lean-square e r r o r  achieved by M (x) operating on J(x) .  10 
0) 
6 
The fact  that e is  independent of 5 suggests that Eq. ( 2 . 6 )  could 
e used to estimate the illuminance of the geometrical image in  the entire 
lane when image-plane data a re  available only over the finite interval. 
'he practical  difficulties in  actually carrying this out a r e ,  however, 
Drmidable. They will be discussed in  section 4. 
3 .  The General Solution and its Mean Square E r r o r  
When the noise term in Eq. (1. 1 )  is not bandlimited to ( - Q / 2 ,  !2/2), 
he se r i e s  of Eq. (2. 1 )  no longer represents  J(x) outside the interval 
-L/2, L/2). It i s ,  however, only the noise t e r m  in Eq. (1. 1) that is  not 
epresented outside the interval ( -L /2 ,  L / 2 ) .  One might suspect f rom 
his that Eq. (2. 6) is still the optimum estimate.  This is  in  fact correct .  
t can be established by showing that Eq. (2 .6 )  does indeed sti l l  minimize 
:q. (1. 3). Now, unfortunately, the mean-square e r r o r  is no longer inde- 
,endent of 5. It contains a x-dependent t e r m  that resul ts  f rom the imper -  
ect representation of N(x)  outside ( -L/2 ,  L /2) .  
The mean square e r r o r  can be written in  the following form,  which 
Lisplays the contribution of the various terms. The function 
s used for notational convenience. 
7 
( 3 . 1 )  
ahere 
igain this resul t  can be derived by substituting Eq. (2 .6)  into Eq. ( 1 .  3 )  
rrith due regard for the different l imits  involved. The f i rs t  t e r m  in Eq. ( 3 .  1) 
s the minimum mean-square e r r o r  achieved when both the geometrical  
mage and the noise a r e  bandlimited to ( - Q / 2 ,  S2/2). 
. e rms  represent  respectively the contributions of the object process 
spectrum and the noise spectrum outside the bandwidth of the spread 
:unction. F r o m  Eq. ( 3 .  2) the third t e r m  i s  f-dependent. 
The second and third 
8 
4. Prac t ica l  Aspects of the Solution 
For  applications we need to know the rate  at which Eq. (2.6) con- 
‘erges ,  since only a finite number of terms can be evaluated. Unfortunately 
his rate is  strongly dependent on 5 even though the %-dependence of the 
nean square e r r o r  Eq. ( 3 .  1) might be small. This can be seen as follows. 
A comparison of Eq. (1. 2) and Eq. (2.6) leads to the conclusion 
hat 
f M ( 5 , x )  is represented with negligible e r r o r  by N t e r m s  in  the series (4. l ) ,  
nus t  achieve a mean-square e r r o r  close to the minimum possible e r r o r  
given by Eq. ( 3 .  1). Hence consider 
9 
03 
2rrL 
n=N+1 
(4. 3 )  
( 5 )  is the mean-square e r r o r  (not to be confused with the statist ical  
J 
ean square e r r o r )  that resul ts  when only N t e r m s  are used in  Eq. (4. 1). 
This i s  determined by W e  now want to know when 6 (5) is small. 
N 
e behavior of M (I). 
nction of w. 
*oportional to 
Now for most  applications q(w) is  a slowly varying n 
F rom Ey. (2. 7 ) ,  then, the gross  behavior of IM (5 )1  is n 
I 
u t  f rom Eq. (1. 5) this i s  proportional to h’(c) I Y  ( c ,  w) I. In Fig. 1 
n n 
(c ,x)  i s  plotted for c = 8 for n = 0 ,  5, and 8. The behavior of Y ( c , x )  
easi ly  understood when we recal l  that h (c)  represents  the fraction of 
1 n 
n 
.e total  energy in the n-th PSWF that is  in  the interval ( -L/2,  L/2). 
11 
where n = C Z C / r r l ,  
6 
.epian and Pollak have shown that for n < n cr i t ’  c r i t  
(c)  1, and for n > n h (c) 0. This behavior i s  c lear ly  displayed 
1 crit’ n 
1 Fig. 1. F r o m  this we conclude that for N > n SN(5)  will be small crit’ 
* 
lr 15 I < L/2. J ( 5 )  will achieve a minimum mean-square e r r o r  close ON 
I the theoretical value given by Eq. ( 3 .  1) for the same N and range of 5. 
10 
)r 151 > L/2 t e rms  of increasing order  will be necessary.  
5. An Example 
To tes t  the solution 1024 data points in  the interval (-1, 1) were 
Imputed by f i r s t  convolving the geometrical  image function 
.th a spread function for  a slit-diffracted image. The Fourier t ransform 
the spread function is 
( 5 . 2 )  
dependent, zero-mean gaussian random variates ,  N and N were then 
lded according to the formula 
k k'  
J = [(Jo,* Sk t Nk) 2 t Nk,  2 ,* 
k (5. 3) 
'he as te r i sk  means convolution. ) 
;ing the Cooley-Tukey Fourier  t ransform algorithm on a CDC 3600 
gital computer. 
mdom number sequence by standard methods. 
olate our initial assumptions of stationarity and additive noise. 
L t a  a r e ,  however, non-negative and serve to represent  a simple geometrical  
The convolution was car r ied  out by 
The normal  var ia tes  were  computed f rom a pseudo- 
Data so constructed 
The 
1 1  
rage. 
The coefficients J were calculated using Simpson's  rule for n = 0 n 
The PSWF' s were computed by using a Legendre 
12  
rough 8 with c = 8. 
ilynomial expansion with tabulated coefficients f rom Stratton e t  al. 
ie PSWF' s were normalized to correspond to Slepian and Pollak' s def- 
itions. 
--
6 13 
The eigenvalues A (c)  a r e  tabulated. n 
To evaluate the weights M (T), cp (w) and cp (u)) must  be specified. In n 0 N 
e absence of any pr ior  information they 
e interval ( - R / 2 ,  02/2). 
can be taken to be constant on 
14 
With these assumptions, by Eq. (2. 5),  
Now 
15 
h e r e  B is a constant that depends on the signal-to-noise ratio. 
-n 
om Eq. (2.6) Mn(5) i s  just  i / A  (c)  times n 
' q ( w ) Y  ( c ,  w). n 
With data on the interval (-1, 1) and c = 
the inverse Fourier  t ransform 
QL/4 = 8 ,  Q = 16. By 
sing these parameters  the weights M ( 5 )  were calculated for n = 0 
lrough 8 using the Cooley-Tukey algorithm. 
n 
The series (4. 2) was then 
raluated with 9 terms and the results plotted on a Calcomp automatic 
.otter. 
L Fig. 2. 
An example with ff = 0.4,  M = 1. 0 and B = 0.5 is reproduced 
12 
6. Discussion 
The example demonstrates that the estimate given by Eq. (4. 2) works 
The reason for this is c lear  i f  we observe el l  for a deterministic function. 
tat with B equal to zero  in  Eq. (5.4), Eq. (2. 5) resul ts  f rom the analytic 
mtinuation of the function J(x) to the infinite interval followed by the 
?convolution of Eq. (1. 1) with N(x) = 0. 
3 J (x) i s  effectively bandlimited to ( - ! 2 / 2 ,  sd/2),  as was the case in the 
cample. 
This can be car r ied  out a s  long 
0 
I f ,  on the other hand, the Fourier  t ransform of J (x) is much 
0 
roader than 61, J (x) will normally be nearly ze ro  outside (-L/2,  L/2). 
' data  were  available on the infinite interval,  we would have the conditions 
)r attempting super-resolution. W e  know, however, that super-resolution 
3 f rustrated by the slightest amount of noise. 3 y  4 y  
y Eq. (4. 2) the resolution is l imited to that allowed by the bandwidth of 
le spread function. 
iemselves a re  bandlimited to ( - Q / Z ,  C2/2)] 
iere is nothing new. 
0 
In  the estimate given 
[Recall section 4 and the fact  that the P S W F '  s 
Hence for super-resolution 
The conclusion seems to be that although super-resolution is not 
ossible when noise is present,  extension of an image a small distance 
utside a finite interval is possible. 
1 3  
Acknowledgments 
I wish to  thank P ro fes so r  C a r l  W. Helstrom for  his comments 
Computer time was  paid for in  pa r t  by instruction and Id cr i t ic ism.  
search  funds made available by the University of California, San Diego, 
Imputing Center and i n  pa r t  by NASA grant  NGR-05-009-079. The author 
current ly  supported by a NASA traineeship. 
14 
Appendix 
The Spheroidal Wavefunction Expansion of 
B andlimi te d Stationary Stochastic P r oc e s s e s 
A stationary stochastic process  X(x) i s  bandlimited i f  i t s  spectral  
nsity i s  non-zero only over a finite interval. 
the process  be K(x). Since 
Let  the covariance function 
e r e  cp(w) is the spectral  density of the process ,  K(x) i s  a bandlimited 
Iction. Because of properties I and I1 of the P S W F '  s (Sec. 1 .  ), the 
16 
W F  expansion 
Le r e  
K(x' - X I ' )  \y (c,x") dx' ' /h (c) 
n n P (XI) = n 
-L/2 
nverges in  the mean of order  two for --03 < x < -03 with x' fixed. For a 
a1 process  K(x-XI) = K(x' -x). Hence we also have convergence in the 
i 
lean of order  two for -m < x' < m with x fixed. 
Now f rom Eq. (1.4) the PSWF' s a r e  eigenfunctions of the Hermitian 
i tegral  operator 
J-L/ 2 
loreover,  the kernel of Eq. (A3) is square integrable over the unbounded 
quare -Q) < x , x l  < 05, and there exists a positive dominating constant M 
uch that 
OD 
or -a < x < w. In fact we have s t r i c t  equality in (A4) with M = 6a/2n. We 
urther observe that any function z(x) bandlimited to ( -$ Sa, 8 0) can be written 
sin Q ~ ( x - x '  ) 
rr(x-x' ) y(x' )dxl z(x) = 
J-L /  2 
Jhere y is a function square integrable over (-L/2,  L/2).  The represen-  
ation in Eq. (A5) and the Hermitian property of the operator L a re  suf- 
icient to establish the pointwise convergence of the se r i e s  in  Eq. (Al)  for 
.a C x < OD and x1 fixed. 
:nsures that the convergence is unifcrm \ for - m < x < 03 and x1 fixed. 
The addi t id id  property, Eq. (A4), of the kernel 7 
17 
-1- 
gain, f rom the symmetry  of K(x-x' ) the convergence is pointwise and 
i form in x1 as  well as  x. 
We can now prove the following theorem. 
IEOREM: The P S W F  expansion of the process  X(x) 
03 
N 
~ ( x )  = 1 a Y ( c , x )  
n n  
n= 0 
le r e  
inverges in  the mean fo r  all x ,  i. e. 
n= 0 
roof: Let  
N 
n= 0 
iii 
ow 
nd 
iut from Eq. (A7) and Eq. (A2) 
L/ 2 
E [a X ( x l ) ]  = ~n K(x' -x")Yn(c, X")dx"/h n ( c )  
-L/2 
Tsing the pointwise convergence of the se r i e s  in Eq. (Al) 
i V  
o evaluate the last t e r m  in  Eq. (A8) consider 
N N  
E, [%N(x)zN(x')] = 1 1 E N [ a  n m  a ] Y  n (c ,x)Y m ( c , x t )  (A12) 
m=O n=O 
sing Eq.. (A7) and Eq. (A10) 
L/2  
[ a  X ( X ~ ~ ) ] Y  (c,x")dx1I/h (c) 
N n m  n m m 
ow 
m=O n=O 
* L/2  
,B (XI!)* (C,X")dx'l/A (c) k,; m m 
'he interchange of summation and integration i s  justified by the uniform 
onvergence of the se r i e s  in Eq. (Al) .  By using Eq. (A2) this can be 
i r ther  simplified to 
V 
(c ,x“)dx”/h  (c) Y (c ,x)  
m=O m m I -  
.nally using the pointwise convergence of the ser,?s in  Eq. (A ) 
3w substituting Eq. (A14), Eq. (All), and Eq. (A9) into Eq. (A8) 
aim 2 
N+ 03 E [ I X ( X ) - ? ~ ( X ) ~  ] = 0 V X G ( - ~ ~ ,  a) 
Q. E. D. 
Dte that f rom Eq. (A13) the coefficients a a re  in general correlated.  n 
vi 
Footnote s 
, C. W. Helstrom, J. Opt. SOC. Am. 57, 297 (1967). 
, D. Slepian, J. Opt. SOC. Am. 57, 918 (1967). 
, G. J. Buck and J. J. Gustinic, IEEE Trans.  Antennas and Propagation 
-
111 
AP-15, 376 (1967). 
, B. R. Frieden, J. Opt. SOC. Am. 57, 1013 (1967). 
Lc 
, C. K. Rushforth and R. W. Har r i s ,  J. Opt. SOC. Am. 58, 539 (1968). 
I
, D. Slepian and H. 0. Pollak, Bell System Tech. J. 40, 43 (1961). - 
6 . Note that Slepian and Pollak define the frequency interval ( -0 ,  S2) 
ra ther  than (-62/2, Q/2) as in  this paper. 
0 should be replaced by hz/2. 
definitions are used directly. 
Hence in  their  formulas 
With this slight change their  
. This is derived f rom Slepian and Pollak Ref. 6 ,  Eq. (29). 
. Ref. 1, Eq. (2. 5). 
. Cf. Ref. 1, Eq. (2.6). 
. 
. J. A. Stratton, P. M. Morse,  L. J. Chu, J. D. C. Litt le,  
The expression [x] means the smallest integer l a r g e r  than x. 
c 
F. J. Corbato, Spheroidal W a v e  Functions (M.I. T. P r e s s ,  
Cambridge, Mass . ,  1956) pp. -203-229. 
. D. Slepian and E. Sonnenblick, Bell System Tech. J. 44, 1745 - 
(1965). 
Figure Captions 
g. 1. Prolate  spheroidal wave functions for  c = 8, n = 0, 5, 8. 
g. 2. Geometrical  image, computed data, and geometrical  image 
estimate with = 0.4 ,  M = 1. 0 ,  and B = 0. 5.  
4 
. 
2a67 


Band limited image restoration by linear mean-square estimation

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server