Surface shape estimation from photometric images

ARTICLE IN PRESS

Optics and Lasers in Engineering 42 (2004) 461–468

Surface shape estimation from photometric images ! Villaa, Juan B. Hurtado-Ramosb,* Jesus a

Laboratorio de Procesamiento Digital de Sen˜ales, Unidad de Ingenier!ıa El!ectrica, ! Universidad Autonoma de Zacatecas, C.P. 98000. Zacatecas, Mexico b Centro de Ingenier!ıa y desarrollo Industrial, Av. Playa Pie de la Cuesta #702, Desarrollo Habitacional San Pablo, C.P. 76130, Quer!etaro, Qro., Mexico Received 16 June 2003; accepted 30 December 2003

Abstract One of the most studied techniques for recovering surface shapes using a computer vision system is the photometric. This method is based on the analysis of one or several images of an object illuminated from a known direction . This kind of images can be considered reflectance maps that give us information of the surface gradient. In computer vision, the problem of recovering 3D information from shaded images is considered an inverse problem. To integrate surface gradient information it is proposed a regularization technique that gives a stable solution of the inverse problem and allows the possibility of reducing errors caused by noise. Results applied on synthetic and real experimental images are presented. r 2004 Elsevier Ltd. All rights reserved. Keywords: Shape from shading; Regularization

1. Introduction Extraction of 3D information from 2D images is a fundamental problem in computer vision. There are different methodologies to realize this task. Among the most known methods is the so-called shape from shading or photometric [1–3]. This technique is based on three main factors: the shape of the object, the reflectance properties of the surface and the illumination, i.e. the light source position. Because physical parameters of reflectance and illumination direction are known, the shape *Corresponding author. E-mail address: [email protected] (J.B. Hurtado-Ramos). 0143-8166/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2003.12.004

ARTICLE IN PRESS 462

! Villa, J.B. Hurtado-Ramos / Optics and Lasers in Engineering 42 (2004) 461–468 Jesus

of the object can be obtained considering that the brightness at every point in the image depends on the surface reflectance and orientation at that point. A number of shape from shading techniques have been developed in the last three decades. Some of them use a single image for recovering 3D information [4,5], however, good results with these methods are restricted to simple objects and ideal experimental conditions, mainly because of the lack of a robust solution methodology. Other researchers have developed techniques using multiple images for more reliable results. Works have also been made using the photometric stereo method for deformation detection [6]. In general, the application of any method mentioned above does not give the height distribution itself but the gradient of the surface. For obtaining the height distribution an integration method is required. However, the gradient integration has a key problem: the orientations are not always consistent, which makes it difficult to define an appropriate relation for slope-to-height conversion. This problem, called the integrability problem [2], can be analyzed from a different perspective. Remember that it is an inverse problem, where it is required to recover 3D information from 2D arrays (i.e., images). Essentially, this inverse problem is considered ill-posed, because observations (i.e., the slopes) does not determine the value of the height distribution in a unique an stable way [7] due to the loss of information during the imaging process that projects a three-dimensional world into two-dimensional images. Lee and Kuo [3] reported work that alleviates in part the integrability problem by using a quadratic cost functional with nodal basis-functions for representation of the surface, however, the method is not robust enough in the presence of noise or in the case in which data is missing or is unreliable because it does not takes into account prior information of the solution. Existence of an ill-posed problem suggests a regularization method to solve it. The key idea to solve an ill-posed problem is to restrict the possible solutions by introducing prior knowledge. In this work, it is proposed a regularization approach to solve the integration problem in the shape from shading technique with three images. This method provides suitable results and allows the possibility of reduce errors caused by noise.

2. Shape from shading with three images Assuming an object whose surface is Lambertian and is illuminated from three different directions (light source is supposed to be sufficiently away from the object in such a way that rays of light illuminating it can be considered parallel). Three different images are then recorded each of them having a different illumination direction. The intensities of the three recorded images are functions proportional to the reflectance maps that provide information of the relationship between the shape gradient and the irradiance. The three images can be represented as [1] Ek ¼ rðs#k  n# Þ;

k ¼ 1; 2; 3;

ð1Þ

ARTICLE IN PRESS ! Villa, J.B. Hurtado-Ramos / Optics and Lasers in Engineering 42 (2004) 461–468 Jesus

463

where r is the albedo of the surface and ðpk ; qk ; 1ÞT s#k ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 þ p2k þ q2k

k ¼ 1; 2; 3;

ð2Þ

ðp; q; 1ÞT n# ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; 1 þ p2 þ q2

ð3Þ

