- Email: [email protected]
Simultaneous estimation of shape and reflectance map from photometric stereo
CVGIP:
IMAGE
UNDERSTANDING
Vol. 55, No. 3, May, pp. 275-286,
1992
Simultaneous Estimation of Shape and Reflectance Map from Photometric Stereo D. TAGARE*
HEMANT Department
of Electrical
and Computer Engineering,
Rice University,
Houston,
Texas 77005
AND
Rur J. P. DEFIGUEIREDO Department
of Electrical
Engineering
and Department Received
of Mathematics,
of California
isolate
California
92717
the Lambertian
component.
A similar
strategy
is
also used in [8]. Silver [9] considered the problem of simultaneously estimating reflectance map parameters and shape without assuming that the components of the reflectance map can be easily separated. Using his technique, the Lambertian part of the reflected light can be removed, thereby isolating the specular
part.
This technique
has not been tested
experimentally, and Silver himself cautioned about the numerical robustness of the algorithm. In this paper, we consider simultaneous estimation of reflectance map parameters and shape in the context of diffusely reflecting surfaces. Such surfaces occur commonly, and reflected light from them has several components. The geometric distribution of light in these components cannot always be easily decomposed, and the techniques developed in [4, 81 cannot be applied. Therefore, it is necessary to develop alternate simultaneous estimation techniques for such surfaces. We begin this development by proposing a class of reflectance maps for such surfaces and then investigate how the parameters of the map and the object shape can be estimated from photometric stereo. Specifically, we consider the following questions of existence, uniqueness, and robustness: 1. The number of light sources: Before launching into specific techniques for obtaining the solution, it is essen-
INTRODUCTION
Shape from shading and photometric stereo techniques estimate the shape of an object from the intensities in an image using knowledge of the light source positions, the viewing direction, and how the object reflects light. The information about light source positions, viewing direction, and reflectance is mathematically described by the reflectance map. Accurate prior knowledge of the reflec-
tance map is required for precise shape reconstruction. Such knowledge is often not available even for controlled conditions, and under these circumstances it is desirable to estimate the reflectance map along with the shape. * Current address: Department of Diagnostic ment of Electrical Engineering, Yale University, icut 06510.
Irvine,
Algorithms for estimating parameters of reflectance maps which are sums of Lambertian and specular components have been proposed before. Most of these are based on the observation that the geometric distributions of reflected light in Lambertian and specular components are so different that the two can be easily separated. Coleman and Jain [4] report a procedure to estimate the shape of an object and the albedo of the Lambertian component of the reflectance map. Four point light sources are used, and the intensities in the images are analyzed to
Press, Inc.
1.
at Irvine,
May 1, 1990; accepted July 11, 1991
This paper investigates whether the shape of an object and certain parameters of its reflectance map can be simultaneously estimated using photometric stereo. This problem has been addressed in the literature for the case where the Lambertian and non-Lambertian components in the image can be easily separated. No such separability is assumed in this paper. A class of reflectance maps for modeling diffusely reflecting surfaces is proposed. This class is based on the physics of scattering from real world surfaces. Next, the problem of joint estimation of some parameters of the map along with the surface shape is analyzed. A bound is obtained on the number of light sources necessary for a unique solution to the problem. The analysis also reveals that some of the estimates can be obtained by a nonparametric method. The behavior of the estimates in the presence of noise is also investigated. It is shown that simultaneous estimation is ill-posed. Regularizing the estimates yields good reconstructions from real world data. 0 1992 Academic
University
Radiology and DepartNew Haven, Connect-
275
All
1049.9hfXVY2 $5.00 Copyright G 1992 by Academic PI-~\\. Iw. rights of reproduction in any form rewrved.
276
TAGARE AND DeFIGUEIREDO
tial to check that the solution exists and is unique. Therefore, we ask, is there a unique solution to simultaneous estimation and how many light sources are needed to guarantee it? 2. Separation of variables: Can the problem be separated into simpler problems? Even though the original problem is complex and might not exhibit a closed form solution for all the unknowns, perhaps some of the unknowns might be available in a simple manner, while the rest may be obtained numerically. 3. Robustness with respect to noise: How robust is the solution with respect to noise? In Section 2, we review the reflectance map model. In Section 3, we formulate the problem of simultaneous estimation of the reflectance map parameters and surface shape. Section 4 addresses the problem of obtaining a unique solution. Section 5 illustrates the ill-posed nature of simultaneous estimation and Section 6 regularizes it. Section 7 contains experimental results. 2. THE REFLECTANCE MAP
Figure 1 illustrates the imaging geometry. The object is illuminated by a point light source at infinity. Bold face n, i, r denote unit vectors along the surface normal, the incident direction, and the receiving direction, respectively. The reference frame is aligned so that r is along the positive z-axis. The plane spanned by i and r is called the “principal plane” and the angle between i and r the “phase angle.” The zenith and azimuth angles of any unit vector are denoted by 0 and 4. Subscripts denote the vector associated with the angles; e.g., the zenith and azimuth angles of n are denoted by 8, and $,, respectively. The zenith angle is measured positively down from the z axis, and the azimuth angle is measured anticlockwise from the x axis.
Light source m
The intensity in the image, Z(x, y), is given by [5] Z(x, Y) = L,R(i,
0,
Y>, r),
where R( > is the reflectance map. For a point light source, R( ) is given by [5] R(i, n, r) = afi(i, n, r) iTn, for iTn, 2 0 = 0 otherwise.
(2)
where fr( ) is the BRDF (the Bidirectional Reflectance Distribution Function) of the surface and (Yis a positive constant called the “albedo” [5]. The BRDF is defined as the ratio of the radiance of light reflected along r to the irradiance of incident light on the surface.
2.1.
The BRDF of Rough Surfaces
In general, the BRDF of any surface is a sum of delta functions and functions of bounded variation. The delta functions are attributed to Specular Rej7ection and the rest to Diffuse Rejection. As a first step toward obtaining a theoretical expression for the BRDF of a solid, the solid can be modeled as a regular array of particles having a planar boundary at the top. Reflection from the boundary is analyzed by considering how individual particles scatter and how their behavior is modified by their being together in a dense array. The Ewald-Oseen theory describes the reflectance behavior of such a configuration and the reader is referred to [lo] for a qualitative discussion of it. The main result of this theory is that perfect solids can only reflect specularly. In this paper, we ignore specular reflection. Real world solids are not perfect. At the scale of the wavelength of visible light their surfaces are rough and their bulk is inhomogeneous. These imperfections cause light to be diffusely reflected in nonspecular directions. Diffuse reflection from most engineering materials is contained in three lobes (as shown in Fig. 2) [12] called the Forescatter Lobe, the Normal Lobe, and the Backscatter Lobe.
Figure 3 is an example of a measured BRDF [ 141which displays all the lobes. This class of surfaces can be described by the reflectance maps:
& /
R(i, n, r) = aI &(i, r
i W,Y)
Object FIG. 1. The imaging geometry.
+
a2fn,,di,
+
a3
n, r) iTn n,
r) iTn
fbdi, n, r) iTn for iTn 2 0
(3)
= 0 otherwise, where the (Y’Sare the albedos and f&J ) is the BRDF due to the forescatter lobe, fnorm() due to the normal lobe, fbsc() due to the backscatter lobe.
SHAPE AND REFLECTANCE Surfacenormal
Reflecting surface FIG. 2. The structure of reflected light.
The forescatter lobe is spread around the specular direction and is a surface phenomenon [l, 2, 141. A rough surface is made up of facets, each facet being a small, pertlectly plane reflector. The facets are inclined randomly about the mean surface. When light is incident on them, each facet reflects like a plane mirror. The distribution of the facets about the mean causes the reflected flux to distribute around the specular direction giving rise to the forescatter lobe [ 141. The normal lobe is spread around the surface normal. (See Fig. 2.) A survey of the data in the literature [7, 1l] shows that the distribution of light in this lobe almost
Forescatter
Normal
always obeys Lambert’s Law, i.e.,&,& ) = l/r. Hence, an explanation of the normal lobe is essentially an explanation of Lambert’s Law. The normal lobe originates from scattering in the bulk of the material. Bulk can be modeled as an optically uniform material filled uniformly and randomly with inhomogeneities. Scattering from bulk is analyzed by radiative transfer theories [3]. The results of radiative transfer theory are reported in [31 and a discussion of them can be found in [12]. From the results of radiative transfer theory it is possible to qualitatively demonstrate the origin and two main properties of Lambert’s Law-the independence of BRDF from the incident and receiving directions and the independence of BRDF to the type of bulk material. The reader is referred to [ 1I, 121for explanations. Radiative transfer theory also predicts that scatterers close to the surface scatter incident radiation back along the incident direction. This creates the backscatter lobe. For a large class of surfaces, the amount of backscatter is negligible. If we use the exact expressions for the BRDF of each of the lobes, the resulting reflectance map obtained from (3) is quite complex. It is shown in [12] that for most photometric stereo geometries (phase angle of about 30”) and for most surfaces, the expression for R(i, n, r) can be simplified to R(i, n, r) = LY~exp{-c2[coss’(nTn)]‘} + a2(iTn) + a3,
-40"
Lobe
Lobe
0.0"
0, (degees)
THE PROBLEM
OF SIMULTANEOUS
ESTIMATION
80”
alO" x 0 -7 0
(4)
where n, = (i + r)/(li + rll is the specular direction and c is a constant that depends on the surface roughness and the wavelength of incident light. The first term in the above expression is due to the forescatter lobe, the second to the normal lobe, and the third to backscatter. Note that n, lies in the principal plane with a zenith angle of OnY= 8i12 and that for a very rough surface typical values for c are around 2.5 [14]. For the surfaces used in this paper, the value of c is experimentally determined to be 2.578 [ 131. For this value of c, the lobes are not clearly separated [12]. Though we use this specific value of c, the conclusions drawn are independent of it and hold for any value of c that causes the lobes not to be clearly separated. We call the reflectance map of (4) the “physical reflectance map.” 3.
-80"
277
FROM STEREO
9; 8; 8; 0;
= = = =
-30° -60 -60’ -750
FIG. 3. An example illustrating the three lobe structure.
Consider the situation in a classical photometric stereo apparatus. Assume that there are k identical light sources illuminating the object from directions il, . . . , ik. Each light source is turned on one at a time and a set of k images is obtained. At every pixel, we have k intensities given by
278
Z(h, n, r)
= 1
TAGARE
AND DeFIGUEIREDO
Z mm
i Z(ik,In, d
Substituting the physical reflectance map of Eq. (4) (without the backscatter term) into the above equation and assuming that the pixel is not in self shadow in any of the images, we get
Two assumptions are needed to make the analysis tractable: 1. The light sources are placed at a constant zenith angle Bi. Hence, Z(i, n, r) is a continuous function of the azimuth angle + and can be represented as Z($i), 2. cos 0 is approximated by its Taylor series expansion cos 8 = 1 - e2/2. Therefore, O*= 2(1 - cos 0), i.e., (toss’ x>* = 2(1 - X). 4. UNIQUE VALUES FOR SIMULTANEOUSESTIMATES
/exp{-c2
In this section, we investigate whether unique simultaneous estimates can be obtained from the data. By considering the continuous function Z(i, n, r) it is possible to establish that & can be obtained independently of the rest of the estimates. In fact, a closed form solution can be found for $,,. The other variables do not seem have a closed form solution, but it is possible to show that they do have a unique solution. Therefore, unique simultaneous estimates can be obtained. The analysis of this section also gives the number of light sources necessary to do so.
[cos-’ nf n12}\
=Zmaxffl I
exp{-c* [toss’ nf, nl*} I
+
Zmax~2
:
0
i;fn
/exp{-c*
=al 1
4.1.
[cos-’ n:, n12}\
4.1.1. Solution from the Continuous Function Z(i, n, r). Consider how the intensity Z( ) varies as a function of +i for fixed 8;) &, and 19,.From Eq. (4),
exp{-c2 [cos-’ nzknl*)1
/iTn\
Solving for &
Z(&) = Ul exp(-c*[cos-’ COS Bi/2 COS 19, + sin 8;/2 sin 8, COS($; - &)12) + U* COS Of COS 8, + a2 sin Bi sin 8, cos(+j - &),
(6)
where, al = Zmax~iand a2 = Zma+z. As indicated in Section 1, the problem of simultaneous estimation is posed as the following series of questions: 1. Is there a unique solution to Eq. (6) for aI, ~2, and n? How many light sources do we need to get this solution? 2. Can we use any property of the equation so that some of the unknowns may be obtained without solving the equations for all the variables? 3. Is the solution numerically robust? Equation (6) is nonlinear and it is not possible to analyze it exactly. We employ the following technique to analyze it by approximations. The left hand side of Eq. (6) is viewed as samples of the continuous function Z(i, n, r) obtained at ii, . . . , ik. For fixed n and r, it is possible to show that Z(i, n, r) is essentially a band limited function and that Z(il, n, r), Z(iz, n, r), . . . , I(&, n, r) are Nyquist samples of it. Then it is possible to establish that simultaneous estimates can be obtained from the continuous function Z(i, n, r). Under the assumption of band limitedness they can also be computed from the samples.
(7)
where the fact that 8, = 8i/2 is used. Z($i) is a periodic function of $i and therefore has a Fourier series expansion. When & = 0, Z($i) is an even function of +i and has a cosine series expansion, I(&)
=
PO*
+
C P
2
Pp
COS(P$i)
P,*
COS(Pd%)*
When r$,, f 0, I(&)
=
po
+
$-
P
C P
ap
WP~i),
where, PO= P& Pp = Pp*C~S(PhA ffp = P,* sin(Mh). Hence, +,, is obtained as the phase of the first harmonic, i.e. & = tan-’ al/Pi.
SHAPE AND REFLECTANCE
4.2.
Solution
to &from
279
FROM STEREO
Samples of Z(i, n, r)
! ! b
p+l
To construct the solution to $n from samples of Z(i, n, r) we first show that Z($i) is essentially band limited.
1,
and the spectrum of g(cos Cpi)falls off. In fact, as shown in Section 4.2.2, the fall-off is very rapid. We approxi4.2.1. The Spectral Behavior ofZ( +J. Using the ap- mate the situation by saying that g(cos +i) is a bandproximation (cos-’ x)~ = 2( 1 - x), I( 4i) can be written as limited function and that it can be sampled at the Nyquist rate without losing any information. Rewriting Eq. (8) in terms of g( ), we get I(+) = al exp(-c2[cos-’ cos 8ii2 cos 8,, + sin 8i/2 sin en COS(+i - $n)]‘) I( C#Ji)^- al exp(-2c2 + 2C2 COS Oil2 COS 0,) + a2 COS 6i COS On X g(cos($; - &)) + a2 cos 8j cos 19~ + a2 sin 8i sin en COS(+i - &) + a2 sin Bi sin en COS($i - +,). = al exp(-2c2 + 2C’ COS Oil2 COS On + 2c2 sin 19;/2sin 9, cos($i - 4,)) When 4,, = 0, I( ) is the sum of g( ), a D.C. term, and a (8) + a2 COS 8i COS en first harmonic term. Hence, the spectral nature of g( ) is + a2 sin Bi sin 0, COS(+i - &) inherited by Z(&) and we can sample Z(+i) at the same rate as we sampled g(cos $i). Sampling Z($i) consists of = L11exp(-2c2 + 2C2 COS di12 COS 0,) uniformly placing as many light sources as needed in one - &)) X exp(2c2 sin 8;/2 sin 0, COS($i period at a constant zenith angle and imaging the surface + a:! cos 0; cos 8, with one light source at a time. + aI sin 0; sin 0, COS($i - 6,). A nonzero & corresponds to a circular shift of Z($J. Therefore, using standard DFT theory, & can be estiThus, the variation of I( ) with $i comes from three mated by terms. The second and the third term have only D.C. and first harmonic content and are completely band limited. Im DFT(Z(@i), 1) The variation of Z( ) with $i in the first term in Eq. (8) & = tan-’ DFT(Z($i), 1)’ Real comes from the factor exp(2c2 sin 8i/2 sin 8, cos(~i 4,)). To investigate the spectral behavior of this factor, where, if k is the number of samples, consider the function g( ) g(x) = exp(cyx) (Y> 0.
Then g(cos 4i) = exp(cr cos $i) and the spectral behavior of the first term can be analyzed by investigating the spectral behavior of g(cos +i) with (Y= 2c2 sin 8i/2 sin 8,. As g(cos $i) is an even function of $i, its Fourier series expansion has only cosine terms. Let
g(COS
DFT(Z($i), 4.2.2.
-bP+l <-
a
p > [4(2C2
sin 8i/2 - I)],
I b 1 -2p+1 = 2c2sin 8i/2 sin en 2ptl * If
where [xl denotes the smallest integer greater than X, then
of Light
Sources for Obtaining
4,,. The number of light sources needed to obtain the solution to & is given by the Nyquist sampling rate for g(cos $i). Denoting the harmonic at which the spectrum of g(cos +i) begins to drop off by p*, we get from the previous analysis
$i) = bo + 2 b, COS(p$i). p=l
In Appendix 1, it is shown that
The Number
m) = $ Z(+il)e-j27im”k.
P * = [&(2C2
sin 8i/2 - l)].
For c = 2.578 and Oi = 40.0”, we get p* = 2. Hence, the spectrum can be expected to start falling after the third harmonic. However, the spectrum falls off much faster than this bound. Figure 4 shows the amplitude spectrum of g(cos 4i) with 8i = 40” and 8, = 20” for different values of c. Note the sharp cutoff. In particular note that for c = 2.578 the function g(cos $i) is essentially band limited to the fourth harmonic and eight or more light sources are sufficient to obtain & from Eq. (9). 4.3.
Solving for en, al, and a2
Having obtained &, consider whether 8,, al, and a2 can be obtained from the continuous function Z(i, n, r).
280
TAGARE AND DeFIGUEIREDO
S,*,@, E*) that give the same Z(i, n, r) for all Xi E [O, I]. This implies that the following equation holds for all Xi E
[O,11: aTexp(c* Xi) - dTexp(c’* Xi) = (6; - at)Xi + h3*- af. (11)
4
2
"0
P
6
8
10
FIG. 4. Amplitude spectrum of g(cos 4,) for different c.
4.3.1. Solution from the Continuous Function Z(i, n, r). Rotate the reference frame in the azimuth plane such that in the new frame c#+,= 0 and all the zenith angles are the same as before. In this frame, Eq. (8) becomes
Z(i, n, r) -L al exp(-2c2 + 2c2 Cos 8i/2 Cos 0,) x exp(2c2 sin 1!&/2cos en cos +i) + a2 sin 8; sin 8, cos +i
The functions on the left and right hand sides of the above equation are linearly independent for Xi E [O, 11. Hence the above equation cannot hold for all Xi E [O, 11. Therefore, two distinct sets (a:, at, a:, c*) and ($, $, 6:, e*) cannot be found to give the same Z(i, n, r) and a?, a?, and a: can be uniquely determined from the continuous function Z(i, n, r). 4.3.2. Solution from Samples of Z(i, n, r). Suppose that light sources are placed so that there are m distinct values for Xi denoted as Xi,, Xi23. . . , Xi,. The intensities in the image due to these lights are I,, Z2, . . . , I,,, given by
+ a2 cos 8i cos 8,
= UTexp(c* Xi) + &Xi + at,
(10)
\Ll
where a; = aI exp(-2c2 + c2 cos 8J2 cos e,), a2*= a2 sin 8i sin e,, af = a2 cos Oi cos e,,
\a: exp(c* Xim)f
U:X, + U3 *I
Suppose that the system of equations (12) cannot be uniquely solved for a?, al, and a:. Then, as before, it is possible to find two distinct sets (a:, a:, a$, c*) and (L$, ti,*, 6:, 6) that give the same left hand side for Eq. (12),
c* = 2c2 sin oil2 cos en, and Xi = COS$+a The variable Xi (and hence Z( )) is an even function of $i in the range (-7~, +r 1. Therefore, we need consider only the range [0, + x] and in this range Xi takes values from 0 to 1. Also, 8,) al, and a2 are related to a?, a:, and a:’ by
LITexp(c*
Xi,)
- hT exp(E* Xi,)
UTexp(c*
Xi*)
- 67 exp(C*
Xi2)
a? exp(c* Xi,) - ciTexp(E* Xi,)
(13) (ti$ -
az*)Xi,
+ ifi: - a$
(ciz*- a:)Xi, + tit - a: =
a1 = exp(-2c2 + c2 cos Bi/2 cos e,y
d
(tit -
a,*)Xim
+ ti: - a3*I.
Equation (13) is Eq. (11) evaluated at Xi,, Xi*>. . . , Xi,. As shown in Appendix 2, Eq. (11) can hold for at most three different values OfXi. Thus, if we set m 2 4, Eq. (13) Thus, if a:, at, and a$ can be uniquely determined from will not hold for some row, thus assuring that Eq. (12) will have a unique solution for a?, at, and a: and thereZ(i, n, r), then so can 8,) al, and a2. Suppose that a?, at, and a: cannot be uniquely deter- fore also for f3,, al, and a2. Due to the symmetry in light mined from Z(i, n, r) as +i goes from 0 to W. Then it is source placement, we need k = 2m light sources to give possible to find two distinct sets (a:, at, a;, c*) and (HT, us m distinct values of Xi. This is shown in Fig. 5. a2 =
sin 8, sin 8i’
SHAPE AND REFLECTANCE
281
FROM STEREO
5.1.
Linear Analysis
The function exp{-c2 [cos-*( )12}can be approximated quite well as a linear function over a large range of its domain. Figure 6 shows the function for c = 2.578 and two linear tits over different ranges. Let us assume that for a particular surface normal n and a set of light sources, the sequence of numbers nfn lies in one of these regions. Then the function 4(x) = exp[-c2{cos-‘}2] can be approximated as 4(x) = +(x0) + (X - x0) C$‘(x0). With this approximation, Eq. (6) can be written as
/ Y FIG. 5.
1
I&, n, r) Eight lights for four different
x,.
i
:
\
= al($(xo) - xo~'(al))
Thus, &, being known, a unique solution can be obtained for f&, al, and a2 with k = 8 (m = 4) light sources. 4.4.
:
n
The Total Number of Light Sources
From the above, it is clear that the number of light sources sufficient to guarantee a unique simultaneous estimate of all unknowns is the maximum of eight and the number obtained from Nyquist sampling of Z( 6;). For c = 2.578 the bandwidth of Z(&) is sufficiently small that eight light sources are sufficient for Nyquist sampling. In general, Fig. 4 shows that the number of light sources needed for Nyquist sampling of I( 4;) is not too different from eight and therefore we expect that the number of light sources needed for a unique answer would not differ very much from eight. All the experiments in this paper are conducted using eight light sources. 5.
i
ILL-POSED
NATURE
OF THE JOINT ESTIMATE
The previous section demonstrated that there is a unique answer to the simultaneous estimates and that this estimate can be obtained with a finite number of light sources. However, the analysis of existence and uniqueness of the estimate does not inform us how robust the estimate is with respect to additive noise in the data. In this section, we explore the robustness to noise of the solution to Eq. (6). As we shall see, the solution is sensitive to noise and therefore the simultaneous estimation problem is ill-posed. Since Eq. (6) is nonlinear it is difficult to demonstrate this rigorously. We demonstrate the ill-posed nature in two ways. First, we linearize Eq. (6) and show that the linear equation is ill-posed. Besides establishing grounds for belief in the ill-posedness, the linear analysis also gives us qualitative insight into its origin. Second, we numerically compute the right hand side of Eq. (6) and show that large changes in n cause small changes in the intensity.
where
and
/sin0;,/2 cos &, P=(
sin I$,/2 sin &
cos &,I2
sin ej,12 sin 4,
cos 8,,12
: \
sin 8i,12 COS +ii
sin Oj,Cos $j,
sin 8i, sin C& cos Bi, 1 1
sin Oii Cos &
. .
sin 0,, sin &A cos 0,, 1 9
Note that the vector u has vertical dimension of 7 (each n has a dimension of 3). Under the assumptions that I$,, B;,, . . . ) 81,are in the range [O, 30”], we have the following: 1. cos 8,,/2 and cos 8, are close to 1, thus columns 3, 6 and 7 in the matrix P are close to being linearly dependent, and
282
0.0
TAGARE AND DeFIGUEIREDO
1:
In Fig. 7, note that although the maximum difference between I,( > and Z,( ) is about 40 grey levels, the mean difference is only 6.1 grey levels. Also, the maximum difference occurs at Bi = 0.0”. For nonzero 0i, the difference falls off quite rapidly. The mean difference of 6.1 grey levels is well within the range of noise generated by most cameras. Hence, we conclude that it is difficult to distinguish between the two surface normals and the two “albedos” even by using a large number of light sources in the region 0 I +i < 360” and 0 5 8i i 30”. The problem is clearly ill-posed.
The function \
0.0
x
If0
FIG. 6. Linear approximation to the forescatter lobe function.
2. sin 04/2 = i sin 8,, hence the pairs of columns 1,4 and 2,5 are close to being linearly dependent. Thus the matrix P has three linearly independent columns and hence only three quantities in the vector u can be estimated with any accuracy at all. Since our ultimate goal is to estimate four quantities (a,, a2, &, &J, the problem is ill-posed. Note that the ill-posed nature of the problem is independent of the vertical dimension of P and hence independent of the number of light sources. Qualitatively, the above argument shows why simultaneous estimation is ill-posed: 1. The part of the observed intensity vector that lies along [l 1 . . . llT cannot be easily decomposed into contributions due to the forescatter and normal lobe, i.e., there is no way of telling how much of the “constant” part of the intensity values was generated by either lobe, and, 2. For the range of incidence angles we are interested in, the variation in intensity due to the change in the zenith angle 8ikcannot be uniquely decomposed into contributions due to either lobe. 5.2.
Numerical
Analysis
For the exact nonlinear equation (7) the ill-posed nature is shown in Fig. 7. The figure shows a plot of ]]Zi(&, 4i) - Z2(0i, +i)JI for 0 5 Bi 5 30” and 0 or +i < 360”. The functions Zi( ) and Z,( ) are specific cases of Eq. (7), where the following parameters were used. For Ii ( ), ai = 140.0, a2 = 70.0, c = 2.578, 8, 33.0”, $,, = O.o”, and for
6. REGULARIZATION
OF SIMULTANEOUS ESTIMATES
Ill-posed problems can be solved reliably by regularization. We regularize Eq. (6) by estimating n, aI, and a2 by minimizing
Q= hl(Z- z)T(Z- Z)+ A& - a:>2+ h&2 - ap, (14)
where Z = vector of observed intensities, / exp{-c2[coss’n$]2}
\
I = zmaxcq I
exp{-c?[cos-1n$]2}
/iT n\ +
Zmax~2
I
:
3
i*ii;fn (16)
hi, Xz, A3> 0, and &and c$are prior estimates of (Y~and CX~. This is a nonlinear optimization problem in n, (Y~and
40.0 \
Z2( ), al = 98.0, a2 = 0.0, c = 2.578, 0, = 10.32”, $,, =
0.0. Thus, Zi( ) corresponds to a surface normal with an azimuth angle of zero degrees and a zenith angle of 33.0”, while Z,( ) is for a surface normal with zenith angle of 10.32” and the same azimuth angle. Note that the values of ui, a2 are also different for both the cases (the values used in this example correspond to the actual grey level values recorded by the camera in our experimental setup).
(15)
FIG. 7. Ill-posed nature of the joint estimate
SHAPE AND REFLECTANCE
FROM STEREO
283
FIG. 8. A sphere.
(Yeand a closed form solution cannot be obtained for it. A standard nonlinear, numerical optimization procedure [6] is used to do this. The entire procedure for joint estimation is: 1. The surface is imaged with eight line sources placed at azimuth angles of O”, +45”, -t90”, ? 135”, and 180” at a constant zenith angle of 25”. 2. We know from Section 4 that these light sources are sufficient to obtain an estimate of &. The estimate is obtained via Eq. (9). 3. The other parameters are obtained by minimizing Q in Eq. (14) with the constraint that the value of the azimuth angle & be that given by the above step. 7.
EXPERIMENTAL
and a; = 0.5. The value of I,,, = 140.0 was used in Eq. (15) while (Y: and a:, the prior guesses, were set to 1.1 and 0.45 respectively. The constants in the minimization function of (15) were Xi = 100, h2 = 1600, and A3= 900. As we shall see, this is not an excessive amount of regularization. The initial values for the optimization procedure were n = [0 0 llT, (Y~= 1.2, and a2 = 0.8. Figure 11 shows the reconstruction of the spherical surface.
A spherical surface was fit to the reconstruction and the r.m.s. error between the best fit surface and the reconstructed surface normals was computed. The error was measured as the mean square difference in the angle
RESULTS
Figures 8, 9, and 10 show the three objects that were used in experiments to validate the theory. It was easily determined that for the illumination conditions of this experiment, the reflectance map for the surface would be best described by I,,,,, = 140.0, CY~ = 1.O,
FIG. 9. A cylinder.
FIG. 11. Regularized reconstruction of the sphere
284
TAGARE AND DeFIGUElREDO
FIG. 12. Histogram for (Y,.
between the best fit surface normals and the reconstructed surface normals. This error had a value of 5.21”. Figures 12 and 13 show the histograms of the value of (~1and (~2as returned by the optimization procedure for the entire surface. The peaks in the distribution occur for aciat 1.04 and (~2at 0.34. From the distribution, note that spread of the values is fairly large and corresponds to a mild amount of regularization. It was also observed that for this range of values of o1 and a2 the contribution from the second and third term in Eq. (14) was quite small compared to the contribution from the first term. This further supports the argument that the amount of regularization was not excessive. Finally, Figs. 14 and 1.5show reconstructions for the cylinder and spoon of Figs. 9 and 10. The exact same procedure was used in these reconstructions as was used for the spherical surface. 8. CONCLUSION
The problem of simultaneous estimation of reflectance map and shape from photometric stereo was posed as a series of three questions about the existence, uniqueness, and robustness of the estimates. We investigated simultaneous estimation for diffuse surfaces whose lobes of the reflectance map were not easily separable. In the first
w
1.0
FIG. 13. Histogram for (Y?.
FIG. 14. Regularized reconstruction of the cylinder.
part of the paper, we showed that a unique answer was available to the simultaneous estimates. Further, it was possible to show that the azimuth angle of the surface normal could be estimated independently of the rest of the unknowns. One of the main contributions of this paper is to establish the number of light sources sufficient to obtain the unique estimate. By an approximate theoretical analysis, it was shown that a small number of light sources yield estimates, but that the estimates are sensitive to noise in the data. From a linear analysis we argue that the sensitivity to noise could not be reduced by adding extra light sources and that prior estimates of some of the parameters of the map are required to regularize the answer.
\
2:o
FIG. 15. Regularized reconstruction of the spoon.
SHAPE AND REFLECTANCE
Experiment with the regularized simultaneous estimation procedure yields good results. Although most of the arguments in this paper are developed in the context of a specific reflectance map, the techniques of Appendix 1 and 2 are fairly general and can be easily used to make similar arguments about other reflectance maps.
285
FROM STEREO
Repeatedly integrating by parts, and using the fact that for 1 I m 5 p - 1, (d”ldy”){(l - Y~)“-“~} = 0 at y = * 1, we get
1’_,g(y)$, ((1- Y*Y-I’*)
dv
= (-l)P j;, g’p’(y)( 1 - Y*)P-“~ dy, APPENDIX
1: THE FOURIER
SERIES FOR g(COS+)
where g@‘(y) denotes the pth derivative of g( ). Thus,
Let g(COS 4i) = bo + C bp COS p$i.
(-1)2P 2pr(p + $>i ‘, (1 -
bp = G
P
Then, for p 2 1,
and
bP= i I,‘”gtcos +i)COS p4i
d4i
b
rep
-= P+’
=-
i
I,” dcos
$i> COS p+i
d$i*
bP=5 j:,
g(y)
VT%
dy,
1
_
2~ + 1
bP
T (y> = (-l)P e P 2p+’ r(p
V’& dP dyp {(I - y*)p-“*}. + 3)
+)
JLl
(I
-
1s”
dP dy”
_
2pcl+l
lh
< -
=-
Evaluating the integral by parts,
1’_,g(y)$ ((1- Y’)~-“~) -j&
a
%+I
=-I
I’ _] g’(y)
(17)
(1
y*)(l
-
(1
I.f-‘I
(1
y*>(l
-
(1 -
-
y2)p-“*
(1
-
-
y21p-“*
Y*)~-“~
y21p-“*
y2Y-“*
s@+“(r) P(r)
P(y>
,@“(u>
P’(y)
41 4
&I 41
&I
/j-i1 (1 - Y~)~-“~ gck)(y) dyl
a
2p + I
= 2c2 sin Oil2 sin 0, 2p+l ’ dv
((1 - y2)“-“2})I, APPENDIX
-
dy’
J-1, (1 - y2)p+‘” g@+“(y) dy j-L1 (1 - y2)p-“* g’p’(y) dy ’
1s”
x ((1 - y*)~-“~) dy
x ((1 - Y~)P-“~} dy.
dy
J’, (1 - y*)~-“~ g’“‘(y)
lj-$
(-1)~ m 6 2”+’ I’(p + +)
g’p+‘)(y)
Since g(y) = exp(cuy), we have g(p+“(y) = ag’p’(y) and since in our case (Y > 0, we have g@‘(y) > 0 and g@+‘)(y) > 0. Further, in the range of integration the term (1 - y*) is always between 0 and 1. Therefore, from Eq. (17) we get
Substituting this into the above equation, g(y)
Y*)P+“~
From this, we get
bP+’ -_-~
where TJy) = cos(p cos-’ y) is the Tchebycheff polynomial of degree p. T,( >are a set of orthonormal polynomials, and are generated by Rodrigues’ formula.
+
2Up + 8)
4.J
Changing the variable of integration to y = cos $i, we get
= g(y)
g@) (Y) dy
Y~)~-“~
s
‘, ‘F’(Y) g’
{(I - Y~)“-“~} dy
2: THE NUMBER OF SOLUTIONS EQUATION (11)
TO
Note that the equation ((1 - y’)p-I’*}
dy. (~1
exp(clx)
-
(~2
exp(c0)
= plx + p2
(18)
286
TAGARE AND DeFlGUElREDO
can have at most three solutions for x given any values of al, QI~,pi, p2 and cl, c2. This is seen as follows: Denote the left hand side of Eq. (18) byf(x) and rewrite it as f(x)
-
PIX
+
P2 =
(19)
0.
Ifx,, x2, . . . ) x, are the solutions of Eq. (18) placed in ascending order, then, between any two consecutive solutions, the derivative of the left hand side of Eq. (19) goes to zero at least once, i.e., if the number of solutions to f’(x)
-
Pl
=
(20)
0
is m, then n I m + 1. Using the same argument on Eq. (20) we have that if f”(x)
= 0
(21)
has k solutions, then m I k + 1. Hence, n I k + 2. But Eq. (21) can be written as cz,c: exp(c,X) = u2c: exp(czx). Taking logarithms, clx = c2x + log (Y~c:- log (Y,c: which is a linear equation and has a single solution. Hence, k = 1 and n = 2. Thus Eq. (18) has at most three roots. (It is easy to find (~1, 02, cl, c2, j31, ~32that cause Eq. (18) to actually have three roots.) Referring back to Eq. (1 l), suppose that for four or more different values of xi, there are two different sets (a:, a:, a$, c”) and (@, @, b:, C*) that give identical values of Z(xJ for some values Of Xi. Then, from Eq. (lo), for all i = 1, 2, 3, 4, . . . , UTexp(c$xJ +
U2*xi
+ U: = LjTeXp(EfXi)
+ Lif,Xi+ hf,
i.e., UT eXp(CfXi)
-
Lif eXp(E*Xi)
=
($
-
Uz*) Xi +
i=
b3* -
Ut,
1,. . . ,4.
This has the same form as Eq. (18) for which there cannot be more that three solutions. Hence, Eq. (11) cannot have more than three solutions. ACKNOWLEDGMENTS We gratefully acknowledge the help we received in regard to Appendix 1 from Dr. Don H. Johnson of the Department of Electrical and Computer Engineering of Rice University.
REFERENCES 1. P. Beckmann and A. Spizzichino, The Scattering of Electromagnetic Wavesfrom Rough Surfaces. Pergamon, Elmsford, NY, 1963. 2. R. C. Birkebak and E. R. G. Eckert, Effect of roughness of metal surfaces on angular distribution of monochromatic reflected light, J. Heat Transfer, February 1965, 85-94. 3. S. Chandrasekhar, Radiative Transfer, Dover, New York, 1960. 4. E. N. Coleman and R. Jain, Obtaining 3-dimensional shape of textured and specular surfaces using four-source photometry, Comput. Graphics Image Process.
18, 1982,309-328.
5. B. K. P. Horn, Robot Vision. M.I.T. Press, Cambridge, MA, 1986. 6. IMSL, ZMSL Math Library, Houston, TX. I. G. Kortum, Reflectance Spectroscopy, Springer-Verlag, Berlin/ New York, 1969. 8. S. K. Nayar, K. Ikeuchi, and T. Kanade, Extracting Shape and Reflectance of Lambertian, Specular and Hybrid Surfaces, Technical Report CMU-RI-TR-88-14, C.M.U., August 1988. 9. W. M. Silver, Determining Shape and Reflectance Using Multiple Images, master’s thesis, MIT., 1980. 10. J. M. Stone, Radiation and Optics, McGraw-Hill, New York, 1963. 11. H. D. Tagare, A Theory of Photometric Stereo for a General Class of Reflectance Maps, Ph.D. thesis, Rice University, Houston, Texas, 1989. 12. H. D. Tagare and R. J. P. deFigueiredo, A framework for the construction of general reflectance maps for machine vision, C.V.G.Z.P., to appear. 13. H. D. Tagare and R. J. P. deFigueiredo. A theory of photometric stereo for diffuse, non-Lambertian surfaces, Z.E.E.E. Trans. Puttern An. Mach. Intelligence, 1991. 14. K. E. Torrance and E. M. Sparrow, Theory of off-specular retlection from roughened surfaces, J. Opt. Sot. Am. 57, 1967, 11051114.
IMAGE
UNDERSTANDING
Vol. 55, No. 3, May, pp. 275-286,
1992
Simultaneous Estimation of Shape and Reflectance Map from Photometric Stereo D. TAGARE*
HEMANT Department
of Electrical
and Computer Engineering,
Rice University,
Houston,
Texas 77005
AND
Rur J. P. DEFIGUEIREDO Department
of Electrical
Engineering
and Department Received
of Mathematics,
of California
isolate
California
92717
the Lambertian
component.
A similar
strategy
is
also used in [8]. Silver [9] considered the problem of simultaneously estimating reflectance map parameters and shape without assuming that the components of the reflectance map can be easily separated. Using his technique, the Lambertian part of the reflected light can be removed, thereby isolating the specular
part.
This technique
has not been tested
experimentally, and Silver himself cautioned about the numerical robustness of the algorithm. In this paper, we consider simultaneous estimation of reflectance map parameters and shape in the context of diffusely reflecting surfaces. Such surfaces occur commonly, and reflected light from them has several components. The geometric distribution of light in these components cannot always be easily decomposed, and the techniques developed in [4, 81 cannot be applied. Therefore, it is necessary to develop alternate simultaneous estimation techniques for such surfaces. We begin this development by proposing a class of reflectance maps for such surfaces and then investigate how the parameters of the map and the object shape can be estimated from photometric stereo. Specifically, we consider the following questions of existence, uniqueness, and robustness: 1. The number of light sources: Before launching into specific techniques for obtaining the solution, it is essen-
INTRODUCTION
Shape from shading and photometric stereo techniques estimate the shape of an object from the intensities in an image using knowledge of the light source positions, the viewing direction, and how the object reflects light. The information about light source positions, viewing direction, and reflectance is mathematically described by the reflectance map. Accurate prior knowledge of the reflec-
tance map is required for precise shape reconstruction. Such knowledge is often not available even for controlled conditions, and under these circumstances it is desirable to estimate the reflectance map along with the shape. * Current address: Department of Diagnostic ment of Electrical Engineering, Yale University, icut 06510.
Irvine,
Algorithms for estimating parameters of reflectance maps which are sums of Lambertian and specular components have been proposed before. Most of these are based on the observation that the geometric distributions of reflected light in Lambertian and specular components are so different that the two can be easily separated. Coleman and Jain [4] report a procedure to estimate the shape of an object and the albedo of the Lambertian component of the reflectance map. Four point light sources are used, and the intensities in the images are analyzed to
Press, Inc.
1.
at Irvine,
May 1, 1990; accepted July 11, 1991
This paper investigates whether the shape of an object and certain parameters of its reflectance map can be simultaneously estimated using photometric stereo. This problem has been addressed in the literature for the case where the Lambertian and non-Lambertian components in the image can be easily separated. No such separability is assumed in this paper. A class of reflectance maps for modeling diffusely reflecting surfaces is proposed. This class is based on the physics of scattering from real world surfaces. Next, the problem of joint estimation of some parameters of the map along with the surface shape is analyzed. A bound is obtained on the number of light sources necessary for a unique solution to the problem. The analysis also reveals that some of the estimates can be obtained by a nonparametric method. The behavior of the estimates in the presence of noise is also investigated. It is shown that simultaneous estimation is ill-posed. Regularizing the estimates yields good reconstructions from real world data. 0 1992 Academic
University
Radiology and DepartNew Haven, Connect-
275
All
1049.9hfXVY2 $5.00 Copyright G 1992 by Academic PI-~\\. Iw. rights of reproduction in any form rewrved.
276
TAGARE AND DeFIGUEIREDO
tial to check that the solution exists and is unique. Therefore, we ask, is there a unique solution to simultaneous estimation and how many light sources are needed to guarantee it? 2. Separation of variables: Can the problem be separated into simpler problems? Even though the original problem is complex and might not exhibit a closed form solution for all the unknowns, perhaps some of the unknowns might be available in a simple manner, while the rest may be obtained numerically. 3. Robustness with respect to noise: How robust is the solution with respect to noise? In Section 2, we review the reflectance map model. In Section 3, we formulate the problem of simultaneous estimation of the reflectance map parameters and surface shape. Section 4 addresses the problem of obtaining a unique solution. Section 5 illustrates the ill-posed nature of simultaneous estimation and Section 6 regularizes it. Section 7 contains experimental results. 2. THE REFLECTANCE MAP
Figure 1 illustrates the imaging geometry. The object is illuminated by a point light source at infinity. Bold face n, i, r denote unit vectors along the surface normal, the incident direction, and the receiving direction, respectively. The reference frame is aligned so that r is along the positive z-axis. The plane spanned by i and r is called the “principal plane” and the angle between i and r the “phase angle.” The zenith and azimuth angles of any unit vector are denoted by 0 and 4. Subscripts denote the vector associated with the angles; e.g., the zenith and azimuth angles of n are denoted by 8, and $,, respectively. The zenith angle is measured positively down from the z axis, and the azimuth angle is measured anticlockwise from the x axis.
Light source m
The intensity in the image, Z(x, y), is given by [5] Z(x, Y) = L,R(i,
0,
Y>, r),
where R( > is the reflectance map. For a point light source, R( ) is given by [5] R(i, n, r) = afi(i, n, r) iTn, for iTn, 2 0 = 0 otherwise.
(2)
where fr( ) is the BRDF (the Bidirectional Reflectance Distribution Function) of the surface and (Yis a positive constant called the “albedo” [5]. The BRDF is defined as the ratio of the radiance of light reflected along r to the irradiance of incident light on the surface.
2.1.
The BRDF of Rough Surfaces
In general, the BRDF of any surface is a sum of delta functions and functions of bounded variation. The delta functions are attributed to Specular Rej7ection and the rest to Diffuse Rejection. As a first step toward obtaining a theoretical expression for the BRDF of a solid, the solid can be modeled as a regular array of particles having a planar boundary at the top. Reflection from the boundary is analyzed by considering how individual particles scatter and how their behavior is modified by their being together in a dense array. The Ewald-Oseen theory describes the reflectance behavior of such a configuration and the reader is referred to [lo] for a qualitative discussion of it. The main result of this theory is that perfect solids can only reflect specularly. In this paper, we ignore specular reflection. Real world solids are not perfect. At the scale of the wavelength of visible light their surfaces are rough and their bulk is inhomogeneous. These imperfections cause light to be diffusely reflected in nonspecular directions. Diffuse reflection from most engineering materials is contained in three lobes (as shown in Fig. 2) [12] called the Forescatter Lobe, the Normal Lobe, and the Backscatter Lobe.
Figure 3 is an example of a measured BRDF [ 141which displays all the lobes. This class of surfaces can be described by the reflectance maps:
& /
R(i, n, r) = aI &(i, r
i W,Y)
Object FIG. 1. The imaging geometry.
+
a2fn,,di,
+
a3
n, r) iTn n,
r) iTn
fbdi, n, r) iTn for iTn 2 0
(3)
= 0 otherwise, where the (Y’Sare the albedos and f&J ) is the BRDF due to the forescatter lobe, fnorm() due to the normal lobe, fbsc() due to the backscatter lobe.
SHAPE AND REFLECTANCE Surfacenormal
Reflecting surface FIG. 2. The structure of reflected light.
The forescatter lobe is spread around the specular direction and is a surface phenomenon [l, 2, 141. A rough surface is made up of facets, each facet being a small, pertlectly plane reflector. The facets are inclined randomly about the mean surface. When light is incident on them, each facet reflects like a plane mirror. The distribution of the facets about the mean causes the reflected flux to distribute around the specular direction giving rise to the forescatter lobe [ 141. The normal lobe is spread around the surface normal. (See Fig. 2.) A survey of the data in the literature [7, 1l] shows that the distribution of light in this lobe almost
Forescatter
Normal
always obeys Lambert’s Law, i.e.,&,& ) = l/r. Hence, an explanation of the normal lobe is essentially an explanation of Lambert’s Law. The normal lobe originates from scattering in the bulk of the material. Bulk can be modeled as an optically uniform material filled uniformly and randomly with inhomogeneities. Scattering from bulk is analyzed by radiative transfer theories [3]. The results of radiative transfer theory are reported in [31 and a discussion of them can be found in [12]. From the results of radiative transfer theory it is possible to qualitatively demonstrate the origin and two main properties of Lambert’s Law-the independence of BRDF from the incident and receiving directions and the independence of BRDF to the type of bulk material. The reader is referred to [ 1I, 121for explanations. Radiative transfer theory also predicts that scatterers close to the surface scatter incident radiation back along the incident direction. This creates the backscatter lobe. For a large class of surfaces, the amount of backscatter is negligible. If we use the exact expressions for the BRDF of each of the lobes, the resulting reflectance map obtained from (3) is quite complex. It is shown in [12] that for most photometric stereo geometries (phase angle of about 30”) and for most surfaces, the expression for R(i, n, r) can be simplified to R(i, n, r) = LY~exp{-c2[coss’(nTn)]‘} + a2(iTn) + a3,
-40"
Lobe
Lobe
0.0"
0, (degees)
THE PROBLEM
OF SIMULTANEOUS
ESTIMATION
80”
alO" x 0 -7 0
(4)
where n, = (i + r)/(li + rll is the specular direction and c is a constant that depends on the surface roughness and the wavelength of incident light. The first term in the above expression is due to the forescatter lobe, the second to the normal lobe, and the third to backscatter. Note that n, lies in the principal plane with a zenith angle of OnY= 8i12 and that for a very rough surface typical values for c are around 2.5 [14]. For the surfaces used in this paper, the value of c is experimentally determined to be 2.578 [ 131. For this value of c, the lobes are not clearly separated [12]. Though we use this specific value of c, the conclusions drawn are independent of it and hold for any value of c that causes the lobes not to be clearly separated. We call the reflectance map of (4) the “physical reflectance map.” 3.
-80"
277
FROM STEREO
9; 8; 8; 0;
= = = =
-30° -60 -60’ -750
FIG. 3. An example illustrating the three lobe structure.
Consider the situation in a classical photometric stereo apparatus. Assume that there are k identical light sources illuminating the object from directions il, . . . , ik. Each light source is turned on one at a time and a set of k images is obtained. At every pixel, we have k intensities given by
278
Z(h, n, r)
= 1
TAGARE
AND DeFIGUEIREDO
Z mm
i Z(ik,In, d
Substituting the physical reflectance map of Eq. (4) (without the backscatter term) into the above equation and assuming that the pixel is not in self shadow in any of the images, we get
Two assumptions are needed to make the analysis tractable: 1. The light sources are placed at a constant zenith angle Bi. Hence, Z(i, n, r) is a continuous function of the azimuth angle + and can be represented as Z($i), 2. cos 0 is approximated by its Taylor series expansion cos 8 = 1 - e2/2. Therefore, O*= 2(1 - cos 0), i.e., (toss’ x>* = 2(1 - X). 4. UNIQUE VALUES FOR SIMULTANEOUSESTIMATES
/exp{-c2
In this section, we investigate whether unique simultaneous estimates can be obtained from the data. By considering the continuous function Z(i, n, r) it is possible to establish that & can be obtained independently of the rest of the estimates. In fact, a closed form solution can be found for $,,. The other variables do not seem have a closed form solution, but it is possible to show that they do have a unique solution. Therefore, unique simultaneous estimates can be obtained. The analysis of this section also gives the number of light sources necessary to do so.
[cos-’ nf n12}\
=Zmaxffl I
exp{-c* [toss’ nf, nl*} I
+
Zmax~2
:
0
i;fn
/exp{-c*
=al 1
4.1.
[cos-’ n:, n12}\
4.1.1. Solution from the Continuous Function Z(i, n, r). Consider how the intensity Z( ) varies as a function of +i for fixed 8;) &, and 19,.From Eq. (4),
exp{-c2 [cos-’ nzknl*)1
/iTn\
Solving for &
Z(&) = Ul exp(-c*[cos-’ COS Bi/2 COS 19, + sin 8;/2 sin 8, COS($; - &)12) + U* COS Of COS 8, + a2 sin Bi sin 8, cos(+j - &),
(6)
where, al = Zmax~iand a2 = Zma+z. As indicated in Section 1, the problem of simultaneous estimation is posed as the following series of questions: 1. Is there a unique solution to Eq. (6) for aI, ~2, and n? How many light sources do we need to get this solution? 2. Can we use any property of the equation so that some of the unknowns may be obtained without solving the equations for all the variables? 3. Is the solution numerically robust? Equation (6) is nonlinear and it is not possible to analyze it exactly. We employ the following technique to analyze it by approximations. The left hand side of Eq. (6) is viewed as samples of the continuous function Z(i, n, r) obtained at ii, . . . , ik. For fixed n and r, it is possible to show that Z(i, n, r) is essentially a band limited function and that Z(il, n, r), Z(iz, n, r), . . . , I(&, n, r) are Nyquist samples of it. Then it is possible to establish that simultaneous estimates can be obtained from the continuous function Z(i, n, r). Under the assumption of band limitedness they can also be computed from the samples.
(7)
where the fact that 8, = 8i/2 is used. Z($i) is a periodic function of $i and therefore has a Fourier series expansion. When & = 0, Z($i) is an even function of +i and has a cosine series expansion, I(&)
=
PO*
+
C P
2
Pp
COS(P$i)
P,*
COS(Pd%)*
When r$,, f 0, I(&)
=
po
+
$-
P
C P
ap
WP~i),
where, PO= P& Pp = Pp*C~S(PhA ffp = P,* sin(Mh). Hence, +,, is obtained as the phase of the first harmonic, i.e. & = tan-’ al/Pi.
SHAPE AND REFLECTANCE
4.2.
Solution
to &from
279
FROM STEREO
Samples of Z(i, n, r)
! ! b
p+l
To construct the solution to $n from samples of Z(i, n, r) we first show that Z($i) is essentially band limited.
1,
and the spectrum of g(cos Cpi)falls off. In fact, as shown in Section 4.2.2, the fall-off is very rapid. We approxi4.2.1. The Spectral Behavior ofZ( +J. Using the ap- mate the situation by saying that g(cos +i) is a bandproximation (cos-’ x)~ = 2( 1 - x), I( 4i) can be written as limited function and that it can be sampled at the Nyquist rate without losing any information. Rewriting Eq. (8) in terms of g( ), we get I(+) = al exp(-c2[cos-’ cos 8ii2 cos 8,, + sin 8i/2 sin en COS(+i - $n)]‘) I( C#Ji)^- al exp(-2c2 + 2C2 COS Oil2 COS 0,) + a2 COS 6i COS On X g(cos($; - &)) + a2 cos 8j cos 19~ + a2 sin 8i sin en COS(+i - &) + a2 sin Bi sin en COS($i - +,). = al exp(-2c2 + 2C’ COS Oil2 COS On + 2c2 sin 19;/2sin 9, cos($i - 4,)) When 4,, = 0, I( ) is the sum of g( ), a D.C. term, and a (8) + a2 COS 8i COS en first harmonic term. Hence, the spectral nature of g( ) is + a2 sin Bi sin 0, COS(+i - &) inherited by Z(&) and we can sample Z(+i) at the same rate as we sampled g(cos $i). Sampling Z($i) consists of = L11exp(-2c2 + 2C2 COS di12 COS 0,) uniformly placing as many light sources as needed in one - &)) X exp(2c2 sin 8;/2 sin 0, COS($i period at a constant zenith angle and imaging the surface + a:! cos 0; cos 8, with one light source at a time. + aI sin 0; sin 0, COS($i - 6,). A nonzero & corresponds to a circular shift of Z($J. Therefore, using standard DFT theory, & can be estiThus, the variation of I( ) with $i comes from three mated by terms. The second and the third term have only D.C. and first harmonic content and are completely band limited. Im DFT(Z(@i), 1) The variation of Z( ) with $i in the first term in Eq. (8) & = tan-’ DFT(Z($i), 1)’ Real comes from the factor exp(2c2 sin 8i/2 sin 8, cos(~i 4,)). To investigate the spectral behavior of this factor, where, if k is the number of samples, consider the function g( ) g(x) = exp(cyx) (Y> 0.
Then g(cos 4i) = exp(cr cos $i) and the spectral behavior of the first term can be analyzed by investigating the spectral behavior of g(cos +i) with (Y= 2c2 sin 8i/2 sin 8,. As g(cos $i) is an even function of $i, its Fourier series expansion has only cosine terms. Let
g(COS
DFT(Z($i), 4.2.2.
-bP+l <-
a
p > [4(2C2
sin 8i/2 - I)],
I b 1 -2p+1 = 2c2sin 8i/2 sin en 2ptl * If
where [xl denotes the smallest integer greater than X, then
of Light
Sources for Obtaining
4,,. The number of light sources needed to obtain the solution to & is given by the Nyquist sampling rate for g(cos $i). Denoting the harmonic at which the spectrum of g(cos +i) begins to drop off by p*, we get from the previous analysis
$i) = bo + 2 b, COS(p$i). p=l
In Appendix 1, it is shown that
The Number
m) = $ Z(+il)e-j27im”k.
P * = [&(2C2
sin 8i/2 - l)].
For c = 2.578 and Oi = 40.0”, we get p* = 2. Hence, the spectrum can be expected to start falling after the third harmonic. However, the spectrum falls off much faster than this bound. Figure 4 shows the amplitude spectrum of g(cos 4i) with 8i = 40” and 8, = 20” for different values of c. Note the sharp cutoff. In particular note that for c = 2.578 the function g(cos $i) is essentially band limited to the fourth harmonic and eight or more light sources are sufficient to obtain & from Eq. (9). 4.3.
Solving for en, al, and a2
Having obtained &, consider whether 8,, al, and a2 can be obtained from the continuous function Z(i, n, r).
280
TAGARE AND DeFIGUEIREDO
S,*,@, E*) that give the same Z(i, n, r) for all Xi E [O, I]. This implies that the following equation holds for all Xi E
[O,11: aTexp(c* Xi) - dTexp(c’* Xi) = (6; - at)Xi + h3*- af. (11)
4
2
"0
P
6
8
10
FIG. 4. Amplitude spectrum of g(cos 4,) for different c.
4.3.1. Solution from the Continuous Function Z(i, n, r). Rotate the reference frame in the azimuth plane such that in the new frame c#+,= 0 and all the zenith angles are the same as before. In this frame, Eq. (8) becomes
Z(i, n, r) -L al exp(-2c2 + 2c2 Cos 8i/2 Cos 0,) x exp(2c2 sin 1!&/2cos en cos +i) + a2 sin 8; sin 8, cos +i
The functions on the left and right hand sides of the above equation are linearly independent for Xi E [O, 11. Hence the above equation cannot hold for all Xi E [O, 11. Therefore, two distinct sets (a:, at, a:, c*) and ($, $, 6:, e*) cannot be found to give the same Z(i, n, r) and a?, a?, and a: can be uniquely determined from the continuous function Z(i, n, r). 4.3.2. Solution from Samples of Z(i, n, r). Suppose that light sources are placed so that there are m distinct values for Xi denoted as Xi,, Xi23. . . , Xi,. The intensities in the image due to these lights are I,, Z2, . . . , I,,, given by
+ a2 cos 8i cos 8,
= UTexp(c* Xi) + &Xi + at,
(10)
\Ll
where a; = aI exp(-2c2 + c2 cos 8J2 cos e,), a2*= a2 sin 8i sin e,, af = a2 cos Oi cos e,,
\a: exp(c* Xim)f
U:X, + U3 *I
Suppose that the system of equations (12) cannot be uniquely solved for a?, al, and a:. Then, as before, it is possible to find two distinct sets (a:, a:, a$, c*) and (L$, ti,*, 6:, 6) that give the same left hand side for Eq. (12),
c* = 2c2 sin oil2 cos en, and Xi = COS$+a The variable Xi (and hence Z( )) is an even function of $i in the range (-7~, +r 1. Therefore, we need consider only the range [0, + x] and in this range Xi takes values from 0 to 1. Also, 8,) al, and a2 are related to a?, a:, and a:’ by
LITexp(c*
Xi,)
- hT exp(E* Xi,)
UTexp(c*
Xi*)
- 67 exp(C*
Xi2)
a? exp(c* Xi,) - ciTexp(E* Xi,)
(13) (ti$ -
az*)Xi,
+ ifi: - a$
(ciz*- a:)Xi, + tit - a: =
a1 = exp(-2c2 + c2 cos Bi/2 cos e,y
d
(tit -
a,*)Xim
+ ti: - a3*I.
Equation (13) is Eq. (11) evaluated at Xi,, Xi*>. . . , Xi,. As shown in Appendix 2, Eq. (11) can hold for at most three different values OfXi. Thus, if we set m 2 4, Eq. (13) Thus, if a:, at, and a$ can be uniquely determined from will not hold for some row, thus assuring that Eq. (12) will have a unique solution for a?, at, and a: and thereZ(i, n, r), then so can 8,) al, and a2. Suppose that a?, at, and a: cannot be uniquely deter- fore also for f3,, al, and a2. Due to the symmetry in light mined from Z(i, n, r) as +i goes from 0 to W. Then it is source placement, we need k = 2m light sources to give possible to find two distinct sets (a:, at, a;, c*) and (HT, us m distinct values of Xi. This is shown in Fig. 5. a2 =
sin 8, sin 8i’
SHAPE AND REFLECTANCE
281
FROM STEREO
5.1.
Linear Analysis
The function exp{-c2 [cos-*( )12}can be approximated quite well as a linear function over a large range of its domain. Figure 6 shows the function for c = 2.578 and two linear tits over different ranges. Let us assume that for a particular surface normal n and a set of light sources, the sequence of numbers nfn lies in one of these regions. Then the function 4(x) = exp[-c2{cos-‘
/ Y FIG. 5.
1
I&, n, r) Eight lights for four different
x,.
i
:
\
= al($(xo) - xo~'(al))
Thus, &, being known, a unique solution can be obtained for f&, al, and a2 with k = 8 (m = 4) light sources. 4.4.
:
n
The Total Number of Light Sources
From the above, it is clear that the number of light sources sufficient to guarantee a unique simultaneous estimate of all unknowns is the maximum of eight and the number obtained from Nyquist sampling of Z( 6;). For c = 2.578 the bandwidth of Z(&) is sufficiently small that eight light sources are sufficient for Nyquist sampling. In general, Fig. 4 shows that the number of light sources needed for Nyquist sampling of I( 4;) is not too different from eight and therefore we expect that the number of light sources needed for a unique answer would not differ very much from eight. All the experiments in this paper are conducted using eight light sources. 5.
i
ILL-POSED
NATURE
OF THE JOINT ESTIMATE
The previous section demonstrated that there is a unique answer to the simultaneous estimates and that this estimate can be obtained with a finite number of light sources. However, the analysis of existence and uniqueness of the estimate does not inform us how robust the estimate is with respect to additive noise in the data. In this section, we explore the robustness to noise of the solution to Eq. (6). As we shall see, the solution is sensitive to noise and therefore the simultaneous estimation problem is ill-posed. Since Eq. (6) is nonlinear it is difficult to demonstrate this rigorously. We demonstrate the ill-posed nature in two ways. First, we linearize Eq. (6) and show that the linear equation is ill-posed. Besides establishing grounds for belief in the ill-posedness, the linear analysis also gives us qualitative insight into its origin. Second, we numerically compute the right hand side of Eq. (6) and show that large changes in n cause small changes in the intensity.
where
and
/sin0;,/2 cos &, P=(
sin I$,/2 sin &
cos &,I2
sin ej,12 sin 4,
cos 8,,12
: \
sin 8i,12 COS +ii
sin Oj,Cos $j,
sin 8i, sin C& cos Bi, 1 1
sin Oii Cos &
. .
sin 0,, sin &A cos 0,, 1 9
Note that the vector u has vertical dimension of 7 (each n has a dimension of 3). Under the assumptions that I$,, B;,, . . . ) 81,are in the range [O, 30”], we have the following: 1. cos 8,,/2 and cos 8, are close to 1, thus columns 3, 6 and 7 in the matrix P are close to being linearly dependent, and
282
0.0
TAGARE AND DeFIGUEIREDO
1:
In Fig. 7, note that although the maximum difference between I,( > and Z,( ) is about 40 grey levels, the mean difference is only 6.1 grey levels. Also, the maximum difference occurs at Bi = 0.0”. For nonzero 0i, the difference falls off quite rapidly. The mean difference of 6.1 grey levels is well within the range of noise generated by most cameras. Hence, we conclude that it is difficult to distinguish between the two surface normals and the two “albedos” even by using a large number of light sources in the region 0 I +i < 360” and 0 5 8i i 30”. The problem is clearly ill-posed.
The function \
0.0
x
If0
FIG. 6. Linear approximation to the forescatter lobe function.
2. sin 04/2 = i sin 8,, hence the pairs of columns 1,4 and 2,5 are close to being linearly dependent. Thus the matrix P has three linearly independent columns and hence only three quantities in the vector u can be estimated with any accuracy at all. Since our ultimate goal is to estimate four quantities (a,, a2, &, &J, the problem is ill-posed. Note that the ill-posed nature of the problem is independent of the vertical dimension of P and hence independent of the number of light sources. Qualitatively, the above argument shows why simultaneous estimation is ill-posed: 1. The part of the observed intensity vector that lies along [l 1 . . . llT cannot be easily decomposed into contributions due to the forescatter and normal lobe, i.e., there is no way of telling how much of the “constant” part of the intensity values was generated by either lobe, and, 2. For the range of incidence angles we are interested in, the variation in intensity due to the change in the zenith angle 8ikcannot be uniquely decomposed into contributions due to either lobe. 5.2.
Numerical
Analysis
For the exact nonlinear equation (7) the ill-posed nature is shown in Fig. 7. The figure shows a plot of ]]Zi(&, 4i) - Z2(0i, +i)JI for 0 5 Bi 5 30” and 0 or +i < 360”. The functions Zi( ) and Z,( ) are specific cases of Eq. (7), where the following parameters were used. For Ii ( ), ai = 140.0, a2 = 70.0, c = 2.578, 8, 33.0”, $,, = O.o”, and for
6. REGULARIZATION
OF SIMULTANEOUS ESTIMATES
Ill-posed problems can be solved reliably by regularization. We regularize Eq. (6) by estimating n, aI, and a2 by minimizing
Q= hl(Z- z)T(Z- Z)+ A& - a:>2+ h&2 - ap, (14)
where Z = vector of observed intensities, / exp{-c2[coss’n$]2}
\
I = zmaxcq I
exp{-c?[cos-1n$]2}
/iT n\ +
Zmax~2
I
:
3
i*ii;fn (16)
hi, Xz, A3> 0, and &and c$are prior estimates of (Y~and CX~. This is a nonlinear optimization problem in n, (Y~and
40.0 \
Z2( ), al = 98.0, a2 = 0.0, c = 2.578, 0, = 10.32”, $,, =
0.0. Thus, Zi( ) corresponds to a surface normal with an azimuth angle of zero degrees and a zenith angle of 33.0”, while Z,( ) is for a surface normal with zenith angle of 10.32” and the same azimuth angle. Note that the values of ui, a2 are also different for both the cases (the values used in this example correspond to the actual grey level values recorded by the camera in our experimental setup).
(15)
FIG. 7. Ill-posed nature of the joint estimate
SHAPE AND REFLECTANCE
FROM STEREO
283
FIG. 8. A sphere.
(Yeand a closed form solution cannot be obtained for it. A standard nonlinear, numerical optimization procedure [6] is used to do this. The entire procedure for joint estimation is: 1. The surface is imaged with eight line sources placed at azimuth angles of O”, +45”, -t90”, ? 135”, and 180” at a constant zenith angle of 25”. 2. We know from Section 4 that these light sources are sufficient to obtain an estimate of &. The estimate is obtained via Eq. (9). 3. The other parameters are obtained by minimizing Q in Eq. (14) with the constraint that the value of the azimuth angle & be that given by the above step. 7.
EXPERIMENTAL
and a; = 0.5. The value of I,,, = 140.0 was used in Eq. (15) while (Y: and a:, the prior guesses, were set to 1.1 and 0.45 respectively. The constants in the minimization function of (15) were Xi = 100, h2 = 1600, and A3= 900. As we shall see, this is not an excessive amount of regularization. The initial values for the optimization procedure were n = [0 0 llT, (Y~= 1.2, and a2 = 0.8. Figure 11 shows the reconstruction of the spherical surface.
A spherical surface was fit to the reconstruction and the r.m.s. error between the best fit surface and the reconstructed surface normals was computed. The error was measured as the mean square difference in the angle
RESULTS
Figures 8, 9, and 10 show the three objects that were used in experiments to validate the theory. It was easily determined that for the illumination conditions of this experiment, the reflectance map for the surface would be best described by I,,,,, = 140.0, CY~ = 1.O,
FIG. 9. A cylinder.
FIG. 11. Regularized reconstruction of the sphere
284
TAGARE AND DeFIGUElREDO
FIG. 12. Histogram for (Y,.
between the best fit surface normals and the reconstructed surface normals. This error had a value of 5.21”. Figures 12 and 13 show the histograms of the value of (~1and (~2as returned by the optimization procedure for the entire surface. The peaks in the distribution occur for aciat 1.04 and (~2at 0.34. From the distribution, note that spread of the values is fairly large and corresponds to a mild amount of regularization. It was also observed that for this range of values of o1 and a2 the contribution from the second and third term in Eq. (14) was quite small compared to the contribution from the first term. This further supports the argument that the amount of regularization was not excessive. Finally, Figs. 14 and 1.5show reconstructions for the cylinder and spoon of Figs. 9 and 10. The exact same procedure was used in these reconstructions as was used for the spherical surface. 8. CONCLUSION
The problem of simultaneous estimation of reflectance map and shape from photometric stereo was posed as a series of three questions about the existence, uniqueness, and robustness of the estimates. We investigated simultaneous estimation for diffuse surfaces whose lobes of the reflectance map were not easily separable. In the first
w
1.0
FIG. 13. Histogram for (Y?.
FIG. 14. Regularized reconstruction of the cylinder.
part of the paper, we showed that a unique answer was available to the simultaneous estimates. Further, it was possible to show that the azimuth angle of the surface normal could be estimated independently of the rest of the unknowns. One of the main contributions of this paper is to establish the number of light sources sufficient to obtain the unique estimate. By an approximate theoretical analysis, it was shown that a small number of light sources yield estimates, but that the estimates are sensitive to noise in the data. From a linear analysis we argue that the sensitivity to noise could not be reduced by adding extra light sources and that prior estimates of some of the parameters of the map are required to regularize the answer.
\
2:o
FIG. 15. Regularized reconstruction of the spoon.
SHAPE AND REFLECTANCE
Experiment with the regularized simultaneous estimation procedure yields good results. Although most of the arguments in this paper are developed in the context of a specific reflectance map, the techniques of Appendix 1 and 2 are fairly general and can be easily used to make similar arguments about other reflectance maps.
285
FROM STEREO
Repeatedly integrating by parts, and using the fact that for 1 I m 5 p - 1, (d”ldy”){(l - Y~)“-“~} = 0 at y = * 1, we get
1’_,g(y)$, ((1- Y*Y-I’*)
dv
= (-l)P j;, g’p’(y)( 1 - Y*)P-“~ dy, APPENDIX
1: THE FOURIER
SERIES FOR g(COS+)
where g@‘(y) denotes the pth derivative of g( ). Thus,
Let g(COS 4i) = bo + C bp COS p$i.
(-1)2P 2pr(p + $>i ‘, (1 -
bp = G
P
Then, for p 2 1,
and
bP= i I,‘”gtcos +i)COS p4i
d4i
b
rep
-= P+’
=-
i
I,” dcos
$i> COS p+i
d$i*
bP=5 j:,
g(y)
VT%
dy,
1
_
2~ + 1
bP
T (y> = (-l)P e P 2p+’ r(p
V’& dP dyp {(I - y*)p-“*}. + 3)
+)
JLl
(I
-
1s”
dP dy”
_
2pcl+l
lh
< -
=-
Evaluating the integral by parts,
1’_,g(y)$ ((1- Y’)~-“~) -j&
a
%+I
=-I
I’ _] g’(y)
(17)
(1
y*)(l
-
(1
I.f-‘I
(1
y*>(l
-
(1 -
-
y2)p-“*
(1
-
-
y21p-“*
Y*)~-“~
y21p-“*
y2Y-“*
s@+“(r) P(r)
P(y>
,@“(u>
P’(y)
41 4
&I 41
&I
/j-i1 (1 - Y~)~-“~ gck)(y) dyl
a
2p + I
= 2c2 sin Oil2 sin 0, 2p+l ’ dv
((1 - y2)“-“2})I, APPENDIX
-
dy’
J-1, (1 - y2)p+‘” g@+“(y) dy j-L1 (1 - y2)p-“* g’p’(y) dy ’
1s”
x ((1 - y*)~-“~) dy
x ((1 - Y~)P-“~} dy.
dy
J’, (1 - y*)~-“~ g’“‘(y)
lj-$
(-1)~ m 6 2”+’ I’(p + +)
g’p+‘)(y)
Since g(y) = exp(cuy), we have g(p+“(y) = ag’p’(y) and since in our case (Y > 0, we have g@‘(y) > 0 and g@+‘)(y) > 0. Further, in the range of integration the term (1 - y*) is always between 0 and 1. Therefore, from Eq. (17) we get
Substituting this into the above equation, g(y)
Y*)P+“~
From this, we get
bP+’ -_-~
where TJy) = cos(p cos-’ y) is the Tchebycheff polynomial of degree p. T,( >are a set of orthonormal polynomials, and are generated by Rodrigues’ formula.
+
2Up + 8)
4.J
Changing the variable of integration to y = cos $i, we get
= g(y)
g@) (Y) dy
Y~)~-“~
s
‘, ‘F’(Y) g’
{(I - Y~)“-“~} dy
2: THE NUMBER OF SOLUTIONS EQUATION (11)
TO
Note that the equation ((1 - y’)p-I’*}
dy. (~1
exp(clx)
-
(~2
exp(c0)
= plx + p2
(18)
286
TAGARE AND DeFlGUElREDO
can have at most three solutions for x given any values of al, QI~,pi, p2 and cl, c2. This is seen as follows: Denote the left hand side of Eq. (18) byf(x) and rewrite it as f(x)
-
PIX
+
P2 =
(19)
0.
Ifx,, x2, . . . ) x, are the solutions of Eq. (18) placed in ascending order, then, between any two consecutive solutions, the derivative of the left hand side of Eq. (19) goes to zero at least once, i.e., if the number of solutions to f’(x)
-
Pl
=
(20)
0
is m, then n I m + 1. Using the same argument on Eq. (20) we have that if f”(x)
= 0
(21)
has k solutions, then m I k + 1. Hence, n I k + 2. But Eq. (21) can be written as cz,c: exp(c,X) = u2c: exp(czx). Taking logarithms, clx = c2x + log (Y~c:- log (Y,c: which is a linear equation and has a single solution. Hence, k = 1 and n = 2. Thus Eq. (18) has at most three roots. (It is easy to find (~1, 02, cl, c2, j31, ~32that cause Eq. (18) to actually have three roots.) Referring back to Eq. (1 l), suppose that for four or more different values of xi, there are two different sets (a:, a:, a$, c”) and (@, @, b:, C*) that give identical values of Z(xJ for some values Of Xi. Then, from Eq. (lo), for all i = 1, 2, 3, 4, . . . , UTexp(c$xJ +
U2*xi
+ U: = LjTeXp(EfXi)
+ Lif,Xi+ hf,
i.e., UT eXp(CfXi)
-
Lif eXp(E*Xi)
=
($
-
Uz*) Xi +
i=
b3* -
Ut,
1,. . . ,4.
This has the same form as Eq. (18) for which there cannot be more that three solutions. Hence, Eq. (11) cannot have more than three solutions. ACKNOWLEDGMENTS We gratefully acknowledge the help we received in regard to Appendix 1 from Dr. Don H. Johnson of the Department of Electrical and Computer Engineering of Rice University.
REFERENCES 1. P. Beckmann and A. Spizzichino, The Scattering of Electromagnetic Wavesfrom Rough Surfaces. Pergamon, Elmsford, NY, 1963. 2. R. C. Birkebak and E. R. G. Eckert, Effect of roughness of metal surfaces on angular distribution of monochromatic reflected light, J. Heat Transfer, February 1965, 85-94. 3. S. Chandrasekhar, Radiative Transfer, Dover, New York, 1960. 4. E. N. Coleman and R. Jain, Obtaining 3-dimensional shape of textured and specular surfaces using four-source photometry, Comput. Graphics Image Process.
18, 1982,309-328.
5. B. K. P. Horn, Robot Vision. M.I.T. Press, Cambridge, MA, 1986. 6. IMSL, ZMSL Math Library, Houston, TX. I. G. Kortum, Reflectance Spectroscopy, Springer-Verlag, Berlin/ New York, 1969. 8. S. K. Nayar, K. Ikeuchi, and T. Kanade, Extracting Shape and Reflectance of Lambertian, Specular and Hybrid Surfaces, Technical Report CMU-RI-TR-88-14, C.M.U., August 1988. 9. W. M. Silver, Determining Shape and Reflectance Using Multiple Images, master’s thesis, MIT., 1980. 10. J. M. Stone, Radiation and Optics, McGraw-Hill, New York, 1963. 11. H. D. Tagare, A Theory of Photometric Stereo for a General Class of Reflectance Maps, Ph.D. thesis, Rice University, Houston, Texas, 1989. 12. H. D. Tagare and R. J. P. deFigueiredo, A framework for the construction of general reflectance maps for machine vision, C.V.G.Z.P., to appear. 13. H. D. Tagare and R. J. P. deFigueiredo. A theory of photometric stereo for diffuse, non-Lambertian surfaces, Z.E.E.E. Trans. Puttern An. Mach. Intelligence, 1991. 14. K. E. Torrance and E. M. Sparrow, Theory of off-specular retlection from roughened surfaces, J. Opt. Sot. Am. 57, 1967, 11051114.