Blur identification by residual spectral matching - Image ... - IEEE Xplore

Report 2 Downloads 34 Views
IEEE I K A N b A C I IUNb ON IMAGE PRoCtSSING, VOL 2 , NO 2. APRIL 1YY3

141

Blur Identification by Residual Spectral Matching Andreas E. Savakis, Member, IEEE, and H . Joel Trussell, Senior Member, IEEE

Abstract-The estimation of the point function (PSF) is often a necessary first step in the restoration of real images. In this work, a new blur identification method is presented. The PSF estimate is chosen from a collection of candidate PSF’s which may be constructed using a parametric model or from experimental measurements. The PSF estimate is selected to provide the best match between the restoration residual power spectrum and its expected value, derived under the assumption that the candidate PSF is equal to the true PSF. Several distance measures were studied to determine which one provides the best match. The a priori knowledge required is the noise variance and the original image spectrum. The estimation of these statistics is discussed and the sensitivity of the method to the estimates is examined analytically and by simulations. The method successfully identified blurs in both synthetically and optically blurred images.

I. INTRODUCTION

I

MAGE restoration deals with the estimation of the original image from a recorded image which is corrupted by blurring and noise. The standard linear model for the description of the degraded image is [l]

h ( k . l ) f ( i- k , j - 1 )

g(i.j) =

+ n,(i.j)

(1)

(k.l)ERS