where n# denotes the unit vector of the gradient of the height z (i.e., p ¼ @z=@x; q ¼ @z=@y) for each object’s surface point and s#k the unit vector defining the slope of an imaginary plane whose normal is the illumination source’s direction and is used to define the light source position. Combining Eqs. (1–3) we rewrite E ¼ rSn# ; where the matrix S contains the unit vectors s#k ; and components of vector E are the three intensities Ek : Thus, assuming that S is not singular, the gradient of the surface is obtained with rn# ¼ S1 E; where

ð4Þ

2

3 s# 1 6 7 S ¼ 4 s# 2 5:

ð5Þ

s# 3

3. Integration of the gradient As explained before, gradient integration is an ill-posed problem that must be regularized to obtain proper results. Regularization of linear problems has been strongly studied. Tikhonov [8] analyzed the solution of the problem of finding z from any given data set g (e.g. gradient data) which are related in the following way: g ¼ Az;

ð6Þ

where A is a noninjective linear operator. In that work it is proposed that the regularized solution can be obtained by finding the value of z that minimizes the functional UðzÞ ¼ jjAz  gjj2 þ mjjPzjj2 ;

ð7Þ

where m is the so-called regularization parameter, and P is the regularization operator that is formed by a linear combination of derivatives. Following the same reasoning on (7), a regularization procedure for obtaining the height distribution in the shape from shading problem can be applied considering the operator A as a first order difference, so that X X jjAz  gjj2 ¼ ½zi;j  zi1;j  pi;j 2 þ ½zi;j  zi;j1  qi;j 2 ; ð8Þ ði;jÞAL

ði;jÞAL

ARTICLE IN PRESS ! Villa, J.B. Hurtado-Ramos / Optics and Lasers in Engineering 42 (2004) 461–468 Jesus

464

where (i,j) represents the coordinates in the image, zi,j the height distribution to be estimated and L a finite regular lattice where the observations are available. As can be seen, for the case m ¼ 0 this problem corresponds to the least-squares procedure of integration [9,10] and height distribution could be obtained minimizing this cost functional with respect to zi,j, however, the least-squares procedure does not always provide suitable results because it assumes a noiseless information and the observations ðp; qÞ does not determine the value of the field zi,j in an unique and stable way, mainly because the algorithm does not include any information about the solution, e.g., smoothness of the estimated signal. In the proposed regularization procedure it is assumed ma0 which restricts the solution to have specific characteristics, i.e., to be smooth. If it is assumed that the operator P is a finite approximation to the second-order partial derivatives, the second term of Eq. (7) becomes X X mjjPzjj2 ¼ m ½ziþ1;j  2zi;j þ zi1;j 2 þ m ½zi;jþ1  2zi;j þ zi;j1 2 : ð9Þ ði;jÞAL

ði;jÞAL

In this way the solution is restricted to have smooth derivatives. Parameter m controls the smoothness of the field to be estimated and its value may depend on the amount of noise. There is not a precise procedure to determine the value of this parameter in a given problem so that it is selected intuitively [11]. Finally, the estimated field zi,j is the minimizer of the cost function (7), so that taking the gradient of U(z) with respect to zi,j and equating it to zero: @UðzÞ ¼ 2ðDULS þ DUR Þ ¼ 0; @zi;j

ð10Þ

where DULS ¼ ½zi;j  zi1;j  pi;j  ½ziþ1;j  zi;j  piþ1;j þ ½zi;j  zi;j1  qi;j  ½zi;jþ1  zi;j  pi;jþ1

ð11Þ

DUR ¼  m½6zi;j  4ðziþ1;j þ zi1;j Þ þ ziþ2;j þ zi2;j  m½6zi;j þ 4ðzi;jþ1 þ zi;j1 Þ þ zi;jþ2 þ zi;j2 :

ð12Þ

and

One can obtain the minimizer zi,j solving the linear system (10) by applying any iterative algorithm from the literature [12]. The minimizer of Eq. (12) smoothes out the observations, hence, the minimization operation can also be considered a low-pass filtering operation. To analyze the frequency response of the minimization it is taken the discrete Fourier transform of Eq. (12) only in the i direction (i.e., considering only the slope p for simplicity) obtaining the transfer function that is expressed as HðoÞ ¼

ZðoÞ 1  expðioÞ ¼ ; PðoÞ 2  2 cosðoÞ þ m½8 cosðoÞ  2 cosð2oÞ  6

ð13Þ

which represents a low-pass filter whose cut-frequency is controlled by parameter m: A graphic of the normalized response for different values of m is shown in Fig. 1. As

