- Email: [email protected]
Deblurring Gaussian blur using a wavelet array transform
t9
Pattern Recognition, Vol. 28, No. 7, pp. 965 976, 1995
Elsevier Science Ltd Copyright @ 1995 Pattern Recognition Society Printed in Great Britain. All rights reserved 0031-3203/95 $9.50+.00
Pergamon
0031-3203(94)00146-4
D E B L U R R I N G GAUSSIAN BLUR U S I N G A WAVELET ARRAY TRANSFORM M A D H U V A I R Y a n d Y. V. V E N K A T E S H * Computer Vision and Artificial Intelligence Laboratory, Department of Electrical Engineering, Indian Institute of Science, Bangalore, India
(Received 8 March 1994; received for publication 30 November 1994) Abstract--Deblurring the Gaussian blur, which is a fundamental problem of signal analysis, has defied a satisfactory solution. The principal reason for the difficulty is the inherently ill-conditioned blur-matrix which poses a challenge to its stable inversion. Most of the literature is concerned with a solution to the problem in restricted domains, and this solution is, in many cases, characterised by inversions that are not stable, We propose a multiscale inversion method based on wavelet arrays which is applicable to a wide class of images, and show that the inversion is stable with respect to noise both in the blurred signal and in the blur variance. We include, for illustration, the result of such a deblurring scheme as applied to a natural image. Gaussian blur Hermite functions
Blur Deblur Wavelet Orthogonal representation Effective spectral width Effective spatial width
inverse o p e r a t o r J ~ - 1 such t h a t f can be recovered from h as
1. I N T R O D U C T I O N
Blurring is h a r d to avoid in any image acquisition system, a n d is caused n o t j u s t by one source, b u t by many, like (i) a t m o s p h e r i c turbulence; (ii) a n out-offocus optical system (which results in the p o i n t spread function of the imaging system n o t being a n impulse function); a n d (iii) a b e r r a t i o n s in the imaging system. Obviously, all these causes c a n n o t be simultaneously c a p t u r e d in a simple model. However, by the central limit t h e o r e m of p r o b a b i l i t y theory, the net result of the c o m p l e x i n t e r p l a y a m o n g t h e s e i n d e p e n d e n t r a n d o m sources can, in general, be a p p r o x i m a t e d by a G a u s s i a n blur. In o r d e r to c o m p e n s a t e for it, we need a technique to restore the image or, m o r e accurately, " d e b l u r " it. T o this end, we employ the G a u s s i a n model of the blurring mechanism. In view of the fact t h a t blurred images are c o m m o n , deblurring the G a u s s i a n blur is a basic p r o b l e m in image restoration, which can be formulated (for simplicity, in one dimension) as follows: Let the function h be defined as the convolution of the two f u n c t i o n s , f a n d g, of a single variable x;
h(x) = (~-f)(x) = f ( x ) * 9 ( x ) = S f ( z ) 9 ( x - ' c ) dr
(1)
-or
where 9 is the G a u s s i a n function defined by g(x)= e -x2/2~ a n d a > 0 is the variance p a r a m e t e r of the Gaussian. The deblurring p r o b l e m then is to determine f from a knowledge of the functions h a n d 9. T h e solution a m o u n t s , theoretically at least, to finding a n
f =~--1 h w i t h o u t numerical instability. Traditionally, the p r o b l e m of inverting the G a u s s i a n blur has been b r a n d e d as ill-conditioned, thereby meaning, in part, t h a t the inverse o p e r a t i o n is sensitive to m i n o r changes in the i n p u t data. In this sense, the inverse o p e r a t i o n would be stable,+ if it were possible to b o u n d some n o r m o n the error A f in the solution f in terms of some (not necessarily the same) n o r m o n the p e r t u r b a t i o n Ah in h HA/I] < CllAhll.
(3)
I n the case of the G a u s s i a n blur, it can be s h o w n t h a t the inverse o p e r a t i o n is n o t stable in the sense described above. In fact, we can find two entirely different functions which o n convolution with the same Gaussian, can result in functions whose n o r m s can be m a d e arbitrarily close, for any choice of the norm. ") 1.1. Relevant results of the literature In the literature, one finds t h a t inverting (or deblurring) the G a u s s i a n blur is modelled as a diffusion process r u n n i n g b a c k w a r d s in time, the Variance par-
~"John t2) invokes a slightly less stringent stability concept, called polynomial stability, whereby the norm of the error A f in the solution is bounded by some power of the norm of the error Ah in h. jlAfl] < CIIAhI[1/k, k finite.
* Author to whom correspondence should be addressed. 965
(2)
966
M. VAIRY and Y. V. VENKATESH
ameter of the Gaussian being treated as proportional to the time elapsed. However, note that this observation does in no way help to solve the inverse problem of deblurring Gaussian blur. One method of solving the deblurring problem consists in assuming a convolution form for the inverse operator. It is found that Fourier inversion techniques are not directly applicable, because the Fourier transform of the inverse kernel is exp ((/)20"), which cannot be the Fourier transform of any function or of a tempered distribution. Even if we attempt to achieve the deconvolution in the frequency domain, the inversion is stable only if we restrict f to belong to the space of bandlimited functions. If the original function (f) is not band-limited, the inversion is extremely sensitive to noise in the blurred signal, and it is not possible to establish the stability of the inversion. In this sense, techniques based on Fourier inversion are grossly inadequate. H u m m e l et al. ta) pose the deblurring problem in a variational framework, and constrain the spaces of the original and the blurred signals, in an attempt to convert it to a well-conditioned problem. In particular, they have shown that it is, indeed, possible to invert the Gaussian blur in a stable manner, if the class of input functions is restricted to polynomials of fixed, finite degree. They assume a convolution form for the inverse operator, so that the inversion may be realized as a filtering operation. To this end, they seek an inverse kernel, which when convolved with the blurred function restores it to the original. The inverse ( o r deblur) kernel DN(x), which involves a polynomial of degree N, restores the blurred function
f (x) = h(x) * DN(x ).
1.2. Motivation for a new framework for deblurring the Gaussian blur The motivation for the present work arose from the following considerations: • The fact the natural images cannot be satisfactorily represented by polynomials of fixed, finite degree; • The results of psycho-physical experiments which indicate that biological (including human) vision employs multi-channel processing, for low-level analysis of sensory data;* • The possibility that images can be considered to have unbounded support in both the spatial and spectral domains; I" and • The possibility of using multi-scale representations to overcome ill-conditioning of the blur-matrix. 1.3. Highlights of the new approach We propose a stable multi-scale inversion technique for a sub-space:~ of ~ ( ~ ) , which is the space of all measurable, real functions O(x) defined on the real line, satisfying
(4)
In fact, they provide an explicit construction procedure for the deblur kernel, and show that their method is a stable inversion for polynomial domains. They show that if f = p + e where p is a polynomial of degree no greater than N and e is a non-polynomial error component in f, then ]]Afll < CNIle[I
the symmetry of the Gaussian and the Toeplitz structure of the blur matrix B. Even though an analytical inversion helps to avoid error propagation in numerical computation, it does not eliminate any of the inherent ill-conditioning of the problem. The inversion may work well with high precision arithmetic only for signals with little or no noise, but is not robust with respect to noise. F o r related restoration results based on Wiener and Kalman filters, see the survey by Biemond et al. (s)
(5)
where the constant CN is explicitly shown to be dependent on N, the dimension of the polynomial domain. Images are more general than polynomials of finite degree. The restriction of the inverse operator to be a convolution (in the polynomial domain) implies that the framework proposed cannot handle natural images. In this case, if we have to realize a reasonable amount of restoration, we need to consider very high degree polynomials Ds(x ). In general, DN(x ) does not converge to any function as N tends to infinity, and the use of a deblur kernel of high degree (large N) produces a ringing effect in the reconstruction of natural images. M o r e recently, Zucker and Kimia (4) attempt an analytical inversion procedure for discrete time signals. In order to minimize the propagation of the error caused by finite precision arithmetic used in the numerical inversion of the blur matrix, the authors in (4) invoke
(1 + x 2 ) e xZJl~(x)[Zdx < 09.
(6)
-o¢
As an immediate observation, we note that LZ(R) is contained in ~ ( ~ ) , and thus the space of functions we consider includes, in the ideal case, all real world signals. The method decomposes the functions f and h into their respective multiscale approximations, and then effects the deconvolution in the frequency domain. Literature contains, what one may call, "single-step" deconvolution techniques, as in Hummel et al.,(3) where it is a convolution operation with the deblur kernel DN, or as in Zucker and Kimia (4) where it amounts to an inversion of the Toeplitz blur matrix B. As a consequence, stability in Ref. (3) has been established for a restricted class of functions (polynomials of a finite,
* Scale space representations found in the literature are based on the idea that different characteristics of a signal reveal themselves at different levels of resolution or, equivalently, in several frequency channels. See Marr. (6~ t As against the standard assumption of band-limitedness in the spectral domain and , hence, infinite spatial spread. :~A more explicit characterization of the class of functions appears to be difficult. This problem will be considered in a separate paper. It is likely that this class contains functions with appropriate smoothness restrictions.
Deblurring Gaussian blur fixed degree) while in Ref. (4) the method is not stable at all. For image representation, we choose a wavelet array composed of the generalized Hermite functions* parametrized by two scale factors (one for each spatial variable). In the proposed approach, the most important point of departure from the existing literature is that deblurring is not attempted in one step. On the contrary, we arrive at a multiscale expansion of f and g, and then achieve deconvolution at each scale stage by stage. The advantage is that we establish a stability result of the type found in Hummel et al.(3) at each scale. We also demonstrate that the deblurring the Gaussian blur is equivalent to solving simple systems of linear, simultaneous equations and is computationally efficient. Further, we can also apply the ideas of Zucker and Kimia (4) and analytically invert the matrix involved in our system of linear equations (although no LDU decomposition is needed due to the simple structure of our matrix). 1.4. Mathematical preliminaries for the new approach In what follows, R denotes the real line and R", the Cartesian n-product of R with itself. All integrals are assumed to be over the real line R, unless explicitly indicated otherwise. Z + denotes the set of non-negative integers. Let ~ ( N ) be the inner product space of all measurable, real functions O(x) defined on R, which are piecewise continuous on every finite sub-interval { - a, a} and which satisfy (1 + x2)e ~2/~10(x)l ~ dx < oo. For two functions f, g in N~(R), the inner product is given by
( f , g) = If(x)g(x) dx. The norm of ~ b ( x ) e ~ ( ~ ) is defined as
967
by ~o. The Fourier transform of f(x) will be denoted by f(o~) and is defined as follows: ~o
f(co)=
~ f(x)e-J~Xdx. ov
IIAxll x~o IIxll
The matrix norm of A is defined by IIA II = sup
2. M U L T I - S C A L E S I G N A L D E C O M P O S I T I O N A WAVELET ARRAY TRANSFORM
USING
In this section, we introduce the basic concept of a wavelet array, which is different from the standard wavelets (see, for instance, Mallat, (8) Daubechies w) and Meyer (1°)) of the literature. In the wavelet array employed here, the basis elements are the Hermite functions (of two variables, x and y) parametrized by the variance measures, al and cr2 (in the x and y directions, respectively). In the present framework, these are actually the scalefaaors (in the terminology of wavelets). The advantage of using a wavelet array representation lies in the ability to control the shape of the resolution cell. By having appropriate Hermite functions as components of a wavelet array, it is possible to simultaneously control [-within restrictions governed by the uncertainty principle(11)], the space-bandwidth product and the space-bandwidth ratio of the basis functions and arrive at fairly arbitrary tilings of the space-frequency plane [-see Refs (12) and (13)]. This advantage, not realizable in the standard wavelet tiling, can be very useful when some information is available a priori, about the space-frequency energy distribution of the signal. Another advantage of the choice of Hermite functions is that their Fourier transforms are also Hermite functions with a change of variable
I4,(x, a) = K( --j)"H,(aco, a)
(7)
h (x) = f(x)* g(x).
where K is a numerical constant. This property simplifies the deconvolution operation in the frequency domain. For simplicity in notation, consider the standard Hermite functions of one variable, H,(x), defined by ?
We write the linear span of the functions f l , f 2 .... ,f,
H,(x) = (-- 1)"e-X2 ~dx (e- x2/2)
II~,(x)II = ~ 10(x)l 2 dx. The convolution h(x), of two functions f(x) and g(x) is written as
d"
(8)
as
[-f~,f2 ..... f,q,
which can be explicitly expressed as
and its closure in ~ ( N ) as
H,(x) = n! k, ( (=~2' "0] ( - ttI(nl)t(2x)~--2t)! at) e-X2/a".
(9)
l-f1, f2 . . . . . f , q. {xJ denotes the largest integer not exceeding x, and I-x] denotes the smallest integer greater than or equal to x. ~ - - 1 is denoted by j. In general, we denote the space (or time) variable by x and the frequency variable
Generalized Hermite functions (parametrized by the scale factor a) are defined by
--Hr X H.(x, or) ~
:n!
,: o
t!(n - 20!
(lO) * For an explanation of the superiority (for certain classes of reconstruction problems) of Hermite functions, in signal representation, over other orthogonal schemes (like Gabor), see Ref. (7).
~ Extension of this notation to two variables follows the Cartesian Product approach, and is not given here.
968
M. VAIRY and Y. V. VENKATESH
More explicitly, the first few generalized Hermite functions are
Ho(x) = e-~/2~;
Hl(x ) =
2x
e ~/2~;
H2(x) = ( 4xa2 - 2 )e-~2/2a;
H3(X)
3. THE NEW DEBLURRING ALGORITHM
In this section, we present the deconvolution method and show that it is well-defined and stable. The key idea behind the deblurring method is based on rewriting equation (1) in terms of the wavelet array representations of both f and h. Consider the multi-scale Hermite expansion of f
=(\ax/~8x3 12X)e_X~/2 ~ x~
N
"
~ #P,°H,(,po )
f(x)=
n=0 N
It is known [see Refs. (14)-(16)] that the functions H~(-, ~) (n = 0, 1, 2,..., oo) for a complete, unconditional, orthogonal basis for ~ ( ~ ) . Any function f ( . ) e ~ ( ~ ) can be written as
N
+ E ~,'H,(.,p,)+...+ E ~;H.(,p~) n-0
n=O r
= f ( x ) + ~ e(x, pfl.
(17)
j-0
f(')=
L )~H,(',a),
(11)
n--0
where Y~ is obtained by an orthogonal projection of f(-) onto H,(', a) in ~ ( ~ ) .
N
( f ( ' ) , H,(', a ) ) 7,~ - < H.(', a), H.(., a) >"
The blurred signal h(x), which, form equation (1), is the blurred version of f(x), can be similarly expanded in a sequence of scales, 2j as follows:
(12)
h(~)= E ~.~°H.(,~o) n=O
The summation on the right hand side of equation (11) converges to f at every point of continuity of f. We can truncate the Hermite series for f , up to a finite number, say N, of terms and obtain an approximation f off, N
f(')= ~
7~n,(.,ao).
(13)
n=O
The error in the approximation, which is explicitly shown to be dependent on the scale a, is given by N
e ( ' , o - ) = f ( - ) - ~ ~/~H,(',a).
(14)
n--0
The approximation error in f at scale ao, 8(-, a0) is orthogortal to + Ho(., ao), H~ (', Zo),..., HN(', %) + . We now, try to approximate e(', ao) by changing the parameter a o to ~rx
N n=0
h(co) = f(o))O(a~o , a) = f(co)e ~,~,2/2.
n--0
/~(o))= f(o)e -"~2/2
~ 7~H,(',crl),
(16)
n--O
(20)
where/~(x) = 5Z,N=oakH,(x , it). Using equations (10), (17) and (18) in equation (20) and choosing 2 = p + a, we get after some simplification,
N
e(',a,)=e(',ao)--
(19)
The first level approximation f o f f i s obtained from:
(15)
where the y~ are obtained by an orthogonal projection operation onto the functions H0(', a0,H~(', ax),..., Hs(', ax). The error in this approximation defined by,
(18)
n=0
Note that, in general, there is no need to retain the same number of basis functions (N) in the Hermite expansions of f and h.§ However, in practice, we find it convenient to do so. No special advantage seems to be gained by varying the number of terms across scales or in the expansions of f and g. With minor modifications, results given below are valid for the general case, too. The convolution equation (1) assumes an algebraic form in the frequency domain
N
~(',a~)= ~ y~'U,(',a~),
N
+ ~ c~'H,(.,Z~)+...+ Z ~H,(',p.)
m-O
t,./21 (_ 1),.(20~X/p)m 2, ~mm! E t o t!(m-- 20!
is orthogonal to the space ~-Ho(',ao),Hl(',cro) .... , HN(', ao)q, as also to the space F Hi(', al) ..... HN(', al)-t. In this manner, we continue to approximate the error function until we achieve the desired degree of accuracy,3~ and arrive at a multi-scale representation for the original function f.
The right hand and the left hand sides of equation (21) are polynomials in co. Equating coefficients of equal powers of ~o, we get a system of equations of the
"~Convergenceconditions are analysedin a separate paper."21
§ In what follows, it is assumed that the blur parameter a in the blur model (l) is known.
t ~ l ( _ 1)k(2O)X/~)k 2~ , k=O
(21)
i=0
Deblurring Gaussian blur form Aix = De: N
969
• Repeat the above steps in an iterative manner, until the residual error is less than the prescribed threshold value.
ix,.m! ( - 1)m(2coxfP) p
3.1. Salient features of the algorithm
(22)
Now, we establish that the inversion at each scale is stable in the sense of equation (3), and determine the constant C in equation (3). But first we need to demonstrate that the inversion is well-defined.
Lemma 1. There exists a unique function f such that/~(x) = f ( x ) * O~(x).
for all p = 0, 1.... , N, where A = [am](N+ a)>, IN+ l) = m ! ( - 1)m0~P)P
for p = 0, 1 , . . . , N and m = 0 , 1 , . . . , N , and D is the upper triangular matrix given by
D = [dpk]lN+l)×(u+l) = k! ( - 1)k(x/2)P
For convenience, we designate by b, the vector De:
( - 1)k(,,/77)"
It can be shown (see Lemma 1 below) that the coefficients ix in the Hermite expansion of f can be uniquely determined from this system of equations. We can now identify with the approximation of the blurred function fi its unique inverse J~ At the second level, the error in the first approximation is now, in turn, approximated by changing the scale parameter in the basis functions. Again, we arrive at an equation similar to equation (22). The new, improved approximation to the solution is given by the sum of the solutions at the first and the second levels, due to the orthogonality of the approximation error to the basis functions of the preceding levels. In this manner, we can continue to obtain successively improving approximations of the solution, until we have recovered f to the desired degree of accuracy. The algorithm for deblurring is as follows: • O b t a i n the Hermite expansion of h at the scale 2j. • Using the linearity of convolution, equate the Hermite expansion of/~ to the Hermite expansion of the f (at scales P1 = 2~ -- a) multiplied by e - 2a°~2. (Due to orthogonality of the error at a particular scale to the basis functions at any of the preceding scales, the equality holds termwise.) • At each scale j, solve the system of linear equations A iu = b, where A and b are as defined above, to arrive at the coefficients # in the Hermite expansion of f at that scale. • Approximate the error in the approximation using a similar Hermite expansion, but at a different scale.
Proof. It suffices to show that the system of equations Aix = b, where A and b are as above, has a unique solution. The polynomial component of H, is even if n is even, and odd if n is odd. Therefore, only alternate Hermite functions of degree greater than or equal to p contribute in the equation involving cop. A is upper triangular and zeros alternate in every row and column. By inspection we can easily see that A has linearly independent columns (or rows), and the proof follows immediately. Alternatively, factor out (,~/p)P from the pth row and then perform the elementary row operation: ap. a m - ap+ 2. For p = 0, 1 , . . . , N - 1 . The matrix A is then diagonalized with 1 and - 1 alternating on the diagonal. det A = ( - 1)N+ 1(x/7)N~u+ 1)i2 # 0.
Remarks 1. The fact that the inversion is well-defined is also obvious from the fact that the r.h.s, of equation (22) is, after all, a polynomial of finite degree, and therefore the 1.h.s. of equation (22) is uniquely defined. Nevertheless, we shall find. Lemma 1 useful in stability and sensitivity analyses, which are now considered. We demonstrate that the proposed inversion sch _eme is stable in the sense that bounded errors in the blurred signal cause only bounded errors in the reconstruction. Theorem 1. B o u n d e d p e r t u r b a t i o n s in h cause bounded perturbations in ix. Proof. Consider a bounded perturbation of Ah in h. That is, n--N
h(.) + Ah(.) = ~, (e, + Ae,)H,(., 2) = b + Ag.
(23)
n=O
Let the corresponding perturbation in ix be Aix. We have, from the set of equations Aix = b, A(ix + Aix) = b + Ab Aix = A - l(Ab).
(24) (25)
Therefore, ]]Aixll = I I a - l A b l l
< IIA-11111Abll < oo
(26)
where ]]A ill, the matrix n o r m of A is given by the modulus of the absolutely largest eigenvalue of A - 1, 12maxI-
970
M. VAIRY and Y. V. VENKATESH
Note that 12m,xl= ~@2,~1,where "~minis the absolutely smallest eigenvalue of A.
Remark 2. While solving the system of linear equations Asu = b, we can use the ideas of Zucker and Kimia 14) and analytically invert the Aj matrix which is easy in the present case because (a) Aj is upper triangular; and (b) even in the super-diagonal entries in the matrix, zeros alternate with non-zeros in every row and every column.
where "~nax is the absolutely largest eigenvalue of
A- 1 dBA. Now, from the system of equations A# = A'~, we have,
I1~11=
IIA-1 dB A~[[ ~> ]'~min[I1~1[•
(31)
F r o m equations (30) and (31), IId#/dallo ~< 2max a
H~ll
(32)
~min
Remark 3. Zucker and Kimia (4) study the effect of changes in the blur variance on the reconstruction. Their empirical results show that when the deblurring variance parameter exactly matches the blur variance, perfect reconstruction is achieved (in the absence of noise) by using high precision arithmetic. However, it is not known if the reconstruction is robust to errors in the blur parameter, namely the variance of the Gaussian blur. We now conduct a sensitivity analysis to study the effect of perturbations in the blur parameter value on the reconstruction accuracy. It turns out that under certain conditions, it is possible to b o u n d the reconstruction error and this b o u n d is independent of the perturbations in the blur parameter o. Hence we have the following theorem:
where 2min is the absolutely smallest eigenvalue of
Theorem 2. For o <
Now, we can b o u n d the relative error IId/zl[ in/z due to an error da in a. Combining equations (32)-(34), we get,
A - 1BA. Since eigenvalues of a matrix are unaltered by a change of basis, 2m,x and 2rain are simply the absolutely largest and the absolutely smallest of the eigenvalues of dB and B, respectively. The sensitivity of p with respect to cr is directly related to the dependence of the eigenvalues of B and dB on a. p ( 1 + a ) (p/q)-I 2max =
p = N da =
/~min = ( 1 - [ - ~ t p/2 =1. P/' Ip=O
,,dp/daU
]1~11
Proof. The solution to the convolution equation is given by the solution, #, to the system of linear equations A/~ = b, with A and b as defined above. The equations in this system are
al,flij= k~o~kal,ktp) j=O
,
(27)
for p = 0 .... , N - 1 and with 2 = p + a. In matrix form, this reduces to
A# = A' ot
(28)
where A'= BA, and B=diag(fli,fl2 ..... fiN), tip=
(~-;')~/~
=
(1 + ~),,/~.
Since det A ¢ 0 (from Lemma 1), we have p = A- 1BA. Now consider a small error of da in the estimate of o. Let the corresponding error in/~ be d/~. We have d p = A - ~dB Aa
(29)
d B = diag (dflp),
(29)
where
and dflp = pfip a da. Therefore, 2p dp = A _ a dodBAzt ~< UA - 1 dffff][dBHail do = I;~La~l I[~11,
(30)
2p (3.3)
- -
where O is a constant, independent of o.
N ( + 1 P/o~(N/2)-1 da
(34)
N a ( I + ~ ) (NI2)-1 <
20
(35)
For a <
IIdp/da Na <~ I1~11
(36) 2p
F o r o <
Ud,uU/ll#N Ida/al
N ~<--. 2
(37)
Remark 4. It can be seen from equations (32) (34) above that the sensitivity of the eigenvalues of B with respect to the blur variance is independent the blur variance or. However, the sensitivity of p with respect to a varies as a polynomial in a, whose degree depends on N, the n u m b e r of basis functions at each scale. It is only possible to establish a conditional independence of the sensitivity of p with respect to a, although a stronger result would be desirable. Since the choice of N ties in our control, we believe that this conditional independence is not a serious deficiency. This shows that for o <
Deblurring Gaussian blur determine the blur parameter: In such cases, the best estimate for the blur parameter has to be adaptively estimated using some a prior~ knowledge of the scene, or by having a known test object in the scene, or by other heuristic means (like methods based on edge localization from zero-crossing contours etc.) and such estimates are prone to be noisy. In such a situation, it is imperative that the inversion be robust to perturbations in the blur parameter. The methods to tackle blur identification and estimation are, beyond the scope of this paper. For an analysis of the effect of noise on the magnitude of the coefficients, see Appendix 1. 4. E X T E N S I O N
TO SEPARABLEKERNELS
IN R a
In this section, we extend the method to separable, positive definite Gaussian blur kernels in higher dimensions. By a positive definite blur kernel in R a, we mean a function g: R d --, R, given by*
g(x) = e xTa~
5. D E B L U R R I N G W I T H A R B I T R A R Y B L U R KERNELS ON R a
In this section, we extend the formulation of the deblurring problem to non-separable kernels in two stages. First, we consider non-separable kernels defined on R. Then, using the tensor product, we extend the formulation to R d. Henceforth, by a slight abuse of terminology, we shall refer to deconvolution of any kernel as deblurring.
Non-separable kernels on R Our starting point is the same convolution equation (19). Following the general theme of the deblurring algorithm of Section 4, we obtain the multiscale Hermite expansions of f, g, and h as follows: N
f ( x ) = ~ It°.°H~(',po) n=O N n=0
R
x ~ f ( z l , z2 ..... zd)e(=' -"')/2"' dzl d'c2"" d'ca.
(40)
a(x)= T. ~°H.f',po) n=0 N
N
+ 12 fl~'H.(',a,)+...+ ~ fl;'H,(',a,) n=0
n=0
(41) N
h(x)= I2 ~.~°H.(,&) n=0 N
N
12 aP.'H.(.,p,) (42)
+ 12 ~,~'H.(',2,)+...+ n=O
n=O
Now, we continue in the same way as in the case of Gaussian deblurring in Section 3, and solve the convolution equation for the approximations off, 9, and h at the first scale. Substituting equations (43), (44) and (45) into equation (19), and setting ).j = pj + aj we get, N
tin/2] I
2., #ram: L m=o ,=o
x
Rd
R
n=0
N
h(x) = f ( x ) * g(x) = ~ f(z) g(x - z) d'c = I e(Xn-rn)2/2..... I e(x2-r2)2/2a2
N
+ 12 ,tlpl ~ . H. c ," pl)+.. + Y~ ~.H.(.,p.)
(38)
where G is a positive definite matrix of rank d. In the degenerate case where G is rank-deficient (i.e. some eigenvalues of G are zero), the blur occurs only in a p-dimensional space, where p = r(A). For remaining part of this section, we may assume without loss of generality that G is of maximal rank. We can obtain the eigenvalue decomposition of G as G = P A P T, where A is a diagonal matrix of the eigenvalues of G and P is an orthogonal matrix formed by the eigenvectors of G. Having determined the eigenvectors, we first eliminate the cross-terms in the quadratic form in the exponential by orienting it along its eigenvectors, which can be done easily by a suitable rotation of the axes. In what follows, x denotes the coordinate vector of the d-dimensional space variable whose coordinates have been measured with respect to the basis formed by the principal axes of G. In this new coordinate system, invoking Fubini's Theorem, the Gaussian blur in d-dimensions is given by
971
=
l~mt'20)
/~m-2t
- - - -
t!(m-- 2t)!
\
e
" l
]
( ~--o ln/2J(--1)k(2X/~)n-2k'~ . /~n!k=o 12 k!(n-2k)! J 2N [a/z] ( _ 1)*(2x/~)q-21 12 ~ q ] E ~. q=o l=o q!(q -- 20[ ,]
(43)
Equating coefficients of cop, we get the following system of equations:
R
(39) Thus, exploiting the separability of the Gaussian, the d-dimensional Gaussian blur can be expressed as d one-dimensionalGaussian blurs. Now, the deblurring method explained in Section 3 can be applied to each of these one-dimensional blurs. With minor modifications, we can easily see that the stability and the sensitivit3, analyses of Section 3 are also valid here.
(-- 1)("- P+~)/zm! (2~/a),-"~
(,,,
,,,,
= ( ~ o /gm(-- I)~-P)/2q'(2V/2)'~ p = 0, 1,...,2N. * In what follows,it is assumed that x is a vector in Rd,with appropriate changes for the other variables.
(44)
972
M. VAIRY and Y. V. VENKATESH
We have N + 1 unknowns in 2N equations. If there exists a unique solution to this system of equations, then N - 1 of these equations will be redundant and can be easily eliminated. Assuming for the moment that a unique inverse exists, we solve equation (47) to arrive at the first approximation of the solution /L Then, we approximate the error in the convolution equation by using a set of basis functions at a different scale, and solve the associated system of linear equations. This solution added to the solution obtained at the first stage yields the next improve estimate of the actual solution to equation (19). We continue in this m a n n e r until we have recovered f to the desired degree of accuracy. Now, we attempt to find the conditions under which an arbitrary convolution kernel is invertible. That is, we have to find out when this system of equations (one system for each scale) has a unique solution. To do this, we could form the associated matrix, as we did in Lemma 1, and find out when it is non-singular. However, the associated matrix turns out to be dense and the algebraic manipulations are very involved. Instead, we provide a simpler argument by borrowing a result from the theory of ordinary differential equations and arrive at necessary and sufficient conditions for the existence of a unique inverse. Let f ( x ) = Y ~ = o # k H k ( a , p ) SO t h a t t f(co)= Y~)'=OlZk(--j)kHk(pco, p). NOW, (19) becomes N
N
ttkHk(pco, p)(__j)k Z fltHt(aco, a)( _j)t = ~(co) k=O
1=0
(45) N
N
2 #kPk(co/4P)(--J) ke-2(p+a)(°2 Z fl~Ut(c°/x/~)(--J) t k--O
1=0
= h(co)
(46)
where Pk is the kth degree Hermite polynomial. The left band side of equation (51) can be written as a linear combination of the first N + 1 Hermite functions at scale p + a with modified coefficients #t N
conditions on f t of the form (D"f*)(Xo) = 0 for n = 0, 1. . . . . N. Hence the following theorem:
Theorem 3. Let SN(R ) be the space of all functions f defined on R and satisfying, (O~f)(xo)=O
c~
(Tf)(x) = h(x) = ~ f(z)g(x - z) dz.
/=0
(47)
We write equation (47) as
(S f*)(x) = h(x)
(50)
-oo
Then T is an isomorphism on SN(R). It is easily seen that the inversion, which we have shown to be well-defined, is also stable. In particular, we can b o u n d the n o r m of the error in h in terms of the n o r m of the error in f and the operator norm of S.
Non-separable kernels on R d The next step is to extend the method to tackle arbitrary blur kernels in R d. We simply generate wavelet bases in R d by tensor product. Before we proceed, we shall find it convenient to introduce some compact notation. In what follows, m, t, n, k, q, 1 are all multiindices in Zd+. To make the summations over multiindices meaningful, we assume an implicit total ordering (say, lexicographic) on these multi-indices. Let m = (mo, m l . . . . . md- O, [m] = ([mo], [m~] . . . . . [me 1], and ~ = o ( ' ) denote the summation of N + 1 terms (,) over the implicit (lexicographic) ordering on m. Further, let m! stand for mo !ml !"" me t !. The length of the multi-index m, Iml = mo + ml + ' " + ma- 1. A multi-index m raised to the power of another multiindex n stands for m~°m'~.... m~_-f. Let x e R a, ncZd+, a e R ~ . We generate the d-dimensional Hermite basis functions as follows:
j-d-a Hn(x, ff) =
~
Hn,(xj, aj).
(51)
j-O
#~Hk( p + a)(O, p + a)( _j)k ~, fltPt(aa), a)( --j)' =/~(co).
(49)
for some Xo ~R. Let g(x)= Z~, ~ f l , H , ( x ) , where H,(x) is the nth degree Hermite function and flo, ill,-.-, flu are constants, with fin ¢ 0. Let T be the operator given by
N
k
n=0,1 ..... N
We can write out multiscale expansion of f, g, and h similar to equations (40), (41) and (42) compactly as follows:
( ~ (48)
where the operator S=~,~_ofltD and f ( x ) = ~ = o # ~ H k ( x , p + tr). Since f * uniquely determines f, we only need to show that S is an isomorphism. S is obviously a homomorphism. To show that S is invertible, we have to establish that its kernel is the zero function. This can be obtained immediately from a set of initial
t The constant K in equation (7) is not important to the argument, and we may assume that it is absorbed in the definition of the Fourier transform of the Hermite functions.
[ m / 2 l t"- -
m=o #=m! Z,=o
l'm'2co ", ) ( NY / P')m - 2 t e-X;/2a] t.~m--2t)!
t;~] ( _ 1 ) k ( 2 ~ ) , -
2N =
/ 2k~
tq/el ( _ l),(2Xf2~)q-2,~
qE=o%Zq!,=o •
q!(q
21)!
e"
(52)
Equation (52) involves multinomials (multivariate polynomials), and leads to systems of linear equations again. In order to equate coefficients corresponding to coY; we assume that the multiscale expansion of f contributes co~and that ofg contributes cop ~. Choosing
Deblurring Gaussian blur = p + a and substituting in (19), we get, (~o
973
it appears that the Hermite functions provide a natural framework to analyse Gaussian deblurring and other related problems in signal and image processing. These functions have already been used successfully for problems in signal reconstruction using partial informationF )
p,,s~o ( - 1)1"- s)/2m' (2xfP)S']
( - 1)(" P+~)/Zm!(2~fa)P-s~
6. EXPERIMENTAL RESULTS
N
This section presents, briefly, the results of the proposed deblurring algorithm on a typical image (see Fig. 1). We demonstrate the effectiveness of the proposed deblurring algorithm on two images obtained
(_ 1)(q-p)/Zq[(2x/~lp~
(53) Naturally, solving this system of equations is more complicated than solving the corresponding system in the one-dimensional case. But apart from an increase in complexity by a factor of d, the associated matrix of equation (A3) is not upper triangular, and is in fact quite dense. In general, solving systems of linear equations has complexity O(N3). While the sparse nature of the matrix of equation (20) implies that solving equation (20) has O(N 2) complexity, the complexity of solving equation (A3) would be O(dN3). The extension to arbitrary kernels has several important applications in the fields of signal and image processing. For instance, a well known task in image processing is edge detection. Edges in an image are points of maximal gradient, and are obtained by setting the Laplacian of the image to zero. To validate the differentiation operation and reduce the effect of noise, the image is first convolved with a suitable Gaussian and the zero-crossings of the Laplacian of this new image (the so-called LoG zero-crossing contours) are points of maximum intensity variation in the original image. It turns out that the Laplacian of the Gaussian is just the second degree Hermite function H2(x ). Thus,
Fig. 1. Original image of GIRL.
Fig. 2(a). Blurred image: ~x = 1.0, ~ry= 1.0.
Fig. 2(b). Restored version of image shown in Fig. 2(a).
974
M. VAIRY and Y. V. VENKATESH
Fig. 3(a). Blurred image: ox = 2.0, ~ry= 2.0.
Fig. 3(b). Restored version of image shown in Fig. 3(a).
Fig. 4(a). Blurred image with noise: crx = 1.0, cry = 1.0. The noise variance is 1000/12 in both the x and the y directions,
Fig. 4(b). Restored version of the blurred and noisy image shown in Fig. 4(a).
by blurring (corresponding to two different values of ~r) the image of Fig. 1. In the implementation, the reconstruction procedure used the first twenty Hermite functions at twenty different scales. The original image in Fig. 1, blurred by an isotropic Gaussian of variance 1.0, is shown in Fig. 2(a). The reconstruction is shown in Fig. 2(b). The second case is a more severe blur with an isotropic Gaussian of variance 2.0, as can be seen in Fig. 3(a). The reconstruction is given in Fig. 3(b). It is observed that multiscale inversion using the proposed deblurring algorithm has restored most of the image sharpness.
The effect of noise is shown in Fig. 4(a), where the original image has been corrupted by uncorrelated, uniformly distributed noise with zero mean. The reconstruction is shown in Fig. 4(b). It can be seen that the reconstruction is stable with respect to noise. 7. C O N C L U S I O N S
The problem of deblurring the Gaussian blur is considered in the framework of a wavelet array transform for representing both the given (blurred) and the unknown images. A new method is proposed for de-
Deblurring Gaussian blur blurring, in a multi-scale framework, based on the use of generalized Hermite functions. The deblurring problem is reduced to solving sparse systems of simultaneous, linear equations, which can be done analytically, and is computationally efficient. The method directly extends to higher dimensions by making use of the separability of tile Gaussian kernel, and the deblurring operation on Euclidean d-space reduces to d independent deblurring operations on R. We have also formulated the problem for nonseparable kernels on R. While earlier attempts in the literature have either suffered from lack of stability or have established stable inversion procedures for restricted domains, we have established the stability of the inversion for perturbations in the blurred signal and also (conditionally)for errors in the choice of the blur parameter for a wide class of input images. The proposed procedure is robust with respect to noise, and seems to work reasonably well even in cases where the blur parameter in not known a priori but has to be estimated by other means. Experimental results for a natural image show that the method is superior to the results of the literature in terms of (i) robustness; and (ii) generality of the assumptions made on the input images. Apart from deblurring for the purpose of signal or image restoration, the m e t h o d can also be used for blur identification, where a known test signal f is fed to the blur mechanism (for which we can arrive at a multiscale Hermite expansion), and the result of the blurring is observed. Then, multiscale deconvolution can be used to estimate the blur kernel and its parameters. Future scope for research in this area includes extensions of d this method to the whole o f L 2 (R), and to non-stationary blurs. Spatially varying blurs are common in camera defocus blurs, where the blur kernel's variance varies with the distance of the object from the optic centre of the camera lens. Given the camera parameters, a deconvolution method for such spatially varying blurs would provide depth information directly, thereby dispensing with the correspondence problem in stereopsis.
975
6. D. Marr, Vision. Freeman, San Francisco, California (1982). 7. Y. V. Venkatesh, Hermite polynomials for signal reconstruction from zero-crossings--part h one dimensional signals, IEF Proc. I 139(6), 587-596 (1992). 8. S. M. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2(R), Trans. Am. Math. Soc. 315(1), 69-87 (1989). 9. I. Daubechies, Ten Lectures on Wavelets. SIAM, Philadelphia (1992). 10. Y. Meyer, Ondelettes et Operateurs--Tomes I e t II. Hermann, Paris (1990). 11. Bruijn, Uncertainty principle in Fourier analysis, in Inequalities, O. Shisha, ed., pp. 57-71. Academic Press, New Yrok (1965). 12. Y. V. Venkatesh, K. Ramani and R. Nandini, Hermite sieve as a wavelet array for 1-D and 2-D signal decomposition, lEE Proc. I Vis. Image Signal Process. 141, 348-356 _(1994). 13. R. Nandini, Wavelet array transform for signal representation using generalized Hermite polynomials, M.S. Thesis, Department of Electrical Engineering, Indian Institute of Science (1992). 14. G. Szego, Orthogonal polynomials, Am. Math. Soc. (1975). 15. N.N. Lebedev, Special Functions and their Applications. Prentice-Hall, Englewood Cliffs, New Jersy (1965). 16. J. R. Higgins, Completeness and Basis Properties of Sets of Special Functions, Cambridge Tracts in Mathematics. Cambridge University Press, Cambridge (1977). 17. Madhu Vairy and Y, V. Venkatesh, Wavelet array representation for deblurring Gaussian blur, Technical Reports Department of Electrical Engineering, Indian Institute of Science, Bangalore, India (1993). 18. Gnedenko, Theory of Probability. MIR, Moscow (1975). APPENDIX
Noise Analysis
We briefly analyse the robustness of the method to noise. Let q(x) be (broad sense) stationary, uncorrelated additive noise with zero mean and variance a 2. At a fixed scale, let 6#k (a random variable)* be the change in the coefficient #k due the additive noise r/(x). go
6/~k= S q(T)Hk(T)dv"
(A1)
~c
If E denotes the expectation operator, then assuming interchangeability of.the integral and the expectation, we get, oD
E[6#k] = Acknowledgements--MV would like to thank Shanmukh, Ganesh Murthy and Chris Heil for several stimulating digcussions. REFERENCES
I. A. N. Tikhonov and V. Y. Arsenin, Solutions of ill-posed problems, F. John, ed., Wiley, New York (1977). 2. F. John, Numerical solution of the equation of heat conduction at preceding times, Ann. Math. 17, 129-142 (1955). 3. R. A. Hummel, B. B. Kimia and S. W. Zucker, Deblurring Gaussian blur, Comput. Vision Graphics Image Process. 38, 66-80 (1987). 4. B. B. Kimia and S. W. Zucker, Analytic inverse of the discrete Gaussian blur, Opt. Engng. 32(1), 166-176 (1993). 5. J. Biemond, R. L. Lagendijk and R. M. Mersereau, Iterative methods for image deblurring, Proc. IEEE 78(5), 856-890 (1990).
1 _ ~ E[q(z)JHk(T)dT=O. IIHkl]~ -
(A2)
Similarly, E[6/t~] = ~
1
~
cc
_~ -~f E[q('Orl('c')JHk('c)Hk('r')d'cd'V
O-2
= iii/4k 1~.
(h3)
Therefore, the variance of the kth coefficient at any scale is decreased by a factor of IIHkII2, where IIHII2 = 2k IIHk- a II2 for k = 1, 2, 3..... and Ilno]]2 = x//~.
* Under certain conditions, the integral in equation (38) converges to a random variable, according to Khinchin's theorem. See Ref. [18] for the proof of Khinchin's theorem.
976
M. VAIRY and Y. V. V E N K A T E S H
About the Author--Y. V. V E N K A T E S H got his Ph.D. from the Indian Institute of Science for a dissertation on the stability analysis of feedback systems. He was an Alexander von H u m b o l d t Fellow at the Universities of Karlsruhe, Freiburg and Erlangen, Germany; National Research Council Fellow (U.S.A.) at the Goddard Space Flight Center, Greenbelt, Maryland; Visiting Fellow at the Australian National University, to name a few. His research m o n o g r a p h on stability and instability analysis of linear-nonlinear time-varying feedback systems has appeared in the Springer Verlag Physics Lecture Notes Series. His present work is on signal representation using wavelet-link arrays, and reconstruction from partial information (like zerocrossings). He is a Professor at the Indian Institute of Science, Bangalore, India, and currently the Chairman of the Department of Electrical Engineering. He is a Fellow of the Indian National Science Academy, Indian Academy of Sciences and the Indian Academy of Engineering.
About the Author--M. VAIRY studied Electrical and Electronic Engineering at the Birla Institute of Technology, Pilani, Rajasthan, India. At the time of submission of the paper, he was a Research Scholar (M.Sc.) at the Indian Institute of Science, Bangalore, India. Now he is a Research Assistant at the Courant Institute, New York University, New York, U.S.A.
Pattern Recognition, Vol. 28, No. 7, pp. 965 976, 1995
Elsevier Science Ltd Copyright @ 1995 Pattern Recognition Society Printed in Great Britain. All rights reserved 0031-3203/95 $9.50+.00
Pergamon
0031-3203(94)00146-4
D E B L U R R I N G GAUSSIAN BLUR U S I N G A WAVELET ARRAY TRANSFORM M A D H U V A I R Y a n d Y. V. V E N K A T E S H * Computer Vision and Artificial Intelligence Laboratory, Department of Electrical Engineering, Indian Institute of Science, Bangalore, India
(Received 8 March 1994; received for publication 30 November 1994) Abstract--Deblurring the Gaussian blur, which is a fundamental problem of signal analysis, has defied a satisfactory solution. The principal reason for the difficulty is the inherently ill-conditioned blur-matrix which poses a challenge to its stable inversion. Most of the literature is concerned with a solution to the problem in restricted domains, and this solution is, in many cases, characterised by inversions that are not stable, We propose a multiscale inversion method based on wavelet arrays which is applicable to a wide class of images, and show that the inversion is stable with respect to noise both in the blurred signal and in the blur variance. We include, for illustration, the result of such a deblurring scheme as applied to a natural image. Gaussian blur Hermite functions
Blur Deblur Wavelet Orthogonal representation Effective spectral width Effective spatial width
inverse o p e r a t o r J ~ - 1 such t h a t f can be recovered from h as
1. I N T R O D U C T I O N
Blurring is h a r d to avoid in any image acquisition system, a n d is caused n o t j u s t by one source, b u t by many, like (i) a t m o s p h e r i c turbulence; (ii) a n out-offocus optical system (which results in the p o i n t spread function of the imaging system n o t being a n impulse function); a n d (iii) a b e r r a t i o n s in the imaging system. Obviously, all these causes c a n n o t be simultaneously c a p t u r e d in a simple model. However, by the central limit t h e o r e m of p r o b a b i l i t y theory, the net result of the c o m p l e x i n t e r p l a y a m o n g t h e s e i n d e p e n d e n t r a n d o m sources can, in general, be a p p r o x i m a t e d by a G a u s s i a n blur. In o r d e r to c o m p e n s a t e for it, we need a technique to restore the image or, m o r e accurately, " d e b l u r " it. T o this end, we employ the G a u s s i a n model of the blurring mechanism. In view of the fact t h a t blurred images are c o m m o n , deblurring the G a u s s i a n blur is a basic p r o b l e m in image restoration, which can be formulated (for simplicity, in one dimension) as follows: Let the function h be defined as the convolution of the two f u n c t i o n s , f a n d g, of a single variable x;
h(x) = (~-f)(x) = f ( x ) * 9 ( x ) = S f ( z ) 9 ( x - ' c ) dr
(1)
-or
where 9 is the G a u s s i a n function defined by g(x)= e -x2/2~ a n d a > 0 is the variance p a r a m e t e r of the Gaussian. The deblurring p r o b l e m then is to determine f from a knowledge of the functions h a n d 9. T h e solution a m o u n t s , theoretically at least, to finding a n
f =~--1 h w i t h o u t numerical instability. Traditionally, the p r o b l e m of inverting the G a u s s i a n blur has been b r a n d e d as ill-conditioned, thereby meaning, in part, t h a t the inverse o p e r a t i o n is sensitive to m i n o r changes in the i n p u t data. In this sense, the inverse o p e r a t i o n would be stable,+ if it were possible to b o u n d some n o r m o n the error A f in the solution f in terms of some (not necessarily the same) n o r m o n the p e r t u r b a t i o n Ah in h HA/I] < CllAhll.
(3)
I n the case of the G a u s s i a n blur, it can be s h o w n t h a t the inverse o p e r a t i o n is n o t stable in the sense described above. In fact, we can find two entirely different functions which o n convolution with the same Gaussian, can result in functions whose n o r m s can be m a d e arbitrarily close, for any choice of the norm. ") 1.1. Relevant results of the literature In the literature, one finds t h a t inverting (or deblurring) the G a u s s i a n blur is modelled as a diffusion process r u n n i n g b a c k w a r d s in time, the Variance par-
~"John t2) invokes a slightly less stringent stability concept, called polynomial stability, whereby the norm of the error A f in the solution is bounded by some power of the norm of the error Ah in h. jlAfl] < CIIAhI[1/k, k finite.
* Author to whom correspondence should be addressed. 965
(2)
966
M. VAIRY and Y. V. VENKATESH
ameter of the Gaussian being treated as proportional to the time elapsed. However, note that this observation does in no way help to solve the inverse problem of deblurring Gaussian blur. One method of solving the deblurring problem consists in assuming a convolution form for the inverse operator. It is found that Fourier inversion techniques are not directly applicable, because the Fourier transform of the inverse kernel is exp ((/)20"), which cannot be the Fourier transform of any function or of a tempered distribution. Even if we attempt to achieve the deconvolution in the frequency domain, the inversion is stable only if we restrict f to belong to the space of bandlimited functions. If the original function (f) is not band-limited, the inversion is extremely sensitive to noise in the blurred signal, and it is not possible to establish the stability of the inversion. In this sense, techniques based on Fourier inversion are grossly inadequate. H u m m e l et al. ta) pose the deblurring problem in a variational framework, and constrain the spaces of the original and the blurred signals, in an attempt to convert it to a well-conditioned problem. In particular, they have shown that it is, indeed, possible to invert the Gaussian blur in a stable manner, if the class of input functions is restricted to polynomials of fixed, finite degree. They assume a convolution form for the inverse operator, so that the inversion may be realized as a filtering operation. To this end, they seek an inverse kernel, which when convolved with the blurred function restores it to the original. The inverse ( o r deblur) kernel DN(x), which involves a polynomial of degree N, restores the blurred function
f (x) = h(x) * DN(x ).
1.2. Motivation for a new framework for deblurring the Gaussian blur The motivation for the present work arose from the following considerations: • The fact the natural images cannot be satisfactorily represented by polynomials of fixed, finite degree; • The results of psycho-physical experiments which indicate that biological (including human) vision employs multi-channel processing, for low-level analysis of sensory data;* • The possibility that images can be considered to have unbounded support in both the spatial and spectral domains; I" and • The possibility of using multi-scale representations to overcome ill-conditioning of the blur-matrix. 1.3. Highlights of the new approach We propose a stable multi-scale inversion technique for a sub-space:~ of ~ ( ~ ) , which is the space of all measurable, real functions O(x) defined on the real line, satisfying
(4)
In fact, they provide an explicit construction procedure for the deblur kernel, and show that their method is a stable inversion for polynomial domains. They show that if f = p + e where p is a polynomial of degree no greater than N and e is a non-polynomial error component in f, then ]]Afll < CNIle[I
the symmetry of the Gaussian and the Toeplitz structure of the blur matrix B. Even though an analytical inversion helps to avoid error propagation in numerical computation, it does not eliminate any of the inherent ill-conditioning of the problem. The inversion may work well with high precision arithmetic only for signals with little or no noise, but is not robust with respect to noise. F o r related restoration results based on Wiener and Kalman filters, see the survey by Biemond et al. (s)
(5)
where the constant CN is explicitly shown to be dependent on N, the dimension of the polynomial domain. Images are more general than polynomials of finite degree. The restriction of the inverse operator to be a convolution (in the polynomial domain) implies that the framework proposed cannot handle natural images. In this case, if we have to realize a reasonable amount of restoration, we need to consider very high degree polynomials Ds(x ). In general, DN(x ) does not converge to any function as N tends to infinity, and the use of a deblur kernel of high degree (large N) produces a ringing effect in the reconstruction of natural images. M o r e recently, Zucker and Kimia (4) attempt an analytical inversion procedure for discrete time signals. In order to minimize the propagation of the error caused by finite precision arithmetic used in the numerical inversion of the blur matrix, the authors in (4) invoke
(1 + x 2 ) e xZJl~(x)[Zdx < 09.
(6)
-o¢
As an immediate observation, we note that LZ(R) is contained in ~ ( ~ ) , and thus the space of functions we consider includes, in the ideal case, all real world signals. The method decomposes the functions f and h into their respective multiscale approximations, and then effects the deconvolution in the frequency domain. Literature contains, what one may call, "single-step" deconvolution techniques, as in Hummel et al.,(3) where it is a convolution operation with the deblur kernel DN, or as in Zucker and Kimia (4) where it amounts to an inversion of the Toeplitz blur matrix B. As a consequence, stability in Ref. (3) has been established for a restricted class of functions (polynomials of a finite,
* Scale space representations found in the literature are based on the idea that different characteristics of a signal reveal themselves at different levels of resolution or, equivalently, in several frequency channels. See Marr. (6~ t As against the standard assumption of band-limitedness in the spectral domain and , hence, infinite spatial spread. :~A more explicit characterization of the class of functions appears to be difficult. This problem will be considered in a separate paper. It is likely that this class contains functions with appropriate smoothness restrictions.
Deblurring Gaussian blur fixed degree) while in Ref. (4) the method is not stable at all. For image representation, we choose a wavelet array composed of the generalized Hermite functions* parametrized by two scale factors (one for each spatial variable). In the proposed approach, the most important point of departure from the existing literature is that deblurring is not attempted in one step. On the contrary, we arrive at a multiscale expansion of f and g, and then achieve deconvolution at each scale stage by stage. The advantage is that we establish a stability result of the type found in Hummel et al.(3) at each scale. We also demonstrate that the deblurring the Gaussian blur is equivalent to solving simple systems of linear, simultaneous equations and is computationally efficient. Further, we can also apply the ideas of Zucker and Kimia (4) and analytically invert the matrix involved in our system of linear equations (although no LDU decomposition is needed due to the simple structure of our matrix). 1.4. Mathematical preliminaries for the new approach In what follows, R denotes the real line and R", the Cartesian n-product of R with itself. All integrals are assumed to be over the real line R, unless explicitly indicated otherwise. Z + denotes the set of non-negative integers. Let ~ ( N ) be the inner product space of all measurable, real functions O(x) defined on R, which are piecewise continuous on every finite sub-interval { - a, a} and which satisfy (1 + x2)e ~2/~10(x)l ~ dx < oo. For two functions f, g in N~(R), the inner product is given by
( f , g) = If(x)g(x) dx. The norm of ~ b ( x ) e ~ ( ~ ) is defined as
967
by ~o. The Fourier transform of f(x) will be denoted by f(o~) and is defined as follows: ~o
f(co)=
~ f(x)e-J~Xdx. ov
IIAxll x~o IIxll
The matrix norm of A is defined by IIA II = sup
2. M U L T I - S C A L E S I G N A L D E C O M P O S I T I O N A WAVELET ARRAY TRANSFORM
USING
In this section, we introduce the basic concept of a wavelet array, which is different from the standard wavelets (see, for instance, Mallat, (8) Daubechies w) and Meyer (1°)) of the literature. In the wavelet array employed here, the basis elements are the Hermite functions (of two variables, x and y) parametrized by the variance measures, al and cr2 (in the x and y directions, respectively). In the present framework, these are actually the scalefaaors (in the terminology of wavelets). The advantage of using a wavelet array representation lies in the ability to control the shape of the resolution cell. By having appropriate Hermite functions as components of a wavelet array, it is possible to simultaneously control [-within restrictions governed by the uncertainty principle(11)], the space-bandwidth product and the space-bandwidth ratio of the basis functions and arrive at fairly arbitrary tilings of the space-frequency plane [-see Refs (12) and (13)]. This advantage, not realizable in the standard wavelet tiling, can be very useful when some information is available a priori, about the space-frequency energy distribution of the signal. Another advantage of the choice of Hermite functions is that their Fourier transforms are also Hermite functions with a change of variable
I4,(x, a) = K( --j)"H,(aco, a)
(7)
h (x) = f(x)* g(x).
where K is a numerical constant. This property simplifies the deconvolution operation in the frequency domain. For simplicity in notation, consider the standard Hermite functions of one variable, H,(x), defined by ?
We write the linear span of the functions f l , f 2 .... ,f,
H,(x) = (-- 1)"e-X2 ~dx (e- x2/2)
II~,(x)II = ~ 10(x)l 2 dx. The convolution h(x), of two functions f(x) and g(x) is written as
d"
(8)
as
[-f~,f2 ..... f,q,
which can be explicitly expressed as
and its closure in ~ ( N ) as
H,(x) = n! k, ( (=~2' "0] ( - ttI(nl)t(2x)~--2t)! at) e-X2/a".
(9)
l-f1, f2 . . . . . f , q. {xJ denotes the largest integer not exceeding x, and I-x] denotes the smallest integer greater than or equal to x. ~ - - 1 is denoted by j. In general, we denote the space (or time) variable by x and the frequency variable
Generalized Hermite functions (parametrized by the scale factor a) are defined by
--Hr X H.(x, or) ~
:n!
,: o
t!(n - 20!
(lO) * For an explanation of the superiority (for certain classes of reconstruction problems) of Hermite functions, in signal representation, over other orthogonal schemes (like Gabor), see Ref. (7).
~ Extension of this notation to two variables follows the Cartesian Product approach, and is not given here.
968
M. VAIRY and Y. V. VENKATESH
More explicitly, the first few generalized Hermite functions are
Ho(x) = e-~/2~;
Hl(x ) =
2x
e ~/2~;
H2(x) = ( 4xa2 - 2 )e-~2/2a;
H3(X)
3. THE NEW DEBLURRING ALGORITHM
In this section, we present the deconvolution method and show that it is well-defined and stable. The key idea behind the deblurring method is based on rewriting equation (1) in terms of the wavelet array representations of both f and h. Consider the multi-scale Hermite expansion of f
=(\ax/~8x3 12X)e_X~/2 ~ x~
N
"
~ #P,°H,(,po )
f(x)=
n=0 N
It is known [see Refs. (14)-(16)] that the functions H~(-, ~) (n = 0, 1, 2,..., oo) for a complete, unconditional, orthogonal basis for ~ ( ~ ) . Any function f ( . ) e ~ ( ~ ) can be written as
N
+ E ~,'H,(.,p,)+...+ E ~;H.(,p~) n-0
n=O r
= f ( x ) + ~ e(x, pfl.
(17)
j-0
f(')=
L )~H,(',a),
(11)
n--0
where Y~ is obtained by an orthogonal projection of f(-) onto H,(', a) in ~ ( ~ ) .
N
( f ( ' ) , H,(', a ) ) 7,~ - < H.(', a), H.(., a) >"
The blurred signal h(x), which, form equation (1), is the blurred version of f(x), can be similarly expanded in a sequence of scales, 2j as follows:
(12)
h(~)= E ~.~°H.(,~o) n=O
The summation on the right hand side of equation (11) converges to f at every point of continuity of f. We can truncate the Hermite series for f , up to a finite number, say N, of terms and obtain an approximation f off, N
f(')= ~
7~n,(.,ao).
(13)
n=O
The error in the approximation, which is explicitly shown to be dependent on the scale a, is given by N
e ( ' , o - ) = f ( - ) - ~ ~/~H,(',a).
(14)
n--0
The approximation error in f at scale ao, 8(-, a0) is orthogortal to + Ho(., ao), H~ (', Zo),..., HN(', %) + . We now, try to approximate e(', ao) by changing the parameter a o to ~rx
N n=0
h(co) = f(o))O(a~o , a) = f(co)e ~,~,2/2.
n--0
/~(o))= f(o)e -"~2/2
~ 7~H,(',crl),
(16)
n--O
(20)
where/~(x) = 5Z,N=oakH,(x , it). Using equations (10), (17) and (18) in equation (20) and choosing 2 = p + a, we get after some simplification,
N
e(',a,)=e(',ao)--
(19)
The first level approximation f o f f i s obtained from:
(15)
where the y~ are obtained by an orthogonal projection operation onto the functions H0(', a0,H~(', ax),..., Hs(', ax). The error in this approximation defined by,
(18)
n=0
Note that, in general, there is no need to retain the same number of basis functions (N) in the Hermite expansions of f and h.§ However, in practice, we find it convenient to do so. No special advantage seems to be gained by varying the number of terms across scales or in the expansions of f and g. With minor modifications, results given below are valid for the general case, too. The convolution equation (1) assumes an algebraic form in the frequency domain
N
~(',a~)= ~ y~'U,(',a~),
N
+ ~ c~'H,(.,Z~)+...+ Z ~H,(',p.)
m-O
t,./21 (_ 1),.(20~X/p)m 2, ~mm! E t o t!(m-- 20!
is orthogonal to the space ~-Ho(',ao),Hl(',cro) .... , HN(', ao)q, as also to the space F Hi(', al) ..... HN(', al)-t. In this manner, we continue to approximate the error function until we achieve the desired degree of accuracy,3~ and arrive at a multi-scale representation for the original function f.
The right hand and the left hand sides of equation (21) are polynomials in co. Equating coefficients of equal powers of ~o, we get a system of equations of the
"~Convergenceconditions are analysedin a separate paper."21
§ In what follows, it is assumed that the blur parameter a in the blur model (l) is known.
t ~ l ( _ 1)k(2O)X/~)k 2~ , k=O
(21)
i=0
Deblurring Gaussian blur form Aix = De: N
969
• Repeat the above steps in an iterative manner, until the residual error is less than the prescribed threshold value.
ix,.m! ( - 1)m(2coxfP) p
3.1. Salient features of the algorithm
(22)
Now, we establish that the inversion at each scale is stable in the sense of equation (3), and determine the constant C in equation (3). But first we need to demonstrate that the inversion is well-defined.
Lemma 1. There exists a unique function f such that/~(x) = f ( x ) * O~(x).
for all p = 0, 1.... , N, where A = [am](N+ a)>, IN+ l) = m ! ( - 1)m0~P)P
for p = 0, 1 , . . . , N and m = 0 , 1 , . . . , N , and D is the upper triangular matrix given by
D = [dpk]lN+l)×(u+l) = k! ( - 1)k(x/2)P
For convenience, we designate by b, the vector De:
( - 1)k(,,/77)"
It can be shown (see Lemma 1 below) that the coefficients ix in the Hermite expansion of f can be uniquely determined from this system of equations. We can now identify with the approximation of the blurred function fi its unique inverse J~ At the second level, the error in the first approximation is now, in turn, approximated by changing the scale parameter in the basis functions. Again, we arrive at an equation similar to equation (22). The new, improved approximation to the solution is given by the sum of the solutions at the first and the second levels, due to the orthogonality of the approximation error to the basis functions of the preceding levels. In this manner, we can continue to obtain successively improving approximations of the solution, until we have recovered f to the desired degree of accuracy. The algorithm for deblurring is as follows: • O b t a i n the Hermite expansion of h at the scale 2j. • Using the linearity of convolution, equate the Hermite expansion of/~ to the Hermite expansion of the f (at scales P1 = 2~ -- a) multiplied by e - 2a°~2. (Due to orthogonality of the error at a particular scale to the basis functions at any of the preceding scales, the equality holds termwise.) • At each scale j, solve the system of linear equations A iu = b, where A and b are as defined above, to arrive at the coefficients # in the Hermite expansion of f at that scale. • Approximate the error in the approximation using a similar Hermite expansion, but at a different scale.
Proof. It suffices to show that the system of equations Aix = b, where A and b are as above, has a unique solution. The polynomial component of H, is even if n is even, and odd if n is odd. Therefore, only alternate Hermite functions of degree greater than or equal to p contribute in the equation involving cop. A is upper triangular and zeros alternate in every row and column. By inspection we can easily see that A has linearly independent columns (or rows), and the proof follows immediately. Alternatively, factor out (,~/p)P from the pth row and then perform the elementary row operation: ap. a m - ap+ 2. For p = 0, 1 , . . . , N - 1 . The matrix A is then diagonalized with 1 and - 1 alternating on the diagonal. det A = ( - 1)N+ 1(x/7)N~u+ 1)i2 # 0.
Remarks 1. The fact that the inversion is well-defined is also obvious from the fact that the r.h.s, of equation (22) is, after all, a polynomial of finite degree, and therefore the 1.h.s. of equation (22) is uniquely defined. Nevertheless, we shall find. Lemma 1 useful in stability and sensitivity analyses, which are now considered. We demonstrate that the proposed inversion sch _eme is stable in the sense that bounded errors in the blurred signal cause only bounded errors in the reconstruction. Theorem 1. B o u n d e d p e r t u r b a t i o n s in h cause bounded perturbations in ix. Proof. Consider a bounded perturbation of Ah in h. That is, n--N
h(.) + Ah(.) = ~, (e, + Ae,)H,(., 2) = b + Ag.
(23)
n=O
Let the corresponding perturbation in ix be Aix. We have, from the set of equations Aix = b, A(ix + Aix) = b + Ab Aix = A - l(Ab).
(24) (25)
Therefore, ]]Aixll = I I a - l A b l l
< IIA-11111Abll < oo
(26)
where ]]A ill, the matrix n o r m of A is given by the modulus of the absolutely largest eigenvalue of A - 1, 12maxI-
970
M. VAIRY and Y. V. VENKATESH
Note that 12m,xl= ~@2,~1,where "~minis the absolutely smallest eigenvalue of A.
Remark 2. While solving the system of linear equations Asu = b, we can use the ideas of Zucker and Kimia 14) and analytically invert the Aj matrix which is easy in the present case because (a) Aj is upper triangular; and (b) even in the super-diagonal entries in the matrix, zeros alternate with non-zeros in every row and every column.
where "~nax is the absolutely largest eigenvalue of
A- 1 dBA. Now, from the system of equations A# = A'~, we have,
I1~11=
IIA-1 dB A~[[ ~> ]'~min[I1~1[•
(31)
F r o m equations (30) and (31), IId#/dallo ~< 2max a
H~ll
(32)
~min
Remark 3. Zucker and Kimia (4) study the effect of changes in the blur variance on the reconstruction. Their empirical results show that when the deblurring variance parameter exactly matches the blur variance, perfect reconstruction is achieved (in the absence of noise) by using high precision arithmetic. However, it is not known if the reconstruction is robust to errors in the blur parameter, namely the variance of the Gaussian blur. We now conduct a sensitivity analysis to study the effect of perturbations in the blur parameter value on the reconstruction accuracy. It turns out that under certain conditions, it is possible to b o u n d the reconstruction error and this b o u n d is independent of the perturbations in the blur parameter o. Hence we have the following theorem:
where 2min is the absolutely smallest eigenvalue of
Theorem 2. For o <
Now, we can b o u n d the relative error IId/zl[ in/z due to an error da in a. Combining equations (32)-(34), we get,
A - 1BA. Since eigenvalues of a matrix are unaltered by a change of basis, 2m,x and 2rain are simply the absolutely largest and the absolutely smallest of the eigenvalues of dB and B, respectively. The sensitivity of p with respect to cr is directly related to the dependence of the eigenvalues of B and dB on a. p ( 1 + a ) (p/q)-I 2max =
p = N da =
/~min = ( 1 - [ - ~ t p/2 =1. P/' Ip=O
,,dp/daU
]1~11
Proof. The solution to the convolution equation is given by the solution, #, to the system of linear equations A/~ = b, with A and b as defined above. The equations in this system are
al,flij= k~o~kal,ktp) j=O
,
(27)
for p = 0 .... , N - 1 and with 2 = p + a. In matrix form, this reduces to
A# = A' ot
(28)
where A'= BA, and B=diag(fli,fl2 ..... fiN), tip=
(~-;')~/~
=
(1 + ~),,/~.
Since det A ¢ 0 (from Lemma 1), we have p = A- 1BA. Now consider a small error of da in the estimate of o. Let the corresponding error in/~ be d/~. We have d p = A - ~dB Aa
(29)
d B = diag (dflp),
(29)
where
and dflp = pfip a da. Therefore, 2p dp = A _ a dodBAzt ~< UA - 1 dffff][dBHail do = I;~La~l I[~11,
(30)
2p (3.3)
- -
where O is a constant, independent of o.
N ( + 1 P/o~(N/2)-1 da
(34)
N a ( I + ~ ) (NI2)-1 <
20
(35)
For a <
IIdp/da Na <~ I1~11
(36) 2p
F o r o <
Ud,uU/ll#N Ida/al
N ~<--. 2
(37)
Remark 4. It can be seen from equations (32) (34) above that the sensitivity of the eigenvalues of B with respect to the blur variance is independent the blur variance or. However, the sensitivity of p with respect to a varies as a polynomial in a, whose degree depends on N, the n u m b e r of basis functions at each scale. It is only possible to establish a conditional independence of the sensitivity of p with respect to a, although a stronger result would be desirable. Since the choice of N ties in our control, we believe that this conditional independence is not a serious deficiency. This shows that for o <
Deblurring Gaussian blur determine the blur parameter: In such cases, the best estimate for the blur parameter has to be adaptively estimated using some a prior~ knowledge of the scene, or by having a known test object in the scene, or by other heuristic means (like methods based on edge localization from zero-crossing contours etc.) and such estimates are prone to be noisy. In such a situation, it is imperative that the inversion be robust to perturbations in the blur parameter. The methods to tackle blur identification and estimation are, beyond the scope of this paper. For an analysis of the effect of noise on the magnitude of the coefficients, see Appendix 1. 4. E X T E N S I O N
TO SEPARABLEKERNELS
IN R a
In this section, we extend the method to separable, positive definite Gaussian blur kernels in higher dimensions. By a positive definite blur kernel in R a, we mean a function g: R d --, R, given by*
g(x) = e xTa~
5. D E B L U R R I N G W I T H A R B I T R A R Y B L U R KERNELS ON R a
In this section, we extend the formulation of the deblurring problem to non-separable kernels in two stages. First, we consider non-separable kernels defined on R. Then, using the tensor product, we extend the formulation to R d. Henceforth, by a slight abuse of terminology, we shall refer to deconvolution of any kernel as deblurring.
Non-separable kernels on R Our starting point is the same convolution equation (19). Following the general theme of the deblurring algorithm of Section 4, we obtain the multiscale Hermite expansions of f, g, and h as follows: N
f ( x ) = ~ It°.°H~(',po) n=O N n=0
R
x ~ f ( z l , z2 ..... zd)e(=' -"')/2"' dzl d'c2"" d'ca.
(40)
a(x)= T. ~°H.f',po) n=0 N
N
+ 12 fl~'H.(',a,)+...+ ~ fl;'H,(',a,) n=0
n=0
(41) N
h(x)= I2 ~.~°H.(,&) n=0 N
N
12 aP.'H.(.,p,) (42)
+ 12 ~,~'H.(',2,)+...+ n=O
n=O
Now, we continue in the same way as in the case of Gaussian deblurring in Section 3, and solve the convolution equation for the approximations off, 9, and h at the first scale. Substituting equations (43), (44) and (45) into equation (19), and setting ).j = pj + aj we get, N
tin/2] I
2., #ram: L m=o ,=o
x
Rd
R
n=0
N
h(x) = f ( x ) * g(x) = ~ f(z) g(x - z) d'c = I e(Xn-rn)2/2..... I e(x2-r2)2/2a2
N
+ 12 ,tlpl ~ . H. c ," pl)+.. + Y~ ~.H.(.,p.)
(38)
where G is a positive definite matrix of rank d. In the degenerate case where G is rank-deficient (i.e. some eigenvalues of G are zero), the blur occurs only in a p-dimensional space, where p = r(A). For remaining part of this section, we may assume without loss of generality that G is of maximal rank. We can obtain the eigenvalue decomposition of G as G = P A P T, where A is a diagonal matrix of the eigenvalues of G and P is an orthogonal matrix formed by the eigenvectors of G. Having determined the eigenvectors, we first eliminate the cross-terms in the quadratic form in the exponential by orienting it along its eigenvectors, which can be done easily by a suitable rotation of the axes. In what follows, x denotes the coordinate vector of the d-dimensional space variable whose coordinates have been measured with respect to the basis formed by the principal axes of G. In this new coordinate system, invoking Fubini's Theorem, the Gaussian blur in d-dimensions is given by
971
=
l~mt'20)
/~m-2t
- - - -
t!(m-- 2t)!
\
e
" l
]
( ~--o ln/2J(--1)k(2X/~)n-2k'~ . /~n!k=o 12 k!(n-2k)! J 2N [a/z] ( _ 1)*(2x/~)q-21 12 ~ q ] E ~. q=o l=o q!(q -- 20[ ,]
(43)
Equating coefficients of cop, we get the following system of equations:
R
(39) Thus, exploiting the separability of the Gaussian, the d-dimensional Gaussian blur can be expressed as d one-dimensionalGaussian blurs. Now, the deblurring method explained in Section 3 can be applied to each of these one-dimensional blurs. With minor modifications, we can easily see that the stability and the sensitivit3, analyses of Section 3 are also valid here.
(-- 1)("- P+~)/zm! (2~/a),-"~
(,,,
,,,,
= ( ~ o /gm(-- I)~-P)/2q'(2V/2)'~ p = 0, 1,...,2N. * In what follows,it is assumed that x is a vector in Rd,with appropriate changes for the other variables.
(44)
972
M. VAIRY and Y. V. VENKATESH
We have N + 1 unknowns in 2N equations. If there exists a unique solution to this system of equations, then N - 1 of these equations will be redundant and can be easily eliminated. Assuming for the moment that a unique inverse exists, we solve equation (47) to arrive at the first approximation of the solution /L Then, we approximate the error in the convolution equation by using a set of basis functions at a different scale, and solve the associated system of linear equations. This solution added to the solution obtained at the first stage yields the next improve estimate of the actual solution to equation (19). We continue in this m a n n e r until we have recovered f to the desired degree of accuracy. Now, we attempt to find the conditions under which an arbitrary convolution kernel is invertible. That is, we have to find out when this system of equations (one system for each scale) has a unique solution. To do this, we could form the associated matrix, as we did in Lemma 1, and find out when it is non-singular. However, the associated matrix turns out to be dense and the algebraic manipulations are very involved. Instead, we provide a simpler argument by borrowing a result from the theory of ordinary differential equations and arrive at necessary and sufficient conditions for the existence of a unique inverse. Let f ( x ) = Y ~ = o # k H k ( a , p ) SO t h a t t f(co)= Y~)'=OlZk(--j)kHk(pco, p). NOW, (19) becomes N
N
ttkHk(pco, p)(__j)k Z fltHt(aco, a)( _j)t = ~(co) k=O
1=0
(45) N
N
2 #kPk(co/4P)(--J) ke-2(p+a)(°2 Z fl~Ut(c°/x/~)(--J) t k--O
1=0
= h(co)
(46)
where Pk is the kth degree Hermite polynomial. The left band side of equation (51) can be written as a linear combination of the first N + 1 Hermite functions at scale p + a with modified coefficients #t N
conditions on f t of the form (D"f*)(Xo) = 0 for n = 0, 1. . . . . N. Hence the following theorem:
Theorem 3. Let SN(R ) be the space of all functions f defined on R and satisfying, (O~f)(xo)=O
c~
(Tf)(x) = h(x) = ~ f(z)g(x - z) dz.
/=0
(47)
We write equation (47) as
(S f*)(x) = h(x)
(50)
-oo
Then T is an isomorphism on SN(R). It is easily seen that the inversion, which we have shown to be well-defined, is also stable. In particular, we can b o u n d the n o r m of the error in h in terms of the n o r m of the error in f and the operator norm of S.
Non-separable kernels on R d The next step is to extend the method to tackle arbitrary blur kernels in R d. We simply generate wavelet bases in R d by tensor product. Before we proceed, we shall find it convenient to introduce some compact notation. In what follows, m, t, n, k, q, 1 are all multiindices in Zd+. To make the summations over multiindices meaningful, we assume an implicit total ordering (say, lexicographic) on these multi-indices. Let m = (mo, m l . . . . . md- O, [m] = ([mo], [m~] . . . . . [me 1], and ~ = o ( ' ) denote the summation of N + 1 terms (,) over the implicit (lexicographic) ordering on m. Further, let m! stand for mo !ml !"" me t !. The length of the multi-index m, Iml = mo + ml + ' " + ma- 1. A multi-index m raised to the power of another multiindex n stands for m~°m'~.... m~_-f. Let x e R a, ncZd+, a e R ~ . We generate the d-dimensional Hermite basis functions as follows:
j-d-a Hn(x, ff) =
~
Hn,(xj, aj).
(51)
j-O
#~Hk( p + a)(O, p + a)( _j)k ~, fltPt(aa), a)( --j)' =/~(co).
(49)
for some Xo ~R. Let g(x)= Z~, ~ f l , H , ( x ) , where H,(x) is the nth degree Hermite function and flo, ill,-.-, flu are constants, with fin ¢ 0. Let T be the operator given by
N
k
n=0,1 ..... N
We can write out multiscale expansion of f, g, and h similar to equations (40), (41) and (42) compactly as follows:
( ~ (48)
where the operator S=~,~_ofltD and f ( x ) = ~ = o # ~ H k ( x , p + tr). Since f * uniquely determines f, we only need to show that S is an isomorphism. S is obviously a homomorphism. To show that S is invertible, we have to establish that its kernel is the zero function. This can be obtained immediately from a set of initial
t The constant K in equation (7) is not important to the argument, and we may assume that it is absorbed in the definition of the Fourier transform of the Hermite functions.
[ m / 2 l t"- -
m=o #=m! Z,=o
l'm'2co ", ) ( NY / P')m - 2 t e-X;/2a] t.~m--2t)!
t;~] ( _ 1 ) k ( 2 ~ ) , -
2N =
/ 2k~
tq/el ( _ l),(2Xf2~)q-2,~
qE=o%Zq!,=o •
q!(q
21)!
e"
(52)
Equation (52) involves multinomials (multivariate polynomials), and leads to systems of linear equations again. In order to equate coefficients corresponding to coY; we assume that the multiscale expansion of f contributes co~and that ofg contributes cop ~. Choosing
Deblurring Gaussian blur = p + a and substituting in (19), we get, (~o
973
it appears that the Hermite functions provide a natural framework to analyse Gaussian deblurring and other related problems in signal and image processing. These functions have already been used successfully for problems in signal reconstruction using partial informationF )
p,,s~o ( - 1)1"- s)/2m' (2xfP)S']
( - 1)(" P+~)/Zm!(2~fa)P-s~
6. EXPERIMENTAL RESULTS
N
This section presents, briefly, the results of the proposed deblurring algorithm on a typical image (see Fig. 1). We demonstrate the effectiveness of the proposed deblurring algorithm on two images obtained
(_ 1)(q-p)/Zq[(2x/~lp~
(53) Naturally, solving this system of equations is more complicated than solving the corresponding system in the one-dimensional case. But apart from an increase in complexity by a factor of d, the associated matrix of equation (A3) is not upper triangular, and is in fact quite dense. In general, solving systems of linear equations has complexity O(N3). While the sparse nature of the matrix of equation (20) implies that solving equation (20) has O(N 2) complexity, the complexity of solving equation (A3) would be O(dN3). The extension to arbitrary kernels has several important applications in the fields of signal and image processing. For instance, a well known task in image processing is edge detection. Edges in an image are points of maximal gradient, and are obtained by setting the Laplacian of the image to zero. To validate the differentiation operation and reduce the effect of noise, the image is first convolved with a suitable Gaussian and the zero-crossings of the Laplacian of this new image (the so-called LoG zero-crossing contours) are points of maximum intensity variation in the original image. It turns out that the Laplacian of the Gaussian is just the second degree Hermite function H2(x ). Thus,
Fig. 1. Original image of GIRL.
Fig. 2(a). Blurred image: ~x = 1.0, ~ry= 1.0.
Fig. 2(b). Restored version of image shown in Fig. 2(a).
974
M. VAIRY and Y. V. VENKATESH
Fig. 3(a). Blurred image: ox = 2.0, ~ry= 2.0.
Fig. 3(b). Restored version of image shown in Fig. 3(a).
Fig. 4(a). Blurred image with noise: crx = 1.0, cry = 1.0. The noise variance is 1000/12 in both the x and the y directions,
Fig. 4(b). Restored version of the blurred and noisy image shown in Fig. 4(a).
by blurring (corresponding to two different values of ~r) the image of Fig. 1. In the implementation, the reconstruction procedure used the first twenty Hermite functions at twenty different scales. The original image in Fig. 1, blurred by an isotropic Gaussian of variance 1.0, is shown in Fig. 2(a). The reconstruction is shown in Fig. 2(b). The second case is a more severe blur with an isotropic Gaussian of variance 2.0, as can be seen in Fig. 3(a). The reconstruction is given in Fig. 3(b). It is observed that multiscale inversion using the proposed deblurring algorithm has restored most of the image sharpness.
The effect of noise is shown in Fig. 4(a), where the original image has been corrupted by uncorrelated, uniformly distributed noise with zero mean. The reconstruction is shown in Fig. 4(b). It can be seen that the reconstruction is stable with respect to noise. 7. C O N C L U S I O N S
The problem of deblurring the Gaussian blur is considered in the framework of a wavelet array transform for representing both the given (blurred) and the unknown images. A new method is proposed for de-
Deblurring Gaussian blur blurring, in a multi-scale framework, based on the use of generalized Hermite functions. The deblurring problem is reduced to solving sparse systems of simultaneous, linear equations, which can be done analytically, and is computationally efficient. The method directly extends to higher dimensions by making use of the separability of tile Gaussian kernel, and the deblurring operation on Euclidean d-space reduces to d independent deblurring operations on R. We have also formulated the problem for nonseparable kernels on R. While earlier attempts in the literature have either suffered from lack of stability or have established stable inversion procedures for restricted domains, we have established the stability of the inversion for perturbations in the blurred signal and also (conditionally)for errors in the choice of the blur parameter for a wide class of input images. The proposed procedure is robust with respect to noise, and seems to work reasonably well even in cases where the blur parameter in not known a priori but has to be estimated by other means. Experimental results for a natural image show that the method is superior to the results of the literature in terms of (i) robustness; and (ii) generality of the assumptions made on the input images. Apart from deblurring for the purpose of signal or image restoration, the m e t h o d can also be used for blur identification, where a known test signal f is fed to the blur mechanism (for which we can arrive at a multiscale Hermite expansion), and the result of the blurring is observed. Then, multiscale deconvolution can be used to estimate the blur kernel and its parameters. Future scope for research in this area includes extensions of d this method to the whole o f L 2 (R), and to non-stationary blurs. Spatially varying blurs are common in camera defocus blurs, where the blur kernel's variance varies with the distance of the object from the optic centre of the camera lens. Given the camera parameters, a deconvolution method for such spatially varying blurs would provide depth information directly, thereby dispensing with the correspondence problem in stereopsis.
975
6. D. Marr, Vision. Freeman, San Francisco, California (1982). 7. Y. V. Venkatesh, Hermite polynomials for signal reconstruction from zero-crossings--part h one dimensional signals, IEF Proc. I 139(6), 587-596 (1992). 8. S. M. Mallat, Multiresolution approximations and wavelet orthonormal bases of L2(R), Trans. Am. Math. Soc. 315(1), 69-87 (1989). 9. I. Daubechies, Ten Lectures on Wavelets. SIAM, Philadelphia (1992). 10. Y. Meyer, Ondelettes et Operateurs--Tomes I e t II. Hermann, Paris (1990). 11. Bruijn, Uncertainty principle in Fourier analysis, in Inequalities, O. Shisha, ed., pp. 57-71. Academic Press, New Yrok (1965). 12. Y. V. Venkatesh, K. Ramani and R. Nandini, Hermite sieve as a wavelet array for 1-D and 2-D signal decomposition, lEE Proc. I Vis. Image Signal Process. 141, 348-356 _(1994). 13. R. Nandini, Wavelet array transform for signal representation using generalized Hermite polynomials, M.S. Thesis, Department of Electrical Engineering, Indian Institute of Science (1992). 14. G. Szego, Orthogonal polynomials, Am. Math. Soc. (1975). 15. N.N. Lebedev, Special Functions and their Applications. Prentice-Hall, Englewood Cliffs, New Jersy (1965). 16. J. R. Higgins, Completeness and Basis Properties of Sets of Special Functions, Cambridge Tracts in Mathematics. Cambridge University Press, Cambridge (1977). 17. Madhu Vairy and Y, V. Venkatesh, Wavelet array representation for deblurring Gaussian blur, Technical Reports Department of Electrical Engineering, Indian Institute of Science, Bangalore, India (1993). 18. Gnedenko, Theory of Probability. MIR, Moscow (1975). APPENDIX
Noise Analysis
We briefly analyse the robustness of the method to noise. Let q(x) be (broad sense) stationary, uncorrelated additive noise with zero mean and variance a 2. At a fixed scale, let 6#k (a random variable)* be the change in the coefficient #k due the additive noise r/(x). go
6/~k= S q(T)Hk(T)dv"
(A1)
~c
If E denotes the expectation operator, then assuming interchangeability of.the integral and the expectation, we get, oD
E[6#k] = Acknowledgements--MV would like to thank Shanmukh, Ganesh Murthy and Chris Heil for several stimulating digcussions. REFERENCES
I. A. N. Tikhonov and V. Y. Arsenin, Solutions of ill-posed problems, F. John, ed., Wiley, New York (1977). 2. F. John, Numerical solution of the equation of heat conduction at preceding times, Ann. Math. 17, 129-142 (1955). 3. R. A. Hummel, B. B. Kimia and S. W. Zucker, Deblurring Gaussian blur, Comput. Vision Graphics Image Process. 38, 66-80 (1987). 4. B. B. Kimia and S. W. Zucker, Analytic inverse of the discrete Gaussian blur, Opt. Engng. 32(1), 166-176 (1993). 5. J. Biemond, R. L. Lagendijk and R. M. Mersereau, Iterative methods for image deblurring, Proc. IEEE 78(5), 856-890 (1990).
1 _ ~ E[q(z)JHk(T)dT=O. IIHkl]~ -
(A2)
Similarly, E[6/t~] = ~
1
~
cc
_~ -~f E[q('Orl('c')JHk('c)Hk('r')d'cd'V
O-2
= iii/4k 1~.
(h3)
Therefore, the variance of the kth coefficient at any scale is decreased by a factor of IIHkII2, where IIHII2 = 2k IIHk- a II2 for k = 1, 2, 3..... and Ilno]]2 = x//~.
* Under certain conditions, the integral in equation (38) converges to a random variable, according to Khinchin's theorem. See Ref. [18] for the proof of Khinchin's theorem.
976
M. VAIRY and Y. V. V E N K A T E S H
About the Author--Y. V. V E N K A T E S H got his Ph.D. from the Indian Institute of Science for a dissertation on the stability analysis of feedback systems. He was an Alexander von H u m b o l d t Fellow at the Universities of Karlsruhe, Freiburg and Erlangen, Germany; National Research Council Fellow (U.S.A.) at the Goddard Space Flight Center, Greenbelt, Maryland; Visiting Fellow at the Australian National University, to name a few. His research m o n o g r a p h on stability and instability analysis of linear-nonlinear time-varying feedback systems has appeared in the Springer Verlag Physics Lecture Notes Series. His present work is on signal representation using wavelet-link arrays, and reconstruction from partial information (like zerocrossings). He is a Professor at the Indian Institute of Science, Bangalore, India, and currently the Chairman of the Department of Electrical Engineering. He is a Fellow of the Indian National Science Academy, Indian Academy of Sciences and the Indian Academy of Engineering.
About the Author--M. VAIRY studied Electrical and Electronic Engineering at the Birla Institute of Technology, Pilani, Rajasthan, India. At the time of submission of the paper, he was a Research Scholar (M.Sc.) at the Indian Institute of Science, Bangalore, India. Now he is a Research Assistant at the Courant Institute, New York University, New York, U.S.A.