Passive Depth From Defocus Using a Spatial ... - Semantic Scholar
Passive Depth From Defocus Using a Spatial Domain Approach D. Ziou Departement de mathematiques et d'informatique Universite de Sherbrooke Sherbrooke, Qc, Canada, J1K 2R1
Abstract
This paper presents an algorithm for a dense computation of the di erence in blur between two images. The two images are acquired by varying the intrinsic parameters of the camera. The image formation system is assumed to be passive. Estimation of depth from the blur di erence is straightforward. The algorithm is based on a local image decomposition technique using the Hermite polynomial basis. We show that any coef cient of the Hermite polynomial computed using the more blurred image is a function of the partial derivatives of the other image and the blur di erence. Hence, the blur di erence can be computed by resolving a system of equations. All computations required are local and carried out in the spatial domain. An algorithm is presented for estimation of the blur in 1D and 2D cases and its behavior is studied for constant images, step edges, line edges and junctions. The algorithm is tested using synthetic and real images. The results obtained are very encouraging.
1 Introduction
In this paper, we are interested in the depth from defocusing approach, which involves calculating the distance between a visible surface of an object and the thin lens of the optical camera by using the degree of blurring. Speci cally, let us consider a camera which generates one or more grey-level images of a static 3D scene. Let us assume the perspective projection of a real world point onto the image. The image intensity at a given pixel (x; y) depends on the distance z between the camera lens and the corresponding point in the scene (Fig. 1). In the case of the thin lens, the image is formed by convolution of the ideal projected image and the camera point spread function (PSF) g(x; y; (x; y)), where (x; y) is the blur parameter at the pixel (x; y). We will consider that (x; y) is constant within a given window, that is (x; y) = . Usually this assumption is not realistic. However, without this assumption we are confronted with a non-shift invariant system. These systems are mathematically
hard to handle. This explains why this assumption has been implicit in most prior work. Many functions are used to approximate the PSF of the camera. The reader can nd more information in [12]. We are interested in the most common case, where PSF is approximated by the Gaussian function: x y (1) g (x; y) = 1 e?
2+ 2 2
2
The relationship between the distance z (Fig. 1) and the blur is given by [7]: ( Fv if z > u (2) z = v?FFv?kf if z < u v?F kf +
where u is the distance between the lens and the position of perfect focus, v the distance between the lens and the image plane, F the focal length, f the f-number, and k the proportionality coecient between the blur circle radius and . This formula shows that estimation of depth z requires the estimation of f; F; v; k, and . The intrinsic parameters of the camera f , F , v, and k are independent of pixel location (x; y). They can be estimated by calibration [8, 13, 2]. The blur at a pixel (x; y) is measured by calculating the PSF function at this pixel. The work of Pentland [6] was the rst research in computer vision to test the use of defocus information as a depth cue. Since Pentland's work, many other schemes for the computation of blur have been proposed. The most common approaches are inverse ltering [6, 7, 2, 1, 13, 14], depth recovery from edges [7, 3, 10], and the S-Transform [11]. Later, we will summarize the idea behind the S-Transform approach because it is directly related to our work. In this paper, we will present an algorithm for computing the di erence in blur between two images. These two images are acquired by varying the intrinsic parameters of the camera. The image formation system is assumed to be passive. Our algorithm involves tting the images by a Hermite polynomial basis. We will show that any coecient of the Hermite
polynomial computed using the more blurred image is a function of the partial derivatives of the other image and the blur di erence. Hence, the blur di erence can be computed by resolving a system of equations. All computations required are local and are carried out in the spatial domain. Furthermore, no image model has been used to derive the algorithm. Consequently, our algorithm yields a dense estimation of blur. In the next section, we will start by describing the blur estimation rule for 1D images. Then we will present the extension of the 1D mathematical model to 2D images. The algorithm for blur estimation in 2D images is described in section 4. Section 5 presents the behavior of the algorithm and the experimental results.
2 Blur Estimation in a 1D Signal
In this section we present the mathematical development behind blur estimation in the case of a 1D signal. Let us consider Ic (x) and Ib (x), two unfocused images of the same scene obtained by varying one or more camera parameters (see Fig. 1). Let us assume that Ib (x) is more blurred than Ic (x). The image which has the smaller grey-level variance is the more blurred one. When the PSF is Gaussian, the relationship between these images is:
Ib (x) = (Ic g )(x)
(3)
where g (x) is a Gaussian of variance and p represents convolution. The parameter = b ? c is the blur di erence between the two images, where b and c are the blur of Ib (x) and Ic (x) respectively. The relationship between and z can be deduced. For example, let us consider the image Ib (x) obtained by changing the parameter v by an amount v in (Eq. 2), which causes a change of magnitude in the blur . The depth z is given by: 2
2
2
8 > ? 1 =< F z > : F? 1
v v (1 + 1
2 +
q
2
2
1 + f F vvv k ) if z > F
q
2 (2 + 2
)
2
2
f v v v v (1 ? 1 + F v k ) if z < F (4) For more details see [15]. This equation shows that the depth can be computed from the intrinsic camera parameters and the blur di erence . It should be recalled that the camera parameters can be obtained by calibration. Now, we are faced with the problem of estimating . To this end, we have developed an image decomposition technique, called the Hermite polynomial transform, in which the image is locally approximated by Hermite polynomials. The scheme involves tting the image, within the Gaussian window g (x) 1
1
2 +
2 (2 + 2
)
2
2
in (Eq. 1), by a Hermite polynomial basis. The nth Hermite polynomial Pn (x) is de ned by:
n x x (5) Pn (x) = (?1)n e d(dx )n e? The sequence fPn (x)g is orthonormal with respect to g (x) on (?1; +1). Our choice of the Hermite poly2 2
2 2
nomial is due to the fact that these polynomials are associated with the Gaussian function. Hence, they have interesting properties and their use in image processing is increasing [5, 4]. In this mathematical model, the complete description of an image requires that the tting process be repeated at a sucient number of pixels within the window. We will consider the case of equidistant sample spacing T . The approximation of image Ib (x) within this window by the polynomial basis fPn (x)g is: 1 X K X
Ib (x) =
n k ?K =0
cn (kT )Pn (x ? kT )
(6)
=
where 2K + 1 is the maximum number of samples in the window and cn (kT ) is de ned by: cn (kT ) = (Ib hn )(kT ) (7) where n x hn (x) = Pn (x)g (x) = (?1)n d(dx )n p1 e? (8) Substituting (Eq. 3) in (Eq. 7), we obtain: cn (kT ) = (Ic g hn )(kT ) = (Ic fn )(kT ) (9) where n x f (x) = (g h )(x) = (?1)n n d p1 e? (10) 2 2
n
n
dxn
2 2
and = + . (Eqs. 7 and 9) show the relationship between image Ib and image Ic . We can use these two equations to estimate the blur. However, we are faced with solving a non-linear equation. Now we want to give a simple form of the relationship between these two images which will make the computation rule of explicit. Let us consider the Taylor expansion of Ic (kT ? x) at kT : 2
2
2
Ic (kT ? x) =
1 (?x)p I p (kT ) X c +
p
( )
=0
p!
(11)
where Icp (x) denotes the pth derivative of Ic (x). Substituting this equation in (Eq. 9), we obtain: ( )
cn (kT ) =
X1 +
p
=0
sn;p Icp (kT ) ( )
(12)
where
Z 1 (?x)p f (x) n dx sn;p = p!
(13)
?1
(Eqs. 7 and 12) show the relationship between image Ib and image Ic . In this relation, a coecient sn;p is proportional to the pth moment of fn (x) and is a function of the blur . In other words, given Ic (x), Ib (x) and fn (x) we are able to estimate . However, in momentbased schemes, calculating the moment directly entails a high computation cost. To overcome this diculty, we propose an ecient computation procedure for sn;p [15]:
8 > < sn;p = > :
p?n
(?1)nn p?n?p?np?n?p?n= if 0 < n < p and p ? n is even (?1)nn if n = p 0 otherwise (
1)(
3)
)!2(
(
3
(14)
This equation is an ecient computation rule for
By using (Eqs. 7, 12, and 14), the relation between the two images Ib (x) and Ic (x) is now given by:
8 P 1 s I p (kT ) = (I h )(kT ) > < c r (kT ) = P b r pr r; p c p 1 > : (cIbr h(kTr ) )(=kT )pr s r ; p Ic (kT ) = (2 )
+
2
2
2
(2 +1)
+
2 +1
2 +1 2 +1
2 +1
(15) We are now able to compute the blur di erence. For the sake of clarity, let us examine the use of coecients c and c in (Eqs. 15) Pas an example. In the case of c , we assume that p> s ; p Ic p (kT ) = 0 and P for c , we assume that p>1 s ; p Ic p (kT ) = 0. According to (Eq. 14), s ; = 1, s ; = ?, s ; = , and s ; = ? . Thus (Eq. 15) becomes: 0
1
(2 )
0
1
0 2
1
1 2 +1
(2 +1)
+
1
0 0
2
(
1 3
2
1 1
0 2
4
4
c (kT ) = Ic (kT0 ) + Ic00 (kT000) = (Ib h )(kT ) c (kT ) = ?Ic (kT ) ? Ic (kT ) = (Ib h )(kT ) 2
0
4
2
1
0
(16) where h (kT ) and h (kT ) are given in (Eq. 8). By resolving the two equations in (16) we obtain respectively: 1
s
= 4 (Ib h )(I 00kT(kT) ?) Ic (kT ) ? c 0
or
s
Ic0 (kT ) ? = 4 (Ib h ?)(IkT000)(+ kT ) 1
c
0
2
0
2
(17)
(18)
+
=0
cn;m (kTx; lTy ) =
=0
=
1X 1 X +
+
p
=0
q
=
sn;p;m;q Icp q (kTx ; lTy ) ( )( )
=0
(21) where sn;p;m;q is de ned by: Z 1 Z 1 (?x)p (?y)q f (x; y) n;m sn;p;m;q = dxdy p ! q ! ?1 ?1 (22) fRn;m(x; y)p is separable, thus s = s s = n;p;m;q n;p m;q 1 ?x fn x dx R 1 ?y q fm y dy. Rapid calculap q ?1 ?1 tion of sn;p;m;q is possible by using (Eq. 14). The coecient cn;m (kTx; lTy ) in (Eq. 21) becomes: +
(
)
( )
(
)
!
cn;m (kTx; lTy ) =
1
4
0
The extension of the 1D scheme for blur estimation to the 2D case is straightforward. In fact, a 2D Hermite polynomial can be created by taking the tensor product of two sets of 1D Hermite polynomials; i.e, fP (x)P (y); ; Pn (x)Pm (y); g, where n is the degree of x and m the degree of y. Since a Gaussian window g (x; y) is also separable, it is easy to show that fPn;m(x; y)g is orthogonal within g (x; y). Let us consider that kTx and lTy are the sample spacings in the window along the x and y axes. The decomposition of image Ib (x; y) into localized polynomials is given by: P P P P Ib (x; y) = n 1 m1 Kk ?K Ll ?L cn;m (kTx; lTy ) Pn;m (x ? kTx ; y ? lTy ) (19) where (2K + 1) (2L + 1) is the size of the sampling lattice and cn;m(kTx ; lTy ) is de ned by: cn;m (kTx; lTy ) = (Ib hn;m)(kTx ; lTy ) (20) where hn;m (x; y) = Pn;m (x; y)g (x; y). Both Pn;m (x; y) and g (x; y) are separable; thus straightforward manipulations, as in the 1D case, give: +
) 2
sn;p . The simpli cation of sn;p also simpli es the computation of cn (kT ) and the approximation of Ib (x).
2
3 Extension to a 2D Image
( )
!
1X 1 X +
+
pn qm
sn;p;m;q Icp q (kTx ; lTy ) ( )( )
which can be simpli ed to obtain:
8 > > > > > > < > > > > > > :
P 1P +
(23)
1
+
c2r;2t (kTx ; lTy ) = pr qt s2r;2p s2t;2q (2p)(2q ) (kTx ; lTy ) Ic +1 +1 c2r;2t+1 (kTx ; lTy ) = pr qt s2r;2p s2t+1;2q+1 (2p)(2q +1) (kTx ; lTy ) Ic +1 +1 c2r+1;2t (kTx ; lTy ) = pr qt s2r+1;2p+1 s2t;2q (2p+1)(2q ) (kTx ; lTy ) Ic +1 +1 c2r+1;2t+1 (kTx ; lTy ) = pr qt s2r+1;2p+1 (2p+1)(2q +1) (kTx ; lTy ) s2t+1;2q+1 Ic
P P
P P
P P
(24)
derivatives Icp q (x; y) can be estimated using di erentiation techniques such as Sobel masks, Gaussian masks, wavelets, polynomial tting, etc. The choice of algorithm parameters Tx , Ty , N , M , P and Q is discussed in [15]. ( )( )
We are now able to compute the blur di erence between two 2D images. For clarity, let us take as an example the coecients c ; , c ; , and c ; in (Eqs. 24). Let us assume that the image is a cubic polynomial; that is Icp q (kTx ; lTy ) = 0; 8p + q > 3. By using (Eqs. 14 and 20), we can rewrite the rst three equations in (24): 0 0
0 1
1 0
( )( )
5 Performance Evaluation
5.1 Algorithm behavior
8 c ; (kTx ; lTy ) = Ic (kTx; lTy ) + rIc (kTx; lTy ) = > > (Ib h ; )(kTx ; lTy ) > > < c ; (kTx; lTy) = ?Ic (kTx; lTy )? > Ic (kTx ; lTy ) = (Ib h ; )(kTx ; lTy ) > > (kTx; lTy )? x ; lTy ) = ?Ic > : c ; (kT Ic (kTx ; lTy ) = (Ib h ; )(kTx ; lTy ) (25) where r is the Laplacian. The rst equation in (25) 2
0 0
4
0 0
(0)(1)
0 1
2
(2)(1)
0 1
4
1 0
(1)(0)
2
(1)(2)
1 0
4
0
is the same as the one implemented by M. Subbarao and G. Surya [11]. These authors propose a spatial approach, called the S-Transform. They assume that the by the cubic polynomial p(x; y) = P image P ?ism tted a x n;m m y n . In this case, they show that m n Ic (x; y) = Ib (x; y) ? rIb (x; y)=4, where r is the Laplacian and the blur di erence. The estimation of , and hence of the depth, is straightforward. It is easy to see that their model for blur estimation is a particular case of our blur estimation model presented in (Eq. 15). 3
=0
0 0
2
2
( )( )
2
4 Algorithm for Blur Estimation
( )( )
We have presented the mathematical developments required for the computation of blur. Our next step is to specify an algorithm for the estimation of the blur di erence , based on these mathematical developments. Let us consider a nite summation in (Eqs. 20 and 24) for a given n and m; that is p 2 [r; P ] and q 2 [t; Q]. The blur is estimated by using (Eq. 23) and any equation in (24). However, there are four equations in (24); thus we estimate by using (23) and each equation in (24). In each window W centered at (x; y), we obtain four values of , from which we choose the one which minimizes the following quadratic error: e2 (x; y ) =
XX r
s
(Ib (x + r; y + s)
( )( )
? (Ic g )(x + r; y + s))
2
(26)
To summarize, given , Tx , Ty , N , M , P , and Q, the algorithm involves three steps: 1) the estimation of cn;m using (Eq. 23); 2) the estimation of four values of using (Eq. 24); 3) the selection of the best according to the quadratic error in (Eq. 26). The image
0
0
3
=0
Up to now, no assumption has been formulated about the images for which our algorithm performs well. In this section, we will examine its behavior in the case constant images, step edges, line edges, and junctions. 1) Constant images: the derivatives of the image are zero and the system of equations in (23) and (24) is reduced to sn; ;m; Ic (kTx; lTy ) = (Ib hn;m)(kTx ; lTy ). According to (Eq. 14), sn; ;m; is one if n = m = 0 and zero elsewhere, thus Ic (kTx; lTy ) = (Ib h ; )(kTx ; lTy ). This equation does not permit computation of blur. In other words, low frequencies are not used by our algorithm to estimate blur. 2) Step edges: we consider the step edge model, I (x; y) R=x ((x ? kTx )cos() + (y ? lTy )sin()), where (x) = ?1 exp(?t )dt and is the edge orientation. I p q (kTx; lTy ) = 0, if p or q is even. In this case, only the last equation in (24) can be used to estimate blur. 3) Line edges: we consider the line model, I (x; y) = e?((x?kTx )cos()+(y?lTy )sin()) , where is the line orientation. I p q (kTx ; lTy ) = 0, if p or q is odd. In this case, only the rst equation in (24) can be used to estimate blur. 4) Junctions: for the L junction model, I (x; y) = (x ? kTx ) (y ? kTy ). I p q (kTx ; lTy ) = 0, if p or q is even. In this case, only the last equation in (24) can be used to estimate blur. This algorithm exhibits the same behavior in the case of a T junction, I (x; y) = (x ? kTx )( (y ? kTy )+1). Let us examine the behavior of the rst equation in (Eq. 25) implemented by Subbarao and Surya, where the image is t by a polynomial basis of degree less than or equal to three. At a step edge, the Laplacian is zero and therefore this equation becomes c ; (kTx; lTy ) = Ic (kTx; lTy ). Thus, it is not possible to compute the blur at a step edge. At a line edge, the Laplacian is not zero and it is possible to compute the blur. For linear junction models (the intersection point of n sectors with constant intensities such as L, V, T, Y, and X), the Laplacian of Gaussian is zero at the junction point. Hence, it is not possible to estimate the blur at the junction point. We conclude that Subbarao and Surya's algorithm cannot be used to estimate blur at step edges and junctions. 0 0
0
5.2 Experimental Results
We have tested the equations in (25) with two real images (Fig. 2.a and 2.b) of a toy house taken using an optical camera with di erent aperture and exposure time settings. The images have 384 320 pixels and 256 grey levels. They were given to us by Y. Xiong of the Robotics Institute at Carnegie Mellon University. The toy house is about two meters away from the camera. It is about 2.5 inches wide, 3 inches long, and 4.5 inches high. The partial derivatives of the image are computed by convolving the image with the appropriate derivatives of the Gaussian function. The scale is 1, the spacings Tx = Ty = 1. Figures 2.c, 2.d, and 2.e present the blur estimated using c ; c ; , and c ; (Eq. 25) respectively. In Fig. 2.c, the black regions are areas where the blur cannot be computed. As we can see in Figs. 2.d and 2.e, the use of c ; and c ; is recommended in the case of vertical and horizontal features respectively. The integration model of the three blur images (Figs., 2.c, 2.d, and 2.e) is given in Fig. 2.f. There are few black regions in this image. This means that the estimation of blur is more dense. Figs. 2.g and 2.h present the blur of Figs. 2.c and 2.f as a function of x and y. It is clear that combination of the blurs estimated using equations c ; , c ; , and c ; leads to more dense and accurate estimation of the blur. The algorithm is fast and allows parallel implementation. We have only tested a particular case of our algorithm (c ; , c ; , and c ; in (Eq. 25)). In further work, we will test other rules for blur computation. 0 0
0 1
1 0
0 1
1 0
0 0
0 1
1 0
0 0
0 1
1 0
Acknowledgments
This work is partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) and le Fonds pour la Formation de Chercheurs et l'Aide a la Recherche (FCAR) du Quebec.
References
[1] J. Ens and P. Lawrence. An Investigation of Methods for Determining Depth from Focus. ieee Trans. on PAMI, 15:97{108, 1993. [2] A. Horii. Depth from Defocusing. Tech. Report ISRN KTH/NA/P{92/16{SE, Royal Institute of Technology, Sweden, 1992. [3] S.H. Lai, C.W. Fu, and S. Chang. A Generalized Depth Estimation Algorithm with a Single Image. ieee Trans. on PAMI, 14:405{411, 1992. [4] J.B. Martens. The Hermite Transform - Applications. ieee Trans. on Acous. Sig. and Speech Proc., 38(9):1607{1618, 1990.
[5] J.B. Martens. The Hermite Transform - Theory. ieee Trans. on Acous. Sig. and Speech Proc., 38(9):1595{1606, 1990. [6] A. P. Pentland. Depth of Scene from Depth of Field. In Proc. Image Understanding Workshop, pp: 253{259, 1982. [7] A.P. Pentland. A New Sense for Depth of Field. ieee Trans. on PAMI, 9:523{531, 1987. [8] A.P. Pentland, S. Scherock, T. Darrell, and B. Girod. Single Range Camera Based on Focal Error. J. Opt. Soc. Am. A, 11:2925{2934, 1994. [9] G. Sansone. Orthogonal Functions. Interscience Publishers, Inc., New York, 1959. [10] M. Subbarao and N. Gurumoorthy. Depth Recovery from Blurred Edges. In Proceedings of IEEE, Int. Conf. on CVPR, pp: 498{503, 1988. [11] M. Subbarao and G. Surya. Depth from Defocus: A Spatial Domain Approach. The Int. Jour. of Comp. Vision, 13(3):271{294, 1994. [12] R.G. White and R.A. Schowengerdt. E ect of Point-Spread Functions on Precision Edge Measurement. J. Opt. Soc. Am. A, 11:2593{2603, 1994. [13] Y. Xiong and S. A. Shafer. Depth from Focusing and Defocusing. In Proceedings of IEEE, Int. Conf. on CVPR, pp: 68{73, 1993. [14] Y. Xiong and S. A. Shafer. Moment and Hypergeometric Filters for High Precision Computation of Focus, Stereo and Optical Flow. Tech. Report CMU-RI-TR-94-28, Carnegie Mellon University, 1994. [15] D. Ziou. Passive Depth From Defocus Using a Spatial Domain Approach. Tech. Report, DMI, Universite de Sherbrooke, 1997.
Lens
Object
Focal
Image kσ
F
u z
v
Figure 1: Image formation in a thin lens.
a
b
c
d
e
f
Blur estimated using c00
Blur estimated using integration model
250
250
200
200
150 Blur
150 Blur
100
100
50
50
0 230
0 230
240
240
250 x
250 x
260
260
270
270
140
150
160 y
170
180 140
g
150
160 y
170
180
h
Figure 2: An example for real images obtained by varying camera parameters. a) and b) Toy house. c), d), and e) describe the blur obtained by equations c ; , c ; , and c ; (Eq. 25) respectively. f) Output of the integration model (Eq. 26). g) and h) The blur of c) and f) within a window, as a function of x and y. The window includes a part of vertical corner of the toy house. 0 0
0 1
1 0
Abstract
This paper presents an algorithm for a dense computation of the di erence in blur between two images. The two images are acquired by varying the intrinsic parameters of the camera. The image formation system is assumed to be passive. Estimation of depth from the blur di erence is straightforward. The algorithm is based on a local image decomposition technique using the Hermite polynomial basis. We show that any coef cient of the Hermite polynomial computed using the more blurred image is a function of the partial derivatives of the other image and the blur di erence. Hence, the blur di erence can be computed by resolving a system of equations. All computations required are local and carried out in the spatial domain. An algorithm is presented for estimation of the blur in 1D and 2D cases and its behavior is studied for constant images, step edges, line edges and junctions. The algorithm is tested using synthetic and real images. The results obtained are very encouraging.
1 Introduction
In this paper, we are interested in the depth from defocusing approach, which involves calculating the distance between a visible surface of an object and the thin lens of the optical camera by using the degree of blurring. Speci cally, let us consider a camera which generates one or more grey-level images of a static 3D scene. Let us assume the perspective projection of a real world point onto the image. The image intensity at a given pixel (x; y) depends on the distance z between the camera lens and the corresponding point in the scene (Fig. 1). In the case of the thin lens, the image is formed by convolution of the ideal projected image and the camera point spread function (PSF) g(x; y; (x; y)), where (x; y) is the blur parameter at the pixel (x; y). We will consider that (x; y) is constant within a given window, that is (x; y) = . Usually this assumption is not realistic. However, without this assumption we are confronted with a non-shift invariant system. These systems are mathematically
hard to handle. This explains why this assumption has been implicit in most prior work. Many functions are used to approximate the PSF of the camera. The reader can nd more information in [12]. We are interested in the most common case, where PSF is approximated by the Gaussian function: x y (1) g (x; y) = 1 e?
2+ 2 2
2
The relationship between the distance z (Fig. 1) and the blur is given by [7]: ( Fv if z > u (2) z = v?FFv?kf if z < u v?F kf +
where u is the distance between the lens and the position of perfect focus, v the distance between the lens and the image plane, F the focal length, f the f-number, and k the proportionality coecient between the blur circle radius and . This formula shows that estimation of depth z requires the estimation of f; F; v; k, and . The intrinsic parameters of the camera f , F , v, and k are independent of pixel location (x; y). They can be estimated by calibration [8, 13, 2]. The blur at a pixel (x; y) is measured by calculating the PSF function at this pixel. The work of Pentland [6] was the rst research in computer vision to test the use of defocus information as a depth cue. Since Pentland's work, many other schemes for the computation of blur have been proposed. The most common approaches are inverse ltering [6, 7, 2, 1, 13, 14], depth recovery from edges [7, 3, 10], and the S-Transform [11]. Later, we will summarize the idea behind the S-Transform approach because it is directly related to our work. In this paper, we will present an algorithm for computing the di erence in blur between two images. These two images are acquired by varying the intrinsic parameters of the camera. The image formation system is assumed to be passive. Our algorithm involves tting the images by a Hermite polynomial basis. We will show that any coecient of the Hermite
polynomial computed using the more blurred image is a function of the partial derivatives of the other image and the blur di erence. Hence, the blur di erence can be computed by resolving a system of equations. All computations required are local and are carried out in the spatial domain. Furthermore, no image model has been used to derive the algorithm. Consequently, our algorithm yields a dense estimation of blur. In the next section, we will start by describing the blur estimation rule for 1D images. Then we will present the extension of the 1D mathematical model to 2D images. The algorithm for blur estimation in 2D images is described in section 4. Section 5 presents the behavior of the algorithm and the experimental results.
2 Blur Estimation in a 1D Signal
In this section we present the mathematical development behind blur estimation in the case of a 1D signal. Let us consider Ic (x) and Ib (x), two unfocused images of the same scene obtained by varying one or more camera parameters (see Fig. 1). Let us assume that Ib (x) is more blurred than Ic (x). The image which has the smaller grey-level variance is the more blurred one. When the PSF is Gaussian, the relationship between these images is:
Ib (x) = (Ic g )(x)
(3)
where g (x) is a Gaussian of variance and p represents convolution. The parameter = b ? c is the blur di erence between the two images, where b and c are the blur of Ib (x) and Ic (x) respectively. The relationship between and z can be deduced. For example, let us consider the image Ib (x) obtained by changing the parameter v by an amount v in (Eq. 2), which causes a change of magnitude in the blur . The depth z is given by: 2
2
2
8 > ? 1 =< F z > : F? 1
v v (1 + 1
2 +
q
2
2
1 + f F vvv k ) if z > F
q
2 (2 + 2
)
2
2
f v v v v (1 ? 1 + F v k ) if z < F (4) For more details see [15]. This equation shows that the depth can be computed from the intrinsic camera parameters and the blur di erence . It should be recalled that the camera parameters can be obtained by calibration. Now, we are faced with the problem of estimating . To this end, we have developed an image decomposition technique, called the Hermite polynomial transform, in which the image is locally approximated by Hermite polynomials. The scheme involves tting the image, within the Gaussian window g (x) 1
1
2 +
2 (2 + 2
)
2
2
in (Eq. 1), by a Hermite polynomial basis. The nth Hermite polynomial Pn (x) is de ned by:
n x x (5) Pn (x) = (?1)n e d(dx )n e? The sequence fPn (x)g is orthonormal with respect to g (x) on (?1; +1). Our choice of the Hermite poly2 2
2 2
nomial is due to the fact that these polynomials are associated with the Gaussian function. Hence, they have interesting properties and their use in image processing is increasing [5, 4]. In this mathematical model, the complete description of an image requires that the tting process be repeated at a sucient number of pixels within the window. We will consider the case of equidistant sample spacing T . The approximation of image Ib (x) within this window by the polynomial basis fPn (x)g is: 1 X K X
Ib (x) =
n k ?K =0
cn (kT )Pn (x ? kT )
(6)
=
where 2K + 1 is the maximum number of samples in the window and cn (kT ) is de ned by: cn (kT ) = (Ib hn )(kT ) (7) where n x hn (x) = Pn (x)g (x) = (?1)n d(dx )n p1 e? (8) Substituting (Eq. 3) in (Eq. 7), we obtain: cn (kT ) = (Ic g hn )(kT ) = (Ic fn )(kT ) (9) where n x f (x) = (g h )(x) = (?1)n n d p1 e? (10) 2 2
n
n
dxn
2 2
and = + . (Eqs. 7 and 9) show the relationship between image Ib and image Ic . We can use these two equations to estimate the blur. However, we are faced with solving a non-linear equation. Now we want to give a simple form of the relationship between these two images which will make the computation rule of explicit. Let us consider the Taylor expansion of Ic (kT ? x) at kT : 2
2
2
Ic (kT ? x) =
1 (?x)p I p (kT ) X c +
p
( )
=0
p!
(11)
where Icp (x) denotes the pth derivative of Ic (x). Substituting this equation in (Eq. 9), we obtain: ( )
cn (kT ) =
X1 +
p
=0
sn;p Icp (kT ) ( )
(12)
where
Z 1 (?x)p f (x) n dx sn;p = p!
(13)
?1
(Eqs. 7 and 12) show the relationship between image Ib and image Ic . In this relation, a coecient sn;p is proportional to the pth moment of fn (x) and is a function of the blur . In other words, given Ic (x), Ib (x) and fn (x) we are able to estimate . However, in momentbased schemes, calculating the moment directly entails a high computation cost. To overcome this diculty, we propose an ecient computation procedure for sn;p [15]:
8 > < sn;p = > :
p?n
(?1)nn p?n?p?np?n?p?n= if 0 < n < p and p ? n is even (?1)nn if n = p 0 otherwise (
1)(
3)
)!2(
(
3
(14)
This equation is an ecient computation rule for
By using (Eqs. 7, 12, and 14), the relation between the two images Ib (x) and Ic (x) is now given by:
8 P 1 s I p (kT ) = (I h )(kT ) > < c r (kT ) = P b r pr r; p c p 1 > : (cIbr h(kTr ) )(=kT )pr s r ; p Ic (kT ) = (2 )
+
2
2
2
(2 +1)
+
2 +1
2 +1 2 +1
2 +1
(15) We are now able to compute the blur di erence. For the sake of clarity, let us examine the use of coecients c and c in (Eqs. 15) Pas an example. In the case of c , we assume that p> s ; p Ic p (kT ) = 0 and P for c , we assume that p>1 s ; p Ic p (kT ) = 0. According to (Eq. 14), s ; = 1, s ; = ?, s ; = , and s ; = ? . Thus (Eq. 15) becomes: 0
1
(2 )
0
1
0 2
1
1 2 +1
(2 +1)
+
1
0 0
2
(
1 3
2
1 1
0 2
4
4
c (kT ) = Ic (kT0 ) + Ic00 (kT000) = (Ib h )(kT ) c (kT ) = ?Ic (kT ) ? Ic (kT ) = (Ib h )(kT ) 2
0
4
2
1
0
(16) where h (kT ) and h (kT ) are given in (Eq. 8). By resolving the two equations in (16) we obtain respectively: 1
s
= 4 (Ib h )(I 00kT(kT) ?) Ic (kT ) ? c 0
or
s
Ic0 (kT ) ? = 4 (Ib h ?)(IkT000)(+ kT ) 1
c
0
2
0
2
(17)
(18)
+
=0
cn;m (kTx; lTy ) =
=0
=
1X 1 X +
+
p
=0
q
=
sn;p;m;q Icp q (kTx ; lTy ) ( )( )
=0
(21) where sn;p;m;q is de ned by: Z 1 Z 1 (?x)p (?y)q f (x; y) n;m sn;p;m;q = dxdy p ! q ! ?1 ?1 (22) fRn;m(x; y)p is separable, thus s = s s = n;p;m;q n;p m;q 1 ?x fn x dx R 1 ?y q fm y dy. Rapid calculap q ?1 ?1 tion of sn;p;m;q is possible by using (Eq. 14). The coecient cn;m (kTx; lTy ) in (Eq. 21) becomes: +
(
)
( )
(
)
!
cn;m (kTx; lTy ) =
1
4
0
The extension of the 1D scheme for blur estimation to the 2D case is straightforward. In fact, a 2D Hermite polynomial can be created by taking the tensor product of two sets of 1D Hermite polynomials; i.e, fP (x)P (y); ; Pn (x)Pm (y); g, where n is the degree of x and m the degree of y. Since a Gaussian window g (x; y) is also separable, it is easy to show that fPn;m(x; y)g is orthogonal within g (x; y). Let us consider that kTx and lTy are the sample spacings in the window along the x and y axes. The decomposition of image Ib (x; y) into localized polynomials is given by: P P P P Ib (x; y) = n 1 m1 Kk ?K Ll ?L cn;m (kTx; lTy ) Pn;m (x ? kTx ; y ? lTy ) (19) where (2K + 1) (2L + 1) is the size of the sampling lattice and cn;m(kTx ; lTy ) is de ned by: cn;m (kTx; lTy ) = (Ib hn;m)(kTx ; lTy ) (20) where hn;m (x; y) = Pn;m (x; y)g (x; y). Both Pn;m (x; y) and g (x; y) are separable; thus straightforward manipulations, as in the 1D case, give: +
) 2
sn;p . The simpli cation of sn;p also simpli es the computation of cn (kT ) and the approximation of Ib (x).
2
3 Extension to a 2D Image
( )
!
1X 1 X +
+
pn qm
sn;p;m;q Icp q (kTx ; lTy ) ( )( )
which can be simpli ed to obtain:
8 > > > > > > < > > > > > > :
P 1P +
(23)
1
+
c2r;2t (kTx ; lTy ) = pr qt s2r;2p s2t;2q (2p)(2q ) (kTx ; lTy ) Ic +1 +1 c2r;2t+1 (kTx ; lTy ) = pr qt s2r;2p s2t+1;2q+1 (2p)(2q +1) (kTx ; lTy ) Ic +1 +1 c2r+1;2t (kTx ; lTy ) = pr qt s2r+1;2p+1 s2t;2q (2p+1)(2q ) (kTx ; lTy ) Ic +1 +1 c2r+1;2t+1 (kTx ; lTy ) = pr qt s2r+1;2p+1 (2p+1)(2q +1) (kTx ; lTy ) s2t+1;2q+1 Ic
P P
P P
P P
(24)
derivatives Icp q (x; y) can be estimated using di erentiation techniques such as Sobel masks, Gaussian masks, wavelets, polynomial tting, etc. The choice of algorithm parameters Tx , Ty , N , M , P and Q is discussed in [15]. ( )( )
We are now able to compute the blur di erence between two 2D images. For clarity, let us take as an example the coecients c ; , c ; , and c ; in (Eqs. 24). Let us assume that the image is a cubic polynomial; that is Icp q (kTx ; lTy ) = 0; 8p + q > 3. By using (Eqs. 14 and 20), we can rewrite the rst three equations in (24): 0 0
0 1
1 0
( )( )
5 Performance Evaluation
5.1 Algorithm behavior
8 c ; (kTx ; lTy ) = Ic (kTx; lTy ) + rIc (kTx; lTy ) = > > (Ib h ; )(kTx ; lTy ) > > < c ; (kTx; lTy) = ?Ic (kTx; lTy )? > Ic (kTx ; lTy ) = (Ib h ; )(kTx ; lTy ) > > (kTx; lTy )? x ; lTy ) = ?Ic > : c ; (kT Ic (kTx ; lTy ) = (Ib h ; )(kTx ; lTy ) (25) where r is the Laplacian. The rst equation in (25) 2
0 0
4
0 0
(0)(1)
0 1
2
(2)(1)
0 1
4
1 0
(1)(0)
2
(1)(2)
1 0
4
0
is the same as the one implemented by M. Subbarao and G. Surya [11]. These authors propose a spatial approach, called the S-Transform. They assume that the by the cubic polynomial p(x; y) = P image P ?ism tted a x n;m m y n . In this case, they show that m n Ic (x; y) = Ib (x; y) ? rIb (x; y)=4, where r is the Laplacian and the blur di erence. The estimation of , and hence of the depth, is straightforward. It is easy to see that their model for blur estimation is a particular case of our blur estimation model presented in (Eq. 15). 3
=0
0 0
2
2
( )( )
2
4 Algorithm for Blur Estimation
( )( )
We have presented the mathematical developments required for the computation of blur. Our next step is to specify an algorithm for the estimation of the blur di erence , based on these mathematical developments. Let us consider a nite summation in (Eqs. 20 and 24) for a given n and m; that is p 2 [r; P ] and q 2 [t; Q]. The blur is estimated by using (Eq. 23) and any equation in (24). However, there are four equations in (24); thus we estimate by using (23) and each equation in (24). In each window W centered at (x; y), we obtain four values of , from which we choose the one which minimizes the following quadratic error: e2 (x; y ) =
XX r
s
(Ib (x + r; y + s)
( )( )
? (Ic g )(x + r; y + s))
2
(26)
To summarize, given , Tx , Ty , N , M , P , and Q, the algorithm involves three steps: 1) the estimation of cn;m using (Eq. 23); 2) the estimation of four values of using (Eq. 24); 3) the selection of the best according to the quadratic error in (Eq. 26). The image
0
0
3
=0
Up to now, no assumption has been formulated about the images for which our algorithm performs well. In this section, we will examine its behavior in the case constant images, step edges, line edges, and junctions. 1) Constant images: the derivatives of the image are zero and the system of equations in (23) and (24) is reduced to sn; ;m; Ic (kTx; lTy ) = (Ib hn;m)(kTx ; lTy ). According to (Eq. 14), sn; ;m; is one if n = m = 0 and zero elsewhere, thus Ic (kTx; lTy ) = (Ib h ; )(kTx ; lTy ). This equation does not permit computation of blur. In other words, low frequencies are not used by our algorithm to estimate blur. 2) Step edges: we consider the step edge model, I (x; y) R=x ((x ? kTx )cos() + (y ? lTy )sin()), where (x) = ?1 exp(?t )dt and is the edge orientation. I p q (kTx; lTy ) = 0, if p or q is even. In this case, only the last equation in (24) can be used to estimate blur. 3) Line edges: we consider the line model, I (x; y) = e?((x?kTx )cos()+(y?lTy )sin()) , where is the line orientation. I p q (kTx ; lTy ) = 0, if p or q is odd. In this case, only the rst equation in (24) can be used to estimate blur. 4) Junctions: for the L junction model, I (x; y) = (x ? kTx ) (y ? kTy ). I p q (kTx ; lTy ) = 0, if p or q is even. In this case, only the last equation in (24) can be used to estimate blur. This algorithm exhibits the same behavior in the case of a T junction, I (x; y) = (x ? kTx )( (y ? kTy )+1). Let us examine the behavior of the rst equation in (Eq. 25) implemented by Subbarao and Surya, where the image is t by a polynomial basis of degree less than or equal to three. At a step edge, the Laplacian is zero and therefore this equation becomes c ; (kTx; lTy ) = Ic (kTx; lTy ). Thus, it is not possible to compute the blur at a step edge. At a line edge, the Laplacian is not zero and it is possible to compute the blur. For linear junction models (the intersection point of n sectors with constant intensities such as L, V, T, Y, and X), the Laplacian of Gaussian is zero at the junction point. Hence, it is not possible to estimate the blur at the junction point. We conclude that Subbarao and Surya's algorithm cannot be used to estimate blur at step edges and junctions. 0 0
0
5.2 Experimental Results
We have tested the equations in (25) with two real images (Fig. 2.a and 2.b) of a toy house taken using an optical camera with di erent aperture and exposure time settings. The images have 384 320 pixels and 256 grey levels. They were given to us by Y. Xiong of the Robotics Institute at Carnegie Mellon University. The toy house is about two meters away from the camera. It is about 2.5 inches wide, 3 inches long, and 4.5 inches high. The partial derivatives of the image are computed by convolving the image with the appropriate derivatives of the Gaussian function. The scale is 1, the spacings Tx = Ty = 1. Figures 2.c, 2.d, and 2.e present the blur estimated using c ; c ; , and c ; (Eq. 25) respectively. In Fig. 2.c, the black regions are areas where the blur cannot be computed. As we can see in Figs. 2.d and 2.e, the use of c ; and c ; is recommended in the case of vertical and horizontal features respectively. The integration model of the three blur images (Figs., 2.c, 2.d, and 2.e) is given in Fig. 2.f. There are few black regions in this image. This means that the estimation of blur is more dense. Figs. 2.g and 2.h present the blur of Figs. 2.c and 2.f as a function of x and y. It is clear that combination of the blurs estimated using equations c ; , c ; , and c ; leads to more dense and accurate estimation of the blur. The algorithm is fast and allows parallel implementation. We have only tested a particular case of our algorithm (c ; , c ; , and c ; in (Eq. 25)). In further work, we will test other rules for blur computation. 0 0
0 1
1 0
0 1
1 0
0 0
0 1
1 0
0 0
0 1
1 0
Acknowledgments
This work is partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) and le Fonds pour la Formation de Chercheurs et l'Aide a la Recherche (FCAR) du Quebec.
References
[1] J. Ens and P. Lawrence. An Investigation of Methods for Determining Depth from Focus. ieee Trans. on PAMI, 15:97{108, 1993. [2] A. Horii. Depth from Defocusing. Tech. Report ISRN KTH/NA/P{92/16{SE, Royal Institute of Technology, Sweden, 1992. [3] S.H. Lai, C.W. Fu, and S. Chang. A Generalized Depth Estimation Algorithm with a Single Image. ieee Trans. on PAMI, 14:405{411, 1992. [4] J.B. Martens. The Hermite Transform - Applications. ieee Trans. on Acous. Sig. and Speech Proc., 38(9):1607{1618, 1990.
[5] J.B. Martens. The Hermite Transform - Theory. ieee Trans. on Acous. Sig. and Speech Proc., 38(9):1595{1606, 1990. [6] A. P. Pentland. Depth of Scene from Depth of Field. In Proc. Image Understanding Workshop, pp: 253{259, 1982. [7] A.P. Pentland. A New Sense for Depth of Field. ieee Trans. on PAMI, 9:523{531, 1987. [8] A.P. Pentland, S. Scherock, T. Darrell, and B. Girod. Single Range Camera Based on Focal Error. J. Opt. Soc. Am. A, 11:2925{2934, 1994. [9] G. Sansone. Orthogonal Functions. Interscience Publishers, Inc., New York, 1959. [10] M. Subbarao and N. Gurumoorthy. Depth Recovery from Blurred Edges. In Proceedings of IEEE, Int. Conf. on CVPR, pp: 498{503, 1988. [11] M. Subbarao and G. Surya. Depth from Defocus: A Spatial Domain Approach. The Int. Jour. of Comp. Vision, 13(3):271{294, 1994. [12] R.G. White and R.A. Schowengerdt. E ect of Point-Spread Functions on Precision Edge Measurement. J. Opt. Soc. Am. A, 11:2593{2603, 1994. [13] Y. Xiong and S. A. Shafer. Depth from Focusing and Defocusing. In Proceedings of IEEE, Int. Conf. on CVPR, pp: 68{73, 1993. [14] Y. Xiong and S. A. Shafer. Moment and Hypergeometric Filters for High Precision Computation of Focus, Stereo and Optical Flow. Tech. Report CMU-RI-TR-94-28, Carnegie Mellon University, 1994. [15] D. Ziou. Passive Depth From Defocus Using a Spatial Domain Approach. Tech. Report, DMI, Universite de Sherbrooke, 1997.
Lens
Object
Focal
Image kσ
F
u z
v
Figure 1: Image formation in a thin lens.
a
b
c
d
e
f
Blur estimated using c00
Blur estimated using integration model
250
250
200
200
150 Blur
150 Blur
100
100
50
50
0 230
0 230
240
240
250 x
250 x
260
260
270
270
140
150
160 y
170
180 140
g
150
160 y
170
180
h
Figure 2: An example for real images obtained by varying camera parameters. a) and b) Toy house. c), d), and e) describe the blur obtained by equations c ; , c ; , and c ; (Eq. 25) respectively. f) Output of the integration model (Eq. 26). g) and h) The blur of c) and f) within a window, as a function of x and y. The window includes a part of vertical corner of the toy house. 0 0
0 1
1 0