where f ( i , j ) is the original image, h ( i . j ) is the spaceinvariant point spread function (PSF) with region of support Rs, and n ( i , j ) is the noise process, assumed signalindependent, white Gaussian. Despite the large variety of image restoration algorithms [ 11, [2], most methods assume a priori knowledge of the PSF. In many practical situations the PSF is not known a priori and must be estimated from the degraded image itself. Therefore, blur identification constitutes an important step in the image restoration process. Early blur identification methods were based on local image characteristics. For example, an isolated bright point in the original image is transformed to the PSF in the blurred image. This approach may be useful with images from astronomy, where bright stars are recorded in a dark background [l], [3]. When an isolated bright point cannot be found, a line or an edge may be used to obtain an estimate of the PSF, provided that the PSF is radially symmetric [l], [3].The estimation of the PSF from local image characteristics is sensitive to noise, and it is not used often [4]. However, the examination of a bright spot in a blurred image can often provide an indication Manuscript received November 23, 1991; revised May 26, 1992. This work was supported in part by a grant from the Eastman Kodak Company. The associate editor coordinating the review of this paper and approving it for publication was Dr. M. Ibrahim Sezan. A.E. Savakis with the College of Engineering and Applied Science, University Rochester, Rochester, NY 14627. H. J. Trussell is with the Deptartment of Electrical Engineering, North Carolina State University, Raleigh, NC 27695. IEEE Log Number 9207227.

of the type of blur which has degraded the image. This information may be used in other blur identification methods, which require knowledge of the blur type to estimate the blur parameters. Another category of blur identification methods is based on the PSF spectral nulls for the estimation of the PSF parameters. These methods are popular because they are easy to implement. Spectral methods estimate the first PSF spectral null from the first local minimum in the recorded image spectrum [ 5 ] . Cepstral methods assume periodic nulls in the PSF power spectrum which result in pronounced spikes in the power cepstrum [6]. These methods do not explicitly take noise into consideration and perform sectioning and averaging to reduce the effects of noise. The bispectrum was proposed as an alternative to the spectral methods [7]. It can be less sensitive to noise, and provided that the data set is large, extensive averaging can suppress the noise effects. All of the above methods are limited to PSF’s which have spectral nulls, such as linear motion and out-of-focus blurs. Recent work in blur identification has focused on simultaneous estimation of the image and blur parameters. The problem has been formulated as a two-dimensional autoregressive moving average (ARMA) parameter estimation problem, where the autoregressive (AR) coefficients represent the image, and the moving average (MA) coefficients represent the blur [SI, [9]. Maximum likelihood estimation has been used for the identification of the ARMA image and blur parameters [lo]. The Expectation-Maximization (EM) algorithm has been applied in the space [ l l ] and frequency domain [12] for iterative maximization of the likelihood function. The advantage of the EM algorithm over other methods is that it requires the solution of linear equations which are decoupled for the image and blur parameters. Another approach is generalized cross-validation (GCV) [13]. It does not require the Gaussian noise assumption, and utilizes the restoration residual for the minimization of the GCV measure over the image and blur parameter space. The methods which simultaneously identify the image and blur parameters require many computations, and using gradient descent may result in convergence to local minima in the search space. In this work, a new blur identification method is presented which does not directly utilize the ARMA modeling approach, but it can incorporate such models if desired. This method can be implemented easily, and has been applied to a variety of PSF types [14]. The a priori knowledge required is the blur type and the second order statistics of the input image and noise, as required in Wiener filtering. The estimation of the blur is decoupled from the estimation of the image by obtaining an estimate of the original image power spectrum

1057-7149/93%03,00 0 1993 IEEE

142

before the identification of the blur. The PSF estimate is chosen from a collection of candidate PSF’s which can be constructed from experimental measurements or using a parametric model. The PSF estimate is chosen to provide the best match based on various distance measures between the residual power spectrum and its expected value, derived under the assumption that the PSF under consideration is equal to the true PSF. In that sense, the method is related to generalized crossvalidation which also uses the restoration residual. However, the matching of the residual spectra is unique to our method. The estimation of the PSF reduces to the minimization of an expression over the set of candidate PSF’s and does not require performing the restoration for each candidate PSF. An exhaustive search is performed throughout the set of candidate PSF’s, so that the one closest to the true PSF can be chosen as the estimate. The exhaustive search approach is preferred over gradient descent to handle experimentally measured PSF’s for which a parametric model is not available, and avoid the problem of convergence to local minima in the search space. The outline of this work is the following. Section I1 discusses parametric PSF models which can be used for the construction of the set of candidate PSF’s. In Section 111, an expression for the power spectrum of the restoration residual is derived under the assumption that the true PSF is used in the restoration. A more accurate form of that expression is obtained by taking into consideration the periodogram spectral estimator that is used. The distance measures which are used to compare the residual spectra are presented. Section IV discusses the estimation of the noise variance and original image spectrum. The sensitivity of the method to the estimates of these quantities is examined. The results for both synthetically and optically blurred images are presented in Section V. A summary is given in Section VI.

A. Out-of-Focus Blur

The PSF for thc out-of-focus blur is usually modeled using geometrical optics. When the aperture is circular, the PSF is modeled as a uniform disk of diameter D.

where K is chosen such that (4) is satisfied. In this case, the PSF estimation is equivalent to the estimation of the blur diameter. The region of support is often circular, but it may also be rectangular, polygonal, or diamond-shaped. The PSF parameters which need to be estimated are defined by the geometry of the region of support. Physical optics have also been used to obtain more accurate PSF representations, but they result in improved restorations only for high SNR [15], [161.

B. Linear Motion Blur The linear motion PSF is approximated by a rectangle, which is one pixel wide, whose orientation and length are determined from the direction and the extent of the blur h(7.1.8) = h(z.1. d. 4)

={

d. i f d w 5 d / 2 0.

and

(j/%) = t a n 4

elvwbere (6)

where d is the length of the camera motion, and 4 is the angle that the motion makes with the horizontal axis. C. Long Term Atmospheric Turbulence Blur Atmospheric turbulence blurs are stochastic blurs, which for long exposures can be approximated by [17]

11. THE SET OF CANDIDATE PSF’S

(7)

Given that the blur type is known and a parametric PSF model is available, blur identification reduces to the estimation of the PSF parameters. Let

h k ( i . j ) = h(i..j.&.)

They are usually encountered in astronomy and remote sensing. The Gaussian PSF is frequently used as an approximation to the PSF for long term atmospheric turbulence.

(4

where 8 k is the vector containing the PSF parameters. The set

’ of candidate parameter values constitutes the search space whose elements are chosen to sample the range of feasible parameter values in the PSF parameter space. The extent of 9 depends on the severity of the blur. The sampling interval depends on the type of parameters that describe the PSF. Examples of blurs which may be parametrized include the out-of-focus blur, the linear motion blur, and the longterm atmospheric turbulence blur. In all cases, it is assumed that due to incoherent illumination the PSF coefficients are nonnegative, and that blurring does not change the mean value of the image:

elsewhere

(8)

09’

where is the variance of the blur, Rh is the region of support, and K is a normalization constant. Thus the Gaussian PSF can be characterized by a single parameter, its variance.

D. Experimentally Generated PSF’s The set of candidate PSF’s is conveniently generated using parametric models, but it can also be constructed experimentally. For example, when the shape of the PSF support changes irregularly with the area of the aperture, as is the case with many inexpensive cameras, a simple parametric model does not accurately represent the PSF. If the camera which was used to take an out-of-focus picture is available, the PSF’s recorded

143

using various aperture shapes and distances of defocus may be used to construct the collection of candidate PSF’s.

The periodogram of the residual is P R , (U.P I ) = I&(’!L,

111. THE POWER SPECTRUM OF THE RESIDUAL

Once the set of candidate PSF’s has been constructed, it remains to select the best estimate from the set. The criterion for selection considered in this work is matching the residual power spectrum with its expected value. These statistics are discussed is this section.

A. Residual Spectrum using Ideal Spectra The restoration residual is expressed in the frequency domain as

where R,(u,v) is the Fourier transform of the residual, and G(u,w ) , P7(u,w),and H,(u.w) are the Fourier transforms of the degraded image, the estimate of the original image, and the PSF obtained by assigning specific values 8, to the PSF parameters. The form of the residual statistics depends on the type of restoration filter. The restoration filter should be chosen so that only when the tested PSF parameters are equal to the true PSF parameters certain terms of the residual are minimized. The choice of restoration filter was restricted to methods for which an analytic expression of the residual could be obtained, such as the inverse filter, the linear maximum a posteriori probability (MAP) filter, and the Wiener filter. The inverse filter and the linear MAP methods resulted in a residual whose form does not allow the true PSF to be differentiated from other PSF estimates. A variation on the Wiener filter provided the most suitable analytical expression of the residual power spectrum. Let f l z ( u , ~be) the Wiener estimate of the original image in the frequency domain. Then,

where H,(u, 7 1 ) is tested under the assumption that it is the true PSF, and Pj(?ilw ) and P,(u, 71) are the power spectra of the original image and noise respectively. Assuming that H , ( U ,? J ) is equal to the true PSF H ( u ,? I ) , and using the relationship between ideal power spectra: Pg(u.?i) = ( H ( ~ . v ) l ~ P f (4~ P71(7i.~) ,v) (11)

the denominator in (10) can be expressed as Pg(ii.71). The Fourier transform of the residual becomes

-

2I)?

I G (U . 7 1 ) I (Pg(U.i i ) - I H , (U,U)I2Pj( U , U)) Pg2(1L.1 ) ) (13)

Given that the ideal power spectra Pg(u.71) and P f ( i i ,TJ) are known, the relationship (11) can be used to expand the numerator of (13) (see (14) at the bottom of the page). If the parameter values tested in H,(u, w) are correct, then H,(?L,1 1 ) = H(7s. w ) , and there is cancellation in the numerator. Thus the expected periodogram of the residual is

The quantities Pg(U. 7 1 ) , P f ( u ,U), and P,(u, U ) are required for the calculation of (13) and (15). The noise is assumed white Gaussian with variance U : , and the noise power spectrum is estimated by the noise variance. The calculation of Pg(u,w ) and Pf(u,71)depends on the data set and the spectral estimation method employed. The periodogram spectral estimator can be used to obtain

P q ( u , , ~=) IG(u.u)I2 P f ( %? I ) = /Fp(?L.?1)(2 denotes estimate, JG(v.7 1 ) l2 is calculated

(16)

(17)

where hat directly from the data, and (Fp(ii. v ) l 2 is an estimate of Pf(ii,7 1 ) using the periodogram of a prototype image; it is possible to use other spectral estimates, such as the AR spectral estimate of P ~ ( uT I.) , as will be discussed in the next section. Using (16) and (17), the periodogram of the residual is calculated as follows

The periodogram of the residual given that JH(U.U ) J 2 becomes

IH,(u,w)~~ =

B. Residual Spectrum using Periodogram Estimates

A more accurate expression than (19) is derived without assuming a priori knowledge of the ideal power spectra Pf and Pg, but instead using the periodogram estimator in be) the periodogram the derivation. Let Q R , = P R , ( P L , w of the residual calculated using periodogram estimates. In

IEEE TRANSACTlONS ON IMAGE PROCESSING, VOL. 2, NO. 2, APRIL 1993

144

(18), I G ( ~. )L I2 . is the only quantity that depends on the The true PSF, and can be expanded in terms of IH1(71,.u)l2. expected residual spectrum Q E , ( U , w) given that the PSF under consideration is the tru? PSF is derived in Appendix A, and is shown to be

This is a more accurate estimate of the expected periodogram of the residual, because it takes into consideration the periodogram spectral estimator employed. The two expressions are equivalent when there is no noise, at the zeros ~ , the noise power dominates over the of I H ( ~ , u )or~ when blurred signal power. For each H , (71. U ) that is considered, QR,( U .71) is calculated independently of Q E ,(U.7 1 ) using (18) and (20), and the two are compared. The H , (71.71) which yields closest agreement between, them is chosen as the blur function that best estimates H ( u . 71). The degree of "closeness" between Q R ,(U.71) and Q E , ( ~ LT I ), is determined using any of several distance measures, which will be discussed shortly. It is of interest to note that the expected value of &n,( 7 1 . 7 1 ) is equal to Q E , ( u .U ) only when the periodogram of the candidate PSF is equal to the periodogram of the true PSF. Let AH be the difference between the periodograms of the It'' candidate PSF and the true PSF:'

C . Distance Measures The distance measures based on the mean square error, the chi square test, the Kolmogorov-Smirnov test, and inner product were considered for comparing QR, and QE,. Some normalization of the spectra is required before carrying out the tests. The purpose of the normalization is twofold: (i) to reduce the dynamic range of the power spectra and (ii) scale the spectra so that they always have the same energy level. The residual spectra that are to be compared exhibit a wide dynamic range, and deviations close to the maxima can dominate the test outcome. The first step in the normalization process is to take the logarithm of the residual spectra. The logarithm function reduces the dynamic range of the two functions, and makes comparisons correspond more to the overall shape of the functions, i.e., deviations at both maxima and minima. For notational convenience let SR(~L U ), = l o g Q ~ ( u . 7 1and ) S E ( ? / .U ) = l o g Q ~ ( u , i i )The . next step is to scale s ~ ( u . 7 1and ) S E ( U . U )so, that the sums of all their elements are equal. If this condition is not met, the largest deviation may occur in cases where the largest variation in energy exists, and not in cases where the largest difference in function shape is found. Since in the cases of measures based on statistical tests the spectra are treated as probability functions, we let

(U.L.)

(11.1.)

The following tests were considered: The Mean Square Error Test

Then

The Chi Square Test The chi square test is a statistical goodness-of-fit test between observed and expected probability density functions:

and

The Kolmogorov-Smirnov Test The Kolmogorov-Smirnov test is a statistical test for comparing cumulative distribution functions:

I

The expected value of Q R , is

(?,~)=(0.0),(0.1):.~.(N-l.N-l).

In general, & ( Q R , )is different than Q E , . They are equal only when 1H,I2 = Iq2.i.e., AH = 0. Therefore, on the average, the best match between QR, and Q E , is obtained for the candidate PSF which is equal to the true PSF. 'The independent variable5 ( I I

I

) are omitted for notational convenience

(28)

The index ( i , j ) indicates the last frequencies in the partial summation of the cumulant. This test depends on the order that the elements are added, which is arbitrary for two-dimensional signals. The order used here was obtained by scanning the image row by row from top to bottom. The Angle Test The angle test determines the angle between the two

.-

n~ r r n

r n r ~ r r 1 r ~ ~ A T In AvN n c C l n T I A 1 CDCPTRAI

spectral vectors by computing their inner product

MATCHING

I45

Sensitivity of Q E to Spectral Estimates The expected periodogram of the residual computed using erroneous estimates is obtained by substituting (30) and (31) into (20). QE,

= QE, f A Q E ,

where

In all tests, a smaller value implies a better match between the power spectra. The effectiveness of the tests was evaluated by simulations and will be discussed in Section V.

IV. ESTIMATING UNKNOWN STATISTICS The calculation of the residual spectra in (18) and (20) requires knowledge of the noise variance and the original image spectrum. These quantities are not known a priori, and must be estimated. The noise variance may be estimated from a low contrast region of the degraded image. The power spectrum of the original image may be estimated from the periodogram of a prototype image, or using AR spectral estimation based on the degraded image or a prototype. The use of a prototype for the estimation of the original image spectrum is based on the assumption that the original image and the prototype are different realizations of the same ensemble, and possess the same statistical properties, one of which is the power spectrum. The periodogram spectral estimator can be used, but it may be affected by local image characteristics which are not present in the original image, such as lines or edges at different orientations. Welch’s spectral estimator [18] may be used to reduce the variations due to local image characteristics. A low-order AR spectral estimator can be used instead of the periodogram in order to obtain a smooth monotonically decreasing spectral estimate. The assumption made in AR spectral estimation is that the power spectrum of the original image is smooth without any nulls. If this is the case, a loworder AR spectral estimate can adequately describe the original image spectrum. The AR spectral estimate may be obtained from a prototype or from the degraded image.

A. Sensitivity Analysis Since the noise variance and the original image power spectrum need to be estimated, it is appropriate to investigate the sensitivity of the method of their estimates. Let the noise variance estimate be

As expected, this term vanishes when 6, = 0 and E f = 0. Let us consider the effects of the errors 6, and E f separately. When there is no error in the estimation of the noise variance, the term due to errors in the original image spectrum is (33) This term is directly proportional to [HI2and becomes negligible at the PSF spectral zeros. It is also directly proportional to E f . For most images, the power spectrum decreases as the frequency increases, and E f is also expected to decrease with frequency, assuming that E , changes proportionally to IFI2. When there are no errors in the estimation of the original image spectrum, the term due to the error in the noise variance is

The error in the noise variance estimate is the same for all frequencies, because of the white noise assumption. The variation of AI- with frequency is due to l/lGI2 and IH121F12. At the PSF spectral zeros, the last term is negligible, but the first two terms assume their maximum value, because 1/(G\’ is large. Since AE is small, the total error is more likely to be dominated by AI-. In the lower frequencies, where l/lGI2 is small, the total error is more likely to be dominated by AE, provided that the error E f is in the same order of magnitude as IF/’.

Sensitivity of Q R to the Spectral Estimate The periodogram of the residual is sensitive to errors in the original image spectrum, but it does not depend on the noise variance. The effect that errors in the estimation of the original image spectrum have on the periodogram of the residual is calculated by substituting (31) into (18).

where 6, is the error in the noise variance estimate. Let IFPI’ be the periodogram of a prototype image, and’

QR,

where

where E.f is the error between the periodograms.

= QR,

Ad)R,

146

IEEE TRANSACTIONS ON IMAGE PROCESSING. VOL. 2, NO. 2, APRIL 1993

*

10

3

g

U2

[cameraman periodograml

10” (chips periodogram]

lom lo9 lo8

10’ Fig. 1.

Original ( 2 . X x 2 j G ) ”camcraman” image.

I

I

0

20

I

I

40 80 Frequency index

I

80

Fig. 3. Radially averaged spectral estimates: periodogram of “cameraman” image, periodogram of “chips” image, and AR spectral estimate of degraded “cameraman” image.

average, the two periodograms are similar, as shown by the radially averaged spectra in Fig. 3. The averaged periodogram spectral estimates are also similar to an AR spectral estimate obtained from the degraded image, blurred with an out-offocus blur with diameter 8 pixels at 40 dB SNR, also shown in Fig. 3. The results for an out-of-focus PSF with diameter 8 pixels Fig. 2. Prototype (256 x 2%) “chips” image. are presented in Table I. Simulations were done with SNR from 40 to 10 dB. The PSF parameter estimates obtained This term vanishes if E f = 0 and at the zeros of IH\*. using the four tests are shown under the columns M for the While both Q E and Q R depend on errors in the estimation mean square error test, C for the chi square test, K the for of the original image spectrum, only Q E depends on errors in Kolmogorov-Smirnov test, and A for the angle test. A (-) the estimation of the noise variance. Simulations have shown indicates that the PSF parameter estimate obtained was at the that the error term AE in (33) may become larger than Q E , boundary of the search space. but the error term Ab- in (34) does not. Therefore, it is The sensitivity to the noise variance estimate was examined expected that the estimation of the original image spectrum by using the true noise variance under heading “Noise variance is much more important than the estimation of the noise known,” and by using the noise variance estimate obtained variance. In the next section, the sensitivity of the method to from a low contrast region of the degraded image under the estimates of the noise variance and original image spectrum heading “Noise variance estimated.” The degraded image is is demonstrated by simulation examples, where it is verified divided in small sections, the sample variance of all sections that the method is more sensitive to the estimate of the original is computed, and the smallest sample variance is chosen as an image spectrum. estimate of the noise variance. The noise variance estimates were usually within 20% of the true noise variance. The V. RESULTS PSF parameter estimates do not significantly change when the The identification of the PSF by residual spectral matching estimate of the noise variance is used instead of the true noise was tested using synthetically blurred images under a variety variance. The sensitivity to the estimate of the original image specof conditions for both one-dimensional signals and images. trum was examined by obtaining results using the periodogram The set of candidate PSF’s was constructed by varying the PSF parameter (focus diameter, motion length, or Gaussian of the original image shown under column “Per. Original,” variance) from a value of two pixels to twice the true PSF the periodogram of the prototype shown under column “Per. parameter value in one pixel increments. Quantization of one Prototype,” and an AR spectral estimate obtained from the pixel of the PSF parameter is usually sufficient for image degraded image shown under column “AR Degraded.” As expected from the sensitivity analysis, the results are affected restoration purposes. Representative results are presented for the 256 x 256 by the accuracy of the estimate of the original image spectrum. “cameraman” image shown in Fig. 1. A prototype which When the periodogram of the original is known, a good was similar to the “cameraman” image was not available, PSF parameter estimate may be obtained for SNR as low as and the one used is the “chips” image shown in Fig. 2 . 10 dB. When the periodogram of the prototype is used, the The periodograms of the “cameraman” and “chips” images estimation of the PSF parameters becomes less accurate. In primarily differ due to edges at different orientations. On the most cases, a PSF parameter estimate within two pixels of

147

TABLE I BLURIDENTIFICAIION RESULTSOUT-OF-FOCUS PSF 8 PIXEL DIAMETER ORIGINAL IMAGE. "CAMERAMAN", PROTOI-YPE IMAGE: "CHIPS" NOISEVARIANCE KNOWN

Per. Prototype

Per. Original SNR

M

C

K

40 dB

8

8

8

20 dB

8

8 8 8

8

30 dB

8 7

7 9

8

x 7

10 dB

7

7

A

9

6

M

AR Degraded

K

A

M

8

8

8

8

8

8

8

h

8 7 7

8

8

8 7

8

8

10

10

9

8 10

4

-

-

-

-

-

7

C

K

C

A

Noise variance estimated

Per. Original SNR

M

C

K

Per. Prototype A

M

C

K

AR Degraded A

M

C

K

A

8

8

40 dB

8

8

8

X

8

8

8

8

30 dB

8

8

8

7

8

8

7

7

8 8

8

8 8

10

10

-

9

-

-

-

20 dB

8

8

7

8

6

6

7

6

10 dB

10

10

12

9

h

6

-

-

TABLE 11 BLURIDENTIFICATION RESULTS. OUT-OF-FOCUS PSF'S, 5 ORIGINAL

AND

~

8

12 PlXtL DIAMETER.

NOISE VARIANCE ESTIMATED IMAGE " C A M E R A M A N PROTOTYPE IMAGE "CHIPS"

~~~

Per Prototype

Per Original SNR

M

C

K

40 dB

5 5 4 7

5 5

30 dB 20 dB 10 dB

A

M

C

K

A

M

C

K

A

5

5

4

4

4

4

6

4

4

3

4

5 4

5 5

4

4

4

6

6

7

7

7

7

4 4

6

7

2 2

5 5 6

5 5

4

5 5

-

-

-

Per Prototype

Per Original SNR

M

AR Degraded

C

K

A

M

C

K

~

AR Degraded A

M

C

K

A

40 dB

12

12

12

12

12

12

12

12

12

12

11

12

30 dB

12

12

12

14

12

12

12

12

12

12

11

12

20 dB

12

12

11

11

11

11

12

9

18

18

14

-

10 dB

14

14

14

16

11

11

8

-

-

-

-

-

the true PSF parameter value is obtained for SNR as low as 10 dB. When the AR spectral estimator is used based on the degraded image, the method is limited to a SNR of 30 dB. The poor performance of the method when using the AR spectral estimator may be attributed to (a) the lack of accounting for the blurring and noise in the data, and (b) the fact that the derivation of the expected residual spectrum is based on the periodogram spectral estimator, but the AR spectrum is not an accurate estimate of the periodogram. Results for other out-of-focus blurs are shown in Table 11, for linear motion blurs in Table 111, and for Gaussian blurs in Table IV. These results are typical of simulations with many images and various blur parameters and noise levels. The simulation results are also used to evaluate the performance of

the tests of distance measure. The MSE test and chi square test performed the most consistently, while the angle test yielded the worst results. Let us call the function obtained by the outcome of a test in the parameter space the error function. There are two characteristics of the error function which are of particular interest: (i) the effects of the SNR on the error function and (ii) the existence of local minima on the error function. As expected, by decreasing the SNR the global minimum in the error function becomes less pronounced. This effect is illustrated in Fig. 4, where the MSE error function is plotted for various SNR levels. The broadness of the error function gives an indication of the degree of confidence in the estimate. As the SNR decreases, the error function becomes flatter, and

148

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 2, NO. 2, APRIL 1993

TABLE 111 BLURIDENTIFICATION RESULTS.HORIZONTAL LINEARMOTIONPSF'S 8 AND 15 PIXEL NOISEVARIANCE ESTIMATED ORIGINAL. IMAGE: "CAMERAMAN"; PROTOTYPE IMAGE "CHIPS"

LL:NGTH

_____

Per Original A

M

C

8

8

8

8

8

8

8

8

9

8

8

9

8

8

9

8

8

7

8

10

9

6

6

-

7 5

M

C

40 dB

8 8 8 8

20 dB 10 dB

Per. Original SNR

M

K

\

M

C

8

6

8

8

A

K

A

6

8

6

8

8

8

-

-

8

-

-

-

-

-

AR Degraded

Per Prototype

C

K

A

M

C

~

AR Degraded

K

SNR

30 dB

Per. Prototype

K

A

M

C

K

A

40 dB

15

15

14

14

15

14

14

14

15

14

15

15

15 15

15

30 dB

14

14

14

15

14

I6

16

15

16

20 dB

14

16

16

16

14

14

15

I6

-

-

15

-

10 dB

16

16

17

18

16

16

13

10

-

-

-

-

TABLE IV BLUR~UENTlFlCATlON RESULTS: GAUSSIA~; PSF'S 8 AND 15 PIXEL VARIANCE NOISEVARIANCE ESTIMAlkD ORIGINAL IMAGE "CAMERAMAN", PROTOTYP~ IMAGE "CHIPS" Per. Original

SNR

M

C

K

A

M

C

K

_______

40 dB

8

30 dB

8

20 dB

9

10 dB

11

8 8 9 11

M

C

40 dB

15

15

30 dB

15

15

20 dB

14

14

10 dB

-

-

A

M

C

K

A

_____

8

8

7

6

6

6

6

6

4

7 9

8

8

7

9

6

8

9

7

6

6

5

6

15

8 15

12

6 5 15

14

10

9

9

5

-

-

-

-

-

Per. Original

SNR

AR Degraded

Per. Prototype

K

AR Degraded

Per. Prototype

C

K

A

13

13

14

14

14

12

12

12

14

11

11

8

11

13

14

14

8

-

A

M

15

14

15

12

16 8

the confidence in the PSF estimate also decreases. It was experimentally found that without errors in the noise variance and original image spectrum, SNR level of 10 dB is limit at which the error function minimum is distinguishable. The second point that deserves attention is that the error function has local minima. It is of interest to note that when there are no errors in the noise variance and original image spectrum, the local minima usually appear for PSF's with singularities, such as the out-of-focus and linear motion PSF's. Local minima are not likely to appear when the PSF considered has a monotonically decreasing power spectrum without spectral nulls, such as the Gaussian PSF. However, errors in the estimate of the original image spectrum may introduce local minima at other points even in the case of Gaussian blur. The existence of local minima in the error function makes the use of gradient descent-type algorithms

C

K

A

14

14

10

13

15

15

13

19

-

-

-

-

-

-

-

-

M

ineffective, since depending on the point of initialization the algorithm could converge to a local minimum. For this reason, the approach of choosing the PSF estimate from a collection of candidate PSF's was preferred. Optically blurred images were recorded using a CCD camera courtesy of the Eastman Kodak Company. Two different images were recorded at various distances of defocus. The images "text" and "objects" recorded at the focal plane (3m imaging distance) are shown in Figs. 5 and 6. The power spectrum of the original image was estimated by the periodogram of the image recorded at a distance of 3m. The blur identification results based on the MSE test are shown in Table V for various degrees of defocus. The results are compared to the results obtained using cepstral analysis. Restorations using both PSF estimates were not significantly different. The identification of real blurs using PSF models

SAVAKIS AND TRUSSLLL BLUK IUEN I l l - I L A I I U N B Y KkhIUUAL b l ’ t L l K A L

”j 70

2

Fig. 4.

149

MATC I l l N G

\ 4

8 1 0 1 2 1 4 1 6 Blur diameter in pixels

6

MSE error function for “cameraman” blurred with out-of-focus PSF 8 pixel diameter at SNR from 40 to 10 dB.

Fig. 6. Target image “objects” recorded at 3117 imaging distance

TABLE V BLUR[DtNTIFICATION RLSLILTS OPTlCAl LY BLURRED IMAGES BLURDIAMETER ESTlMArlON FOR “TEXT“ AND “OBJECTS” SEQUENrES In-camera blurred images Focal plane at 3m Estimated SNR 20 dB Imaging Dist.

Cepstrum

Res. Stat. (MSE test)

5.00 m

Y

10

2.00 m 1.60 m

12 24

25

1.20 m

31

39

13

Estimated SNR 30 dB Fig. 5.

Target image “text” recorded at 3111 imaging distance.

Imaging Dist.

Cepstrum

Res. Stat. (MSE test)

1.75 m

16

18

5.25 m

16

16

~

based on geometrical optics raises the question of whether PSF models based on physical optics can result in more accurate representation of the PSF and better restorations. Based on the results of [15], improved restorations can be obtained by using PSF’s modeled by physical optics when the SNR is above 30 dB. In the above examples, the SNR was not high enough to expect better restorations due to the use of the physical PSF.

spectral estimate obtained from the degraded image were used to estimate the original image spectrum. The results indicated that the method can be used to estimate the PSF’s in both synthetically and optically blurred images at SNR as low as 20 dB.

VI. SUMMARY

A new blur identification method based on residual spectral matching was developed and tested. The new method does not require parametric PSF models or AR image models, but can incorporate such approaches if desired. The PSF is chosen from a collection of candidate PSF’s to provide the best match between the restoration residual spectrum and the expected residual spectrum given that the candidate PSF is the true PSF. The method requires a priori knowledge of the original image spectrum and the noise variance. Both the analysis and simulations showed that the results are more sensitive to the estimate of the original image spectrum than the estimate of the noise variance. The periodogram of a prototype or an AR

APPENDIXA In this Appendix, the expected value of the residual power spectrum is derived by taking into consideration the periodogram spectral estimate in (16). The periodogram estimates (16) and 17 are used in (13) to obtain (18). In that expression, (G(u.u)1 IS the only quantity that depends on the true PSF. By expanding IG(u. in terms of IH,(w,?))I2, letting lHt(u. .)I2 = I H ( u ,v)I2. and taking expectations, Q E , ( u ,7)) can be derived from Q R , ( u . 7 j ) The . quantity I G ( ~ , Vis) ~ ~ expanded as follows:

I?

())I2

IG1’ = / H F + NI’

150

IEEE TRANSACTIONS ON IMAGE PROCESSING. VOL. 2, NO. 2, APRIL 1993

= [HF

+ N ] [ H F+ NI”

= IHFI’

+ ]NI2+ 2 [ N r ( f W r+

~

1

,

4

w

,

]

where * denotes complex conjugate, and the subscripts T and m denote real and imaginary parts. Thus (see equation at the top of the page) The expected value of the periodogram of the residual is obtained by taking the expected value and using the following identities

€ ( N , ( u . v ) )= 0

(36)

E(Nm(?l,?l)) =0

(37)

€ ( N r ( U ,7l)N,(U. €(N;(?L.71))

Pi))

=0

.7‘.

= 2

€(N,2(U,1/))=

4

-

2

(38) (39)

(40) (41)

-

-

2a;IHF\’

PI2

IGI’

0:

+

which is not the same as P E ~ ( ?1 1L) .in (19). This expression has an extra term compared to PE? of (19). REFERENCES H. C. Andrews and B. R. Hunt, Digital Image Restoration, Englewood Cliffs, NJ: Prentice-Hall, 1977. M. 1. Sezan and A.M. Tekalp, “Survey of recent developments in digital image restoration,” Optical Engineering, vol. 29, pp. 393-404, May 1990. A. Rosenfeld and A. Kak, Digital Picture Processing. New York: Academic, 1982. M. I. Sezan, G. Pavlovic, A.M. Tekalp, and A. T. Erdem, “On modeling the focus blur in image restoration,” in Proc. IEEE In/. Conf: Acoust., Speech, and Signal Proc., Toronto, Canada, May 1991, pp. 2485-2488. D. B. Gennery, “Determination of optical transfer function by inspection of frequency-domain plot,”J. Opt. Soc. Amer., vol. 63, pp. 1571-1577, Dec. 1973. T. M. Cannon, “Blind deconvolution of spatially invariant image blurs with phase,’’ IEEE Trans. Acous., Speech, Signal Proc., vol. A S S - 2 4 , pp. 58-63, Feb. 1976. M. M. Chang, A. M. Tekalp, and A. T. Erdem, “Blur identification using the bispectrum,” IEEE Trans. Acous., Speech, Signal Proc., vol. 39, pp. 2323-2325, Oct. 1991. A.M. Tekalp, H. Kaufman, and J. W. Woods, “Identification of image and blur parameters in the restoration of noncausal blurs,’’ IEEE Trans. Acoust., Speech, Signal Proc., vol. ASSP-34, pp. 963-972, Aug. 1986. A.M. Tekalp and H. Kaufman, “On statistical identification of a class of linear space-invariant image blurs using non-minimum phase ARMA models,” IEEE Trans. Acoust., Speech, Signal Proc., vol. 36, pp. 1360-1363, Aug. 1988. R. L. Lagendijk, A. M. Tekalp, and J. Biemond, “Maximum likelihood image and blur identification: a unifying approach,” Optical Engineering, vol. 29, May 1990. R. L. Lagendijk, J. Biemond, and D. E. Boekee, “Identification and restoration of noisy blurred images using the expectation-maximization algorithm,” IEEE Trans. Acoust., Speech, Signal Proc., vol. 38, pp. 1180-1191, July 1990. K. T. Lay and A. K. Katsaggelos, “Image identification and restoration based on the EM algorithm,” Optical Engineering, vol. 29 pp. 436-445, May 1990. S.J. Reeves and R.M. Mersereau, “Blur identification by the method of generalized cross-validation,” IEEE Trans. Image Proc., vol. 1, pp. 301-311, July 1992.

SAVAKIS AND TRUSSELL: BLUR IDENTIFICATION BY RESIDUAL SPECTRAL MATCHING

[14] H. J. Trussell and A. E. Savakis, “Blur identification by statistical analysis,” in Proc. IEEE In/. Conf: Acoust., Speech, and Signal Proc., Toronto, Canada, May 1991, pp. 2493-2496. [15] A. E. Savakis and H. J . Trussell, “On the accuracy of PSF representation in image restoration,” IEEE Trans. Image Proc., this issue. [ 161 A. E. Savakis, Blur Identification by Statistical Analysis, Ph.D. dissertation, North Carolina State Univ., 1991. [17] R. E. Hufnagel and N. R. Stanley, “Modulation transfer function associated with image transmission through turbulent mcdia,” J . Opt. Soc. Amer., vol. 54, pp. 52-61, 1964. [lS] P. D. Welch, “The use of the Fast Fourier Transform for the Estimation of Power Spectra,’’ IEEE Trans. Audio Electroacoust., vol. AE-15, pp. 70-73, June 1967.

Andreas Savakis wdS born in Athens, Greece He received the B S degree (Summa Cum Laude) and the M.S degree from Old Dominion University, Norfolk, VA, in 1984 dnd 1986, re\pectively, and the Ph.D. degree from North Carolina State University, Raleigh, NC, in 1991, all in Ekctricdl Engineering He is currently a Postdoctoral Research Fellow at the University of Rochester, Rochester, NY, working in the Center for Electronic lmaging and the Rochester lmdging Consortium His resedrch interests include image processing, pattern recognition, and infrared imaging He is a member of Alpha Chi, Tau Beta Pi, dnd the Signdl Processing dnd Computer Societies of the IEEE

151

Joel Trussell (S’75-M’76-SM’91) was born in Atlanta, CA, on February 3, 1945. He received the B.S. degree from Georgia Tech in 1967, the M.S. degree from Florida State in 1968, and the Ph.D. degree from the University of New Mexico in 1976. He joined the Los Alamos Scientific Laboratory, Los Almos, NM, in 1969, where he began working in the image and signal processing dept. in 1971. During 1978-79, he was a visiting professor at Heriot-Watt University, Edinburgh, Scotland. In 1980, he joined the Electrical and Computer Engineering Department at North Carolina State University, Raleigh, NC. During 1988-89, he was a visiting scientist at the Eastman Kodak Company in Rochester, NY. His current interests include signal restorationireconstruction, color reproduction, mathematical methods. He won the IEEE-ASSP Society Senior Paper Award (with M.R. Civanlar). He has been an associate editor for the Transactions on ASSP and is currently chairman of the image and Multidimensional Digital Signal Processing Committee of that Society. He is a member of Sigma Xi, Tau Beta Pi, Phi Kappa Phi, and the Signal Processing Society.

t