IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-4, NO. 2, MARCH 1982
Nonlinear Restoration of Noisy Images
JEAN-FRAN§SOIS ABRAMATIC, MEMBER, IEEE, AND LEONARD M. SILVERMAN, FELLOW, IEEE
Abstract-The restoration of images degraded by an additive whitenoise is performed by nonlinearly filtering a noisy image. The standardWiener approach to this problem is modified to take into account theedge information of the image. Various filters of increasing complexityare derived. Experimental results are shown and compared to thestandard Wiener filter results and other earlier attempts involving non-stationary filters.
Index Terms-Image restoration, nonlinear filter, nonlinear restora-tion.
I. INTRODUCTIONIMAGE restoration is often defined as the process of recover-
ing an original image from a distorted version. Wheneverexact restoration is not feasible (this happens, for example,when a random noise is involved in the distorting process),the restoration problem becomes an approximation problem.One has to define a distance between two images and try tominimize the distance between the original image and itsrestored version. Our paper is related to the case where thedistorted image is the sum of the original image and a sta-tionary white noise process. By extending the well-knownWiener and Kalman filtering techniques to the case of 2-Dsignals, Helstrom [10], Pratt [16], Nahi [13], Habibi [9],Attasi [4], Aboutalib et al. [1], Woods and Radewan [18],and Hunt [3] have proposed various algorithms to performthe restoration. These algorithms have in common two basicassumptions. First, the image is considered to be the realiza-tion of a 2-D stationary random process. Second, the dis-tance between two images is an L2-distance, generally a meansquare error between two random processes. As a conse-quence of these assumptions, the restored image is obtainedas the output of a 2-D linear filter whose input is the image.The algorithms differ in the type of filters chosen (recursiveor nonrecursive), but all of them are essentially low-passfilters. This is due to the fact that the noise is supposed to bewider band than the signal. They are also stationary filters,due to the stationarity assumption for the image model. Un-fortunately, these assumptions and their consequences pro-duce unsatisfactory results when dealing with images. Edgeswhich are often of great interest are smeared by the low-passfilters. As a matter of fact, edges can be viewed as high-frequency information of the original image that is not taken
Manuscript received April 3, 1980; revised May 8, 1981. This workwas supported in part by the National Science Foundation under GrantENG77-1883 1.
J.-F. Abramatic is with the Institut National de Recherche en Infor-matique et Automatique, Le Chesnay, France.L. M. Silverman is with the Department of Electrical Engineering,
University of Southern California, Los Angeles, CA 90007.
into account in the model. Another way to look at it is tointerpret edges as locations in the image where the modelchanges, contrary to the stationarity assumption. Regardingthese issues, more recent work in restoration of noisy imageshas tried to derive nonstationary and/or nonlinear algorithms.Nahi and Habibi [14] first dealt with the case of two models.Ingle et al. [11] allowed the parameters describing the imagemodel to change inside the image, deriving then an identifica-tion-estimation algorithm. Ingle and Woods [12] consideredthe case of five models (describing edges of different orienta-tions or isotropic situation), deriving then a multiple-channelestimation algorithm. Using a totally different approach,Anderson and Netravali [2], then Peyrovian and Sawchuk[15], Barba [6], and Clara [7] derived nonstationary non-recursive filters that give a compromise between the loss ofresolution around the edges, and the effect of the noise in theflat regions. Most of the latter papers used an approach firstintroduced by Backus and Gilbert [5] for the restoration ofgeophysical data. Following an approach first suggested byFrieden [8], and generalizing it, we will try to fill the gapbetween these new attempts and the more classical Wienertechnique.We propose a general formulation of the problem as a non-
linear problem. We show how the "Backus and Gilbert" ap-proach is a particular case of it and how it can be viewed as anonlinear generalization of the parametric Wiener filter intro-duced by Hunt [3]. We propose another simplified modelthat leads to an easier design for the nonstationary filter.We show also how a nonisotropic nonlinear filter can producebetter results due to the nonisotropy nature of an edge. Wepropose finally suboptimal separable filters that are very easyto design and give satisfactory results.
II. THE PROBLEM
Stationary Case- Wiener ApproachThe Wiener formulation of the problem can be posed in the
following way. We try to find an estimate z(i, j) of a signalz(i, j) knowing the measurements y(i, j), where
y(i, j) = z(i, j) + w(i, j) (1)
and w(i, j) is assumed to be a white noise with known variancea2; Z(i, j) is taken to be the realization of a stationary randomprocess with known covariance A(k, 1):
IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-4, NO. 2, MARCH 1982
z(i, j) and w(i, j) are assumed to be independent processes.shall look for a smoothing linear estimate, that is,
z(i,j)=Z h(k,l)y(i - k,j- 1)k, I
where k and 1 can take positive and negative finite valThe estimate will be said to be optimal if it minimizes
= E [z(i, j) - Z(i,j)]2.
The error function ew can also be written as
e = a2 E h2(k, 1) + L 3 (6(k, 1) - h(k, 1))k,l k,l k,l
*A(k - k', 1 - 1') (5(kII 1') - h(k', 1')).The first term of the right-hand side of expression (6) wilcalled the noise term, the second will be called the resolu,term. We see that the Wiener formulation weights tlequally. The solution is then given by
H(zI, z2) = h(k,l) Z-kz = S(Z1,Z2)k,lISZ15Z)+a
Our approach to the restoration of noisy images will corin modifying the model of (1) and the criterion of (6'adjust them to the local characteristics of the image. Thisbe done by taking into account the masking effect whichdescribe in the text.
Earlier attempts to deal with the tradeoff between nreduction and resolution are described in [3]. They appl:the deconvolution of noisy blurred images. In our case treduce to compromising between the Wiener filter andidentity filter. The parametric Wiener filter is defined as
H(z1,z2) = SZ(Z1, Z2)SH(z= , Z2) + -yu'
One of our filters will be of the same form, but with y a fttion of image location, generalizing this filter in a nonlirway. The choice of y is determined by the masking funct:
Masking EffectThe masking effect derives from the fact that noise is
ceived by a human eye differently depending upon whethcoccurs close to or far from an edge. The eye is quite sensito a small amount of noise in a flat field, but is able to tcate a large amount of noise in the surroundings of an esThis effect was first taken into account in the restoratiornoisy images by Anderson and Netravali [2]. As a quanitive measure of spatial detail (or edge content), they cho"masking function" M(*, *) of the form
M(i,j) =E C MH(-k2+ -12(Im(k, 1)lk, I
+ m (k,I)|)where mH and mv are the horizontal and vertical gradi4of the original image. They defined then an equivalent w:noise
v(i, j) = f(M(i,j)) w(i, j)
by subjective tests, where they presented first, images withnoise w added to the points where the masking function takesthe value M, and second, images where the noise v(i,j) hasbeen added to the whole picture. Observers then decidedwhen the noise effects were equally disturbing. As a resultof this procedure, the visibility function f(-) was completelydetermined from the original image.Later on, Barba [6] and Clara [7] in the case of color images
used a simpler visibility function f(-), taking advantage of aremark in [2] arguing that a coarsely tuned decreasing func-tion will be a good guess, and will produce satisfactory experi-mental results.
Backus-Gilbert ApproachOnce the visibility function is determined, all of the authors
previously mentioned used the same procedure. At each point(i, j), the coefficients of the nonstationary filter h"i(k, 1) thatgives the restored image
(7) z(i,j) = E h1j(k, 1) y(i - k, j - 1)k, I
(11)
are defined by solving a nonlinear optimization problem de-rived from the Backus-Gilbert [5] approach to the deconvolu-tion of geophysical data.
* Minimize
J(h"(k, 1)) = oa'i L hj(k, 1)2k, I
+ ( - a(i)E hj(k, 1)2 (k2 + 12 )k, I
under the constraint
(8) E3 hi(k, 1) = 1.k, I
anc- * The coefficients cx4' are then adjusted to satisfytear , h"(k, 1)2 . f)'(M(i, j)) = 0ion. k,I
(12)
(13)
(14)
where y, 0 are tuning parameters.
per- This method and some variants trying to take into accounter it the nonisotropy of the problem gave good experimentalitive results, but one can see that the process of defining the proper
)ler- filter at each point is not easy. Moreover, nothing was sug-dge gested as to estimating the masking function from a noisyn of image. The results were then either academic (in the case
tita- where the masking function was calculated from the originalse a image, which is not supposed to be available), or much poorer
(in the case where the masking function was calculated di-rectly on the noisy image). Finally, one has to point out thatthe filter will not depend on the signal to noise ratio whichseems not to be a good feature.
(9) Our approach will try to make the design of the nonsta-tionary filter much easier. We shall also take into account the
ents nonisotropy of the problem in our basic formulation. We shallhite propose a method to estimate the masking function from the
noisy image, making the whole problem more realistic. Finally,(10) our filters will depend upon the signal-to-noise ratio.
142
ABRAMATIC AND SILVERMAN: NONLINEAR RESTORATION OF NOISY IMAGES
III. THE NONLINEAR PROBLEM FORMULATIONAs already indicated, the approach to be taken here is a
modification of the Wiener equations to take into account themasking effect. First, we modify the measurement of (1) as
y(i,j) = V(i,j) + w*(i,j) (15)
Rij(h'i(k, 1), v(-, *), Av(- - ))= R W(Oi(k, 1), v(-, -), Av(, -))- E E (6(k, 1) - h"i(k, 1)) Av(k - k', 1 - 1')
k,l k,l
* (5(k', I') - h'j(k', 1')).where v(i,j) will be termed the visible signal and w*(i, j) theequivalent noise, with
v(i, j) = z(i, j) + (w(i, j) - w*(i, )) (16)
w*(i,j) = W11(z(, *), w(i, j)). (17)
Equation (16) defines v(i, j) as a noisy version of z(i, j) whichcannot be distinguished from z(i, j) because of the maskingeffect. The masking effect is measured through Wij, whichlocally depends upon the image and noise. Second, we modifythe criterion (6) as
eij = Ni (h 1(k, 1), z( .
+ Rij(h"j(k, 1), z(, Az *)) (18)
where Nij and Rij are the noise and resolution terms, respec-tively. The problem can then be formulated in the followingway.
* Given the measurement y(i, j)
y(i,j) = v(i,j) + Wi1(z(I, *), W(i,j)) (19)* find a linear estimate u(i,j) of v(i,j)
Depending upon the choices made for Wij, N0i, and Rij, weshall show more precisely how the masking effect can be takeninto account. In all cases we shall choose a suboptimal ap-proach which consists in first estimating the values taken byWij, N\i, and Ri,, and then in a second step, solve the resultinglinear nonstationary problem. For clarity, we will present thesteps in the reverse order. We will assume that the estimationsof Wj, N0j, and Rij have been accomplished, and then presentthe solution of the nonstationary process.
Generalized Backus-Gilbert ApproachUsing the above formulation one can fill the gap between
the Wiener and Backus-Gilbert approach by choosing
But in that case, z(i, j) = v(i, j). Therefore, compared to theWiener approach the noise term in the criterion is the onlyterm to be changed, scaled by the visibility function. M(i, j)is defined by (9), and f(-) is a monotonic decreasing functionfrom 1 to 0 as M goes from 0 to oo. The optimal solution isthen given by
where fi0 denotes f(M(i, j)).This approach does not use the concept of visible signal, but
just weights in a nonstationary fashion the noise and resolu-tion terms in the criteria. We suggest now a "signal equiva-lent" approach which leaves the criterion unchanged andmodifies the measurement equation.
The term fii need not be the same in the two approaches;however, it plays the same role introducing the adaptivity,and we shall keep the same notation.In this case the optimal solution minimizes the criterion
(21). However, this criterion is just the standard Wienercriterion applied to the signal v(i,j) instead of z(i,j). There-fore, the solution will be given by
k Z2= E S y(z ,Z2)H0I(zl,)ZZZ2)>=3 hUl(k,)zzjI _____ 2___k,l1 SZ(z 1,z2) + uT
Hij(zI Z2) = fii S (Z Z2)2 + (1- fi).SJZI, Z2) + cJ
(31)
(32)
(33)
Comparison Between the Two ApproachesIn both cases, as fi0 varies from 0 to 1, the filter chosen is a
compromise between the identity filter and the Wiener filter.When fif is equal to zero, the equivalent noise is zero, and themeasurement itself is considered to be equivalent to the image
143
111(h','(-, -) v(-, -), Av(-, -)).
IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-4, NO. 2, MARCH 1982
with respect to the human eye. When fii is equal to 1, theequivalent noise is equal to the real noise, the visible signal isthe signal itself, and the Wiener filter is implemented. How-ever, the way the tradeoff is made in the case where fii is be-tween 0 and 1 is easier to handle in the second approach. Therestored image is formed by taking a linear combination be-tween the standard Wiener filter output, and the noisy imageitself. The coefficients of the linear combination are given bythe visibility function. In the generalized Backus-Gilbertapproach, one has to implement a different filter for eachdifferent value of the visibility function. The "signal equiva-lent" approach appears to produce a very straightforwardmethod for adaptively restoring a noisy image.To enhance the procedure, we can take advantage of the
edge information contained in the masking function, and use anonisotropic filter. This will be taken into account by modify-ing the resolution term in the criterion (21).
IV. ESTIMATION OF THE MASKING FUNCTIONAs pointed out earlier, starting from our nonlinear formula-
tion of the restoration problem, we have assumed so far thatthe local characteristics of the image have been estimatedpreviously, and we have used these estimates in our designprocedures. We now propose a method for estimating themasking function. Thinking of the nonisotropic extensions ofour procedures, we will generalize the concept of a maskingfunction. We will define
MD(i, j) = Z exp [- (i - k)2 + (j - l)2/do] ImD(k, 1)|k, 1
(34)where MD will be calied the masking function in direction Dand mD(k, 1) is the gradient of z in direction D at point (k, 1).Dealing with discrete sets, D will take practically a finite num-ber of values, basically horizontal and vertical, and possiblyfirst and second diagonal orientations. Our estimation prob-lem can then be stated as follows.Given the measurement
y(i,j) = z(i,j) + w(i,j) (35)find an estimate
fD(i j) of MD(i, j).
In practice, we shall be interested in the four followinggradients:
mH(k, l) = z(k, l + 1) - z(k, l - 1)
mV(k, 1) = z(k + 1,1) - z(k - 1,1)
mDl(k,l)=z(k+ 1,1+ l)-x(k- 1,1-1)
mD2(k, 1) = z(k + 1, 1- 1) - z(k - 1,1 + 1). (36)As we do not want to look at these estimations, but just usethem as a guide for the restoration process, we do not need avery accurate estimation. This is the reason why the standardWiener procedure will be chosen for these particular estima-tion problems. This will be a suboptimal procedure, since we
use linear filtering where nonlinearity is introduced by themagnitude in (34). More precisely, the estimation of the mD7sis the linear problem solved by means of the Wiener filter. Inthe horizontal case, for example, the masking function will beestimated in the following way.Given the measurements ,uH(k, 1)
,uH(k, l) =y(k, I + 1) - y(k, 1- 1)
find an estimate
mH(k 1) of mH(k, 1) (37)
where mH(k, 1) is defined by (36).From (35) and (37) we can rewrite
pH(k, 1) = mH(k, 1) + v(k, 1)
and
v(k,l)=w(k+ 1,1)- w(k- 1,1).
(38)
(39)This problem can be seen as the estimation of a signal dis-torted by a colored noise. The estimator will then be given by
m'(k, I) = , af p, q) ,utH(k - p, I - q)p, q
where
A(z1, Z2) = E a(p, q) z-Pzjqp, q
S'm(Zl , Z2)Sm(Zl, Z2) + Sv(Z1, Z2)
(40)
(41)
Sm and S are the power spectra of the signals m(k, 1) and(k, 1), given by
Sm(ZI, Z2) = (Z2 - Z21) SZ(zl, Z2)
SU,(Z1, Z2) = (Z2 - Z21) c42 .
(42)
(43)It is clear that A(-, -) will be the standard Wiener filter whichwill be of great interest in the design of the "signal equivalent"filter. We shall then get an estimation of the masking functionMH(k, 1) through
(44)It is clear that similar results can be realized for masking func-tions associated with any orientation. We want to emphasizeat this point that the implementation of the "signal equiva-lent" filter, including the estimation of the masking function,is quite appealing. Basically, the output of the same standardWiener filter will be used for estimating both the maskingfunction and for calculating the final estimate. This scheme isdescribed in Fig. 1.
V. NONISOTROPIC FILTERING
Any of the approaches can be made nonisotropic by choos-ing a resolution term different from Ri), which is the general-ization to the nonlinear case of the Wiener resolution term.We chose the form
144
ABRAMATIC AND SILVERMAN: NONLINEAR RESTORATION OF NOISY IMAGES
=E (8(k, I) - 0i(k, I)) Atz(k - k' I - I')k, I k, I
* (5 (k', I') - h'j(k', I')) rij (z( ), k- k',I - 1'). (45)This can be sought as weighting, unequally, the measure-
ments available in the surroundings of the point (i, j).Our weighting function rij will then be expressed as
rij(z(-, *), k - k', 1 - 1') = pij(Mk4;), k - k', I - 1'). (46)
This modification of the criterion can be interpreted as amodification of the autocorrelation function A(k, 1) into anonstationary "autocorrelation function" Aij(k, 1)
A1j(k, 1) = A(k, 1) * pi(Mi(;), k, 1).
Defining
Si"(Z1, Z2) =': A,1(k, 1) Z-kZ-lk, I
the optimum filter is then given by
)~ Sz'(z 1 Z2)H" (zI, Z2)
(47)
Separable Nonisotropic FiltersIn the case where the autocorrelation function is separable
such that
A(k, 1) = h(k) * A (l) (51)the Wiener filter is not separable in general. A suboptimalfilter will be the separable filter built with the two 1 -D Wienerfilters
H1(Z , Z2) =Hi/(z) * Hi'i(z2). (52)We shall therefore obtain for the Backus-Gilbert approach
H5?(z) = Sij(Z)andforth sigl euiv,al
and for the "signal equivalent" approach
(53)
(54)(48) H1f(z)=f4' S (2) + -
where
fD =-f(MD(i, j))(49)
in the case of the Backus-Gilbert approach. In the "signalequivalent" case we find
2 (1 -ii,Sz'(z1, Z2) +Or (50)
At this point we want to emphasize that although we callMD a masking function, as it is a direct extension of the oneintroduced before, we do not use it as a quantitative measureof the masking effect, which is related to the noise visibility.We only take advantage of the edge information available toprevent the filtering process from blurring the edges.We want to point out also that the design of the filters
remains easier than the previous procedures suggested inAnderson and Netravali [21, Barba [6], and Clara [7]. No opti-mization procedure is involved in calculating the coefficientsof the filters as the results are given by analytic expressions.However, as the filter has to be calculated at each point, it ispossible that good suboptimal procedures might be found toreduce the amount of computation while still giving goodrestoration results. We now present two such suboptimal pro-cedures based upon two different assumptions on the auto-correlation function A(k, 1).
(55)and f is again a monotonic decreasing function from 1 to 0 asM(i, j) goes from 0 to oo. This design procedure leads to verysimple calculations. It is useful to consider the very oftenused case where the autocorrelation function is supposed to befirst order, that is,
1 - p2S(Z)-= (z- p)(z1 -p) (6
where p is the correlation coefficient. We can then introducethe nonisotropy only through f4 and f1V In this case thedesign reduces to solving two second degree polynomial equa-tions per pixel (i, j), which is a rather light burden.
Four-Direction Nonisotropic FiltersIn the case where the autocorrelation function is supposed
to be circular symmetric instead of separable, it may be usefulto introduce directions other than horizontal and vertical.Basically, one must calculate the 1-D filters associated witheither one of our two approaches, that is through (53) or (54).Depending upon the number of directions for which 4'J has
been calculated, one can then reconstruct the 2-D filter bybilinear interpolation, for example. Although the use of Ndirections is theoretically feasible, it seems reasonable from apractical point of view to limit ourselves to the use of four
145
(56)
IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-4, NO. 2, MARCH 1982
w(i,j)
z(i,i) A Y(i,i) Direction #1 _ _ Masking #1
IDirection # _.w7recaton#2mD(ai,ni) ExpotmDi,aj)
Estimation Maqnitude Exponential
Imoot I ng
VisibilityFunction
H_(z1 z(i2j)
Fig. 2. Implementation of the nonisotropic filters for four directions.
basic directions, adding to the horizontal and vertical direc-tions the orientations of 450 and 1350. In this case anothersuboptimal approach, suggested by Clara [7], is to use eithera separable or a separable rotated of 450 2-D filter, dependingupon the respective values of the visibility function.An implementation of a nonisotropic nonstationary linear
filter is suggested in Fig. 2.
VI. EXPERIMENTAL RESULTS
Several experiments have been performed that illustrate ourapproach to image restoration. All of them assume the imageto be a realization of a 2-D random process with separablecovariance
Az(k, 1) = AH(k) . Av(l) (57)where AH and AV are covariances of Markov random pro-cesses of order 1. The correlation coefficient p is taken to0.9 or 0.96. The filters h11(k, 1) we consider are separable,as previously discussed, and their impulse response is truncatedto various sizes. They are scaled to satisfy
L3 h1(k, 1) = 1. (58)k, I
The emphasis has been made on the comparison of the "gen-eralized Backus-Gilbert" and "signal equivalent" approaches,and on the effect of calculating the masking function, eitheron the original image, or on the noisy image. This makespossible a comparison of our procedure with earlier algo-rithms. Following a suggestion of Anderson and Netravali[2], we have chosen as a visibility function a simple decreas-ing function
f(M(ij)) = ____59__aM(i,j) + 1
that decreases from 0 to 1 as the masking function goes from0 to oo. The term a is a tuning parameter which controls thetradeoff between the Wiener filter and the identity filter forthe entire image. If a is equal to 0, we use the stationaryWiener filter, whereas a large a gives the identity filter.
ResultsFig. 3 illustrates the variations of the impulse response of the
1-D filter for different values of the masking function and
fixed parameter a. Results of the algorithms are then shownin Figs. 4, 5, 6, and 7. The original (top left), noisy (topcenter) and Wiener restored (top right) images are presentedin the first row of each figure, for various signal to noise ratioand correlation coefficient p. In Figs. 4 and 5 the second rowillustrates the results of the "signal equivalent" approach, andthe third row the results of the "generalized Backus-Gilbert"approach, for a masking function computed on the originalimage (left), the noisy observation (center), or the Wienerrestored image (right). Figs. 6 and 7 present results of bothapproaches, for a masking function evaluated from the Wienerrestored image.
Finally, Fig. 8 displays the masking functions of the original,noisy, and Wiener restored images, and Fig. 9 illustrates theeffect of the tuning parameter a. The left bottom picture isobtained for a close to zero, whereas the bottom right picturecorresponds to a large value of a. In all cases the impulseresponse of the filters have been truncated to 9 pixels. Similarexperiments have been performed on wider impulse responses,but the resulting quality was not sensitively improved. Also,in all images except Fig. 9, the parameter a was tuned sub-jectively to its "optimal" value to give the best possiblequality.
VII. POSSIBLE EXTENSIONS OF THE NONLINEARRESTORATION PROCESS
The technique we have presented can be extended in variousways to more complicated restoration problems. Our ap-proach to the restoration problem by means of a generalizedWiener formulation makes it applicable to the deconvolutionproblem. In this case, the distortion is supposed to have twocauses. First, the signal is blurred. Second, some noise isadded to the blurred image. A Wiener stationary solution isavailable for this problem for the case where the blur processis known. If B((zi, z2)) denotes the Z-transform of the -im-pulse response of the blur, the Wiener solution is given by
(60)HB*(zS,Z(Zl, Z2)HW(Z1~Z2 B(zi, z2)|z(, Z2) + AWThe signal equivalent approach will give the following
ABRAMATIC AND SILVERMAN: NONLINEAR RESTORATION OF NOISY IMAGES
0
0
A
M = .5: Edge of High Constrast
M
N
= .05 Edge of Low Contrast
= 0: No Edge
-4 -3 -_z - a
Fig. 3. 1-D filter impulse response.
Fig. 4. Results of both approaches for various masking functions.SNR = 5 (variance), p = 0.9.
The use of the inverse filter, even in a nonstationary way,does not seem to be a good feature. On the other hand, thegeneralized Backus-Gilbert approach will give
H(zlSZ2)= B*(ZZ 2) Sz(Zl Z2) . (62)
The inverse is only the limit of the filter when fij goes to 0.This approach appears at first glance more appealing for thedeconvolution case. Further studies on this point are underinvestigation.Another possible extension of our work can be the use of
recursive filters to perform restoration. In that case the signalequivalent approach appears more promising, at least for thenoisy case. The nonstationary filter will simply be the nonsta-
ri g. Z. Results o0 both approaches for vanous masking fUtnctions.SNR = 2 (variance), p = 0.9.
tionary linear combination of the Kalman filter and theidentity filter. For the Backus-Gilbert approach, the filterchanges at each point. Problems of initial conditions investi-gated in Rajala [17] will have to be solved here too.
Finally, one can think of applying our technique to therestoration of color images. Clara [7] showed that using anadequate trichromatic representation for the color signal, themasking effect for each of the three components is indepen-dent of the other two, simplifying then the extension of theprocedure developed for black and white images.
VIII. CONCLUSIONWe want to emphasize that our new approach to the prob-
lem of restoring noisy images has allowed us to derive a varietyof algorithms. The range goes from the basic "signal equiva-
147
* 8
IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-4, NO. 2, MARCH 1982
Fig. 6. Results of both approaches. SNR = 5, p = 0.96.
Fig. 8. )m right,Wiener bottom left.
lent" filter that can be implemented with "real-time" hard-ware, to more sophisticated nonisotropic filters that willrequire a software implementation.
REFERENCES[1] A. 0. Aboutalib, M. S. Murphy, and L. M. Silverman, "Digital
restoration of images degraded by general motion blur," IEEETrans. Automat. Contr., vol. AC-22, no. 3, pp. 294-302, 1977.
[2] G. L. Anderson and A. N. Netravali, "Image restoration based onsubjective criterion," IEEE Trans. Syst., Man, Cybern., vol.SMC-6, pp. 845-853, Dec. 1976.
[3] H. C. Andrews and B. R. Hunt, Digital Image Restoration,Prentice-Hall Signal Processing Series. Englewood Cliffs, NJ:Prentice-Hall, 1977.
[4] S. Attasi, "Modelling and recursive estimation for double in-dexed sequences," in System Identification. Advances and CaseStudies, R. K. Mehra and D. G. Laniotis, Eds. New York:Academic, 1976.
[5] G. Backus and F. Gilbert, "Uniqueness in the inversion of inac-curate gross earth data," Phil. Trans. Roy. Soc. London, vol. 266,pp. 123-192, 1970.
[6] B. D. Traitement, "Lineaire d'images degradees par filtrageadaptatif avec critere psychovisuel de qualite," presented atColloque Nat. sur le Traitement du Signal et Appl., Nice, France,1977.
[7] F. J. Clara, "Filtrage adaptatif d'images couleur avec criterepsychovisuel," Docteur Ing. dissertation, Paris VI, 1980.
Fig. 9. Effect of the tuning parameter a.
[8] B. R. Frieden, "Image enhancement and restoration," in Topicsin Applied Physics, T. S. Huang, Ed. New York: Springer, 1975.
[9] A. Habibi, "Two-dimensional Bayesian estimates of images,"Proc. IEEE, vol. 60, pp. 878-883, July 1972.
[10] C. W. Helstrom, "Image restoration by the method of leastsquares," J. Opt. Soc. Amer., vol. 57, no. 3, pp. 297-303, 1967.
[11] V. K. Ingle, A. Radpour, J. W. Woods, and H. Kaufman, "Recur-sive estimation with non-homogeneous image models," in Proc.IEEE Pattern Recog. Image Processing Conf., Chicago, IL, 1978,pp. 105-108.
[12] V. K. Ingle and J. W. Woods, "Multiple model recursive estima-tion of images," in Proc. IEEE Int. Conf. Acoust., Speech, SignalProcessing, Washington, DC, 1976, pp. 642-645.
[131 N. E. Nahi, "Role of recursive estimation in statistical image en-hancement," Proc. IEEE, vol. 60, pp. 872-877, July 1972.
[14] N. E. Nahi and A. Habibi, "Decision-directed image enhance-ment," IEEE Trans. Circuits Syst., vol. CAS-22, pp. 286-293,1975.
[15] M. J. Peyrovian and A. A. Sawchuk, "Restoration of noisy blurredimages by a smoothing spline function," Appl. Opt., vol. 16,pp. 3147-3153, 1977.
[16] W. K. Pratt, "Generalized Wiener filtering computation tech-niques," IEEE Trans. Comput., vol. C-21, pp. 636-641, July1972.
[17] S. A. Rajala, "Adaptive non-linear image restoration by a modi-fied Kalman filtering approach," Ph.D. dissertation, Rice Univ.,1979.
[18] J. W. Woods and C. H. Radewan, "Kalman filtering in two dimen-
148
IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-4, NO. 2, MARCH 1982
1.. X1 ~..Jean-Frangois Abranatic (M'78) graduatedfrom the Ecole des Mines de Nancy in 1971.He received the Docteur-Ing. degree from theUniversity Paris XI, Orsay, in 1975, and theDocteur es Sciences degree from the UniversityParis VI in 1980.Currently, he is Ingenieur de Recherche at
INRIA (Institut National de Recherche enInformatique et Automatique), Le Chesnay,France. where he has been working in theImage Processing Group since 1974. During the
year 1979 he was on sabbatical leave at the Image Processing Instituteof the University of Southern California, Los Angeles. His currentinterests are signal processing, image processing, and analysis.
Leonard M. Silverman (S'60-M'66-SM'77-...........F 79) received the A.B., B.S., M.S., and Ph.D.
degrees in electrical engineering, all from Col-umbia University, New York, NY, in 1961,1962, 1963, and 1966, respectively.From June 1966 to June 1968 he was an
Acting Assistant Professor in the Departmentj of Electrical Engineering and Computer Science
A;|s at the University of California, Berkeley. Sincethat time he has been at the University ofSouthern California, Los Angeles, where he is
currently a Professor in the Department of Electrical Engineering. Hisresearch interests are generally in multivariable system theory, withapplication to control, estimation, image processing, and circuit theory.Dr. Silverman has been Associate Editor for Linear Systems of the
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, Vice President of the IEEEControl Systems Society, and Chairman of the IDC. He was the recipi-ent of 1968 Prize Paper Award of the IEEE Circuit Theory Group.
Digital Straight Lines and Convexity ofDigital Regions
CHUL E. KIM AND AZRIEL ROSENFELD, FELLOW, IEEE
Abstract-It is shown that a digital region is convex if and only ifevery pair of points in the region is connected by a digital straight linesegment contained in the region. The midpoint property is shown to bea necessary but not a sufficient condition for the convexity of digitalregions. However, it is shown that a digital region is convex if and onlyif it has the median-point property.
Index Terms-Convexity, digital geometry, straightness.
I. INTRODUCTIONC ONVEXITY is an important property of planar figures. It
is well defined and understood for the case of continuousplanar figures. However, when a planar figure is convertedinto a digital region for processing by a computer, it is nolonger clear what is meant by a digital region being convex.Many different definitions of convexity of digital regions havebeen suggested and studied [3], [5], [6], [10]-[12]. Re-cently, a new definition was proposed in [4]. A consequenceof the definition was that a regular digital region is convex ifand only if there exists a convex planar figure whose digitiza-tion is the given digital region. Thus, the definition turns out
Manuscript received May 30, 1980. This work was supported in partby the National Science Foundation under Grants MCS-7822313 andMCS-7623763.C. E. Kim was with the Department of Computer Science, University
of Maryland, College Park, MD 20742. He is now with the Departmentof Computer Science, Washington State University, Pullman,WA 99163.A. Rosenfeld is with the Computer Vision Laboratory, Computer
Science Center, University of Maryland, College Park, MD 20742.
to be equivalent to the definition of Sklansky [10] for thecase of regular digital regions. It was further shown that thedefinition is also equivalent to the definition of Minsky andPapert [5]. Hence, the definition satisfies several existingcriteria for digital convexity.A planar figure is said to be convex if for every pair of points
in the figure, the line segment connecting them lies entirely init. In [9] digital straight line segments were characterized, andit was suggested that convexity of digital regions might bedefined in terms of digital straight line segments. In this paperwe accomplish this by showing that a digital region is convex ifand only if every pair of points in the digital region is con-nected by a straight line segment lying in it. This resultfurther justifies the definition of digital convexity proposed in[4].A planar figure is convex if and only if it has the midpoint
property, that is, the midpoint of any pair of its points alsobelongs to it. Such a property is discussed in [3] ; however,the definition of convexity of digital regions used in [3] isdifferent from ours. We define a midpoint property for digitalregions and show that this property is a necessary conditionfor convexity, but not a sufficient condition. The medianpoint property is defined and shown to be a necessary andsufficient condition for convexity of digital regions.In the next section we introduce our notation and defini-
tions. The relation between the convexity of a digital regionand the digital straight line segment connecting a pair of its