ARTICLE IN PRESS ! Villa, J.B. Hurtado-Ramos / Optics and Lasers in Engineering 42 (2004) 461–468 Jesus

465

Fig. 1. Frequency response of the minimization operation for different values of m:

can be seen the analytic response has an asymptote for o ¼ 0 which is an expected result for p=0. On the other hand, in the particular case of least squares (i.e, m ¼ 0) the transfer function is represented by 1  expðioÞ : ð14Þ HðoÞ ¼ 2  2 cosðoÞ 4. Experimental results Performance of the proposed technique was tested in two experiments. The first one was a numerical estimation with synthetic images of reflectance maps using the simulated sphere shown in the contour map of Fig. 2(a). The images were corrupted with –5 dB of uniformly distributed noise then, the gradient (p,q) was calculated using Eq. (4). The height distribution obtained applying least-squares (m ¼ 0) is shown in Fig. 2(b). Fig. 2(c) shows the results applying the proposed technique with m ¼ 3: In a second experiment an estimation of the height distribution of a pottery object was made. In the experiment the values p1=0.36, p2=0, p3=0.36 and q1=0.05, q2=0.36, q3=0.05 for the positions of the light source were used. The resulting images of 250 250 pixels are shown in Fig. 3(a)–(c). The contour map of the height distribution applying the proposed technique is shown in Fig. 3(d).

5. Conclusions Recovering of height distribution from gradients using the shape from shading technique is an ill-posed problem which requires the introduction of prior

ARTICLE IN PRESS 466

! Villa, J.B. Hurtado-Ramos / Optics and Lasers in Engineering 42 (2004) 461–468 Jesus

Fig. 2. (a) Contour map of the height distribution of a simulated sphere. Estimated shape using (b) the simple least-squares (m ¼ 0) and (c) the proposed method for m ¼ 3:

ARTICLE IN PRESS ! Villa, J.B. Hurtado-Ramos / Optics and Lasers in Engineering 42 (2004) 461–468 Jesus

467

Fig. 3. Shaded images of a pottery object using the source positions described by the slopes (a) p1=0.36, q1=0.05; (b) p2=0, q2=0.36; (c) p3=0.36, q3=0.05. (d) contour map of the estimated height distribution (all dimensions are in pixels).

information about the behavior of the solution to solve it, i.e., regularizing the solution. In this work it was presented a procedure for gradient integration in shape from shading technique with three images. As shown, this inverse problem can be easily solved by using a regularization technique that improves the estimations obtained in the particular case of the least-squares method. It was explained and shown that the technique allows the possibility of reduce errors caused by noise. Results of the proposed algorithm were presented with simulated and real experiments.

Acknowledgements Authors would like to acknowledge the financial support of the Centro de Ingenier!ıa y Desarrollo Industrial (Quere! taro, Me! xico) to complete the present work.

ARTICLE IN PRESS 468

! Villa, J.B. Hurtado-Ramos / Optics and Lasers in Engineering 42 (2004) 461–468 Jesus

References [1] Horn BKP. Robot vision. Cambridge, MA: MIT Press; 1986 [Chapters 10,11]. [2] Horn BKP. Height and gradient from shading. Int J Comput Vision 1990;5:584–95. [3] Lee KM, Kuo J. Surface reconstruction from photometric stereo images. J Opt Soc Am 1993;10: 855–68. [4] Zheng Q, Chellapa R. Estimation of illumination direction, albedo, and shape from shading. IEEE Trans Pattern Anal Mach Intell 1991;13:680–702. [5] Szeliski R. Fast surface interpolation using hierarchical basis functions. IEEE Trans Pattern Anal Mach Intell 1990;12:513–28. [6] Carver AE, Schalkoff RJ. Deformation detection on composite reflectance surfaces using two image photometric stereo. Opt Eng 2001;40:295–302. [7] Hadamard J. Sur les problems aux derives partielles et leur signification physique. Princeton University Bulletin 13. Princeton, NJ: Princeton University; 1902. [8] Tikhonov AN. Solution of incorrectly formulated problems and the regularization method. Sov Math Dock 1963;4:1035–8. [9] Fried DL. Least-squares fitting a wave front distortion estimate to an array of phase difference measurements. J Opt Soc Am 1977;67:370–5. [10] Hudgin RH. Wave-front reconstruction for compensated imaging. J Opt Soc Am 1977;67:375–8. [11] Bertero M, Boccacci P. Introduction to inverse problems in imaging. London: Institute of Physics; 1998 [Chapter 5]. [12] Gollub GH, Van Loan CF. Matrix computations. Baltimore, MD: Johns Hopkins U. Press; 1990.