Shape and Source from Shading

I. Introduction

Given an image, E, and a reflectance map, R, the shape­ from-shading problern may be regarded as that of recover­ ing a smooth surface, z(x, y), satisfying the image irradiance equation

over the domain 0 of E. Any given boundary conditions on z( x, y} should also be satisfied. The problern takes the form of a first-order partial differential equation (Horn, 1970 & 1975).

Implicit in this formulation are a number of assump­ tions, the principal one being that the brightness of a sur­ face patch does not depend on its position in space. Another is that an image depicts a smooth surface of homogeneaus reflectance. Several algorithms have been devised to tackle the problem, notably those of Horn (1975), Strat (1979), and Ikeuchi & Horn (1981).

*

On leave from:

Flinders University of South Australia,

Bedford Park, SA 5042, Australia.

SHAPE AND SOURCE FROM SHADING

Michael J. Brooks* Berthold K.P. Horn

Artificial Intelligence Labaratory Massachusetts Institute of Technology Cambridge, Massachusetts, 02139, U.S.A.

One of the many difficulties these schemes face in prac­ tice is that the reflectance map is typically not known. The reflectance map specifies how the brightness of a surface patch depends on its orientation under given circumstances. It therefore encodes information about the reflecting prop­ erties of the surface, and information about the distribution and intensity of the light sources. In fact , the reflectance map can be computed from the bidirectional reflectance­ distribution function and the light source arrangement, as shown by Horn & Sjoberg (1979) .

When encountering a new scene, the information re­ quired to determine the reflectance map is usually not avail­ able. Yet without this information, the shape-from-shading problern cannot be formulated , much less solved. We may resolve this dilemma if a calibration object of known shape appears in the scene, since the reflectance map can be com­ puted from its image. Here we wish to consider the situa­ tion where we are not that fortunate.

It is interesting to evaluate how some basic assumptions can resolve this impasse. Pentland (1984} has looked at the problern of recovering shape from shading under the assumption that the image depicts a Lambertian surface illuminated by a point source whose direction is unknown. Under the additional assumption that the surface is locally umbilical, surface normals are shown to be recoverable by a local operation. This method does not depend on the iterative propagation of information across the image.

There are some serious drawbacks to the local approach, however. One problern is that the umbilical assumption is very restrictive. In fact, spheres are the only surfaces whose points are all umbilical. So this method naturally computes incorrect normals for other shapes, although the errors for approximately spherical surfaces, such as ellip­ soids of low eccentricity, may be acceptable. Further, the constraining effect of known occluding boundary normals cannot be incorporated into the local method. This is un­ fortunate because these normals provide powerful boundary conditions on the shape-from-shading problem, as shown by Bruss (1983}. Finally, because a local method does not take into account neighbors when calculating the normal at a point, nearby normals may differ a great deal, particularly in the presence of noise.

We now present an alternative approach that does not suffer from these disadvantages.

ß. An iterative scheme for shape and source

The task is to recover the shape of a smooth surface depicted in an image, E, that is defined over a region 0 in the xy-plane. Let the shape of the surface be characterized by the function, n, that associates a unit normal with each point in 0. The problern is therefore to find n{x, y) over 0. Assurne for now that the object has a Lambertian surface, and that it is illuminated by a single point source. If the vector s points to the source, and n{x, y) is the unit nor­ mal of a surface patch, then the apparent brightness of the patch is given by the reflectance map

We do not, by the way, force 8 to be a unit vector; this allows for the possibility that the intensity or the source may be unknown. This way we can also deal with unknown sensor sensitivity and unknown surface albedo, provided it is uniform.

Our problern is now to find a smooth shape, n, and source direction, 8, satisfying the image irradiance equation

In practice, brightness cannot be determined with per­ fect accuracy, and so it appears reasonable to transform this into a minimization problern {lkeuchi & Horn, 1981; Horn & Brooks, 1985). There is another reason to con­ sider this as a minimization problem: if we simply try to solve the image irradiance equation as it stands, we obtain a set of differential equations equivalent to the character­ istic strip equations. Here, however, we seek an iterative scheme Jending itself to a parallel implementation on a grid, as originally suggested by Horn {1970). Further, a shape­ from-shading problern that has noisy image data may weil not have a theoretical solution. A minimization approach will, however, enable the recovery of a shape that fits the given data best, in a sense determined by the functional chosen.

We seek a smooth shape, n, and a source direction, 8,

If a solution exists, and there are no errors in brightness measurements, then the image irradiance equation will have been satisfied (although there is no guarantee that the re­ sulting n will be integrable; see Horn & Brooks, 1985). We adopt a regularizing component (Poggio & Torre, 1984) by incorporating the expression

which is intended to select a particularly smooth solution from a possibly infinite set of candidates. Note that a sub­ script here denotes partial differentiation, and that squaring a vector is equivaJent to taking the dot-product with itself. Finally, we wish to insist that normals have unit length. This is accomplished with the constraint

Combining the three terms gives the composite functional

which is to be minimized with respect to n and 8. Here, >. is a scalar that weights the relative importance of the regularization term, while ~J(x, y) is a Lagrangian multiplier function used to impose the constraint that n(x, y) be a unit vector (see Horn & Brooks, 1985).

Minimizing I is a problern in the calculus of variations. First assume that s is known and that minimization is to be achieved by a suitable choice of n . Extrema of tunetionals are obtained by examining solutions to the associated Euler equations (see Courant & Hilbert, 1953). The functional

has the Euler equation

So, by substitution, it follows that I has the Euler equation

where

Now, a discrete approximation to the Laplacian operator is given by

in which E' is the distance between adjacent picture cells in the image, and Dij is the local average of normals

Hence we may translate the Euler equation into the discrete form

Rearranging this in order to isolate n;j on one side yields the iterative scheme

which computes shape given direction of the light source. Other approximations for the Laplacian may Iead to im­ proved results, at the cost of increased computation. For example, if we use the more accurate 9-point approximation for the Laplacian, in which

then twice as many array accesses are needed (and the con­ stant multiplier E" 2 /4>. becomes 3E'2 /10>.). The simple 5point approximation was adequate for our purposes.

Note that we have yet to solve for IJ, the Lagrangian multiplier. This can be avoided, however, by observing that the division of the right band side by {1 + ~J(E' 2 /4>.)) does not change the direction of the vector being computed.

Since 1-' is intended to ensure that the result is normalized, we can simply do this explicitly, as in

Now consider the problern of minimizing I with respect to s, given that n is known. This is a problern in conven­ tional calculus. Computing the partial derivative of I with respect to s, we have

and so

Noting that

then by substitution we have

where (nnT) , and also the integral of (nnT), are 3 x 3 ma­ trices. From this we finally obtain the desired equation

Here -l denotes the inverse of a matrix. A discrete version of this formula, in which the integrals are replaced by sums, is easily obtained. An iterative scheme in both n and s is then within grasp:

This takes advantage of the fact that s need not be a unit vector. If the brightness of the source is known , we can normalize s so that it is a unit vector. Then the determi­ nation of s becomes slightly more complex, since it involves a constrained minimization.

111. Properties and performance of the scheme

The iterative scheme has two components: one con­ cerned with the recovery of shape, the other concerned with the determination of the source direction. The shape­ recovery component has an intuitively satisfying form. In essence, a new normal is computed by taking a local av­ erage, and adjusting this either toward or away from the source. The magnitude and sign of the adjustment is de­ termined by the brightness error of the current estimate.

For a given shape, a new source direction is computed by a single pass through the image; unlike shape-recovery, no iteration is necessary. The 3x 3 matrix. (nnT) is summed across the image, as is the vector (En). The source direc­ tion can then be computed using Gaussian elimination or even Cramer's method (see Korn and Korn, 1Q68). The source-recovery component has been tested on a number of images and shapes. When the data are free of noise, the estimate of source direction is extremely accurate. Fur­ thermore, estimates remain very good in the face of signif­ icant noise. For example, a synthetic image was generated of a sphere illuminated by a point source in the direction ( -4, 3, 8)T. The image was quantized at 255 irradiance Iev­ els, and the correct surface normal was given for each of the 1250 image points. Gaussian noise was added to the image giving an average perturbation in irradiance values of 34. Despite this, the source-finder computed an estimate of source direction that was only 2.7° in error. Further trials gave similar results.

The source-direction estimates are robust because the whole image is used. Theoretically, the problern is highly over-determined as source· direction can be recovered from brightness values of only three different surface orienta­ tions. Using the whole image ensures, however, that noise effects are significantly reduced.

The shape-and-source-from-shading problern for a point light source and Lambertian surface has a natural twcrway ambiguity. If the image irradiance equation is satisfied over 0 by the shape n1 and the source direction s1, it will also be satisfied by the dual shape n2 and source direction s2 where

and z is the viewing direction. Here, both source direc­ tion and surface normals are reflected about . the viewing direction. This is easily verified by observing that

Given some initial values for the normals, the shape-and­ source scheme will head for one or the other of these scr lutions. The dual shape and source direction can then be determined immediately using the equations given above.

We now present two examples of the program at work. Synthetic images were used, each of which depicted a Lam­ bertian surface illuminated by a point source in the di­ rection (3, 2, Q)T. The images each contained more than 1000 points at which normals were to be determined. Nor­ mals were assigned an initial value of (0, 0, l)T , as was the source direction. Occluding boundary normals were given. The equations for n;j could be solved sequentially using the Gauss-Seidel algorithm. Since we are ultimately inter­ ested in parallel implementation on a grid, we used the Ja­ cobi method instead ( despite the fact that the Gauss-Seidel method has slightly better convergence properties).

The first image portrayed a hemisphere viewed from directly overhead. After 100 iterations with A = 0.005, the average angular difference between estimated and correct normals was less than 3°. The maximum such deviation was less than 2.5 times the average value. The estimate of the light source direction, at this time, bad errors in azimuth and zenith angles of 1.4° and 1.6° respectively.

The second image depicted a cylinder with rounded, hemispherical ends, viewed from a direction perpendicular to its axis. After 60 iterations, this time with A = 0.003, the average angular error in surface-normal was less than 5°. A further 30 iterations brought this value down to 4°. The maximum error remained somewhat !arger, however, due to the scheme's tendency to smooth the intersections of the cylinder and the hemispheres. The errors in azimuth and zenith angles for the source were 7.3° and 1.1° respec­ tively, achieved after QO iterations. These, too, improved slowly with further processing.

The scheme was sometimes slow in converging. After rapid initial improvements, the rate of progress would de­ crease appreciably. However , one might expect the scheme to be slower than some of the current iterative methods given the disadvantage of not knowing the light source di­ rection. Convergence could be accelerated by employing multigrid techniques that propagate information across the image more quickly {see Terzopoulos, 1Q84). Interestingly, in the examples considered, a reasonable estimate for the source direction was obtained after only a few iterations. Subsequent processing just improved the estimate.

IV. Iterative schemes for other reßectance maps

We now present two more new iterative schemes: the first extends the shape-and-source finder to cover situations in which a simple model of the sky is also included; the sec­ ond uses methods developed above to find shape from shad­ ing, given a general reflectance map, but does not recover source direction.

A. Incorporating a sky component

The reflectance map

captures the situation in which a Lambertian surface is il­ luminated by a hemispherical source of uniform radiance (Brooks, 1Q78; Horn & Sjoberg, 1Q7Q). A point source may be added to the map to give

Here, a controls the relative intensity of the sun and the sky. We can now generate a method of shape-and-source recovery, under the assumption that the imagewas formed in accordance with the reflectance map R,. .

We seek to minimize

with respect to both n and 8. Fixing 8 for the time being, we are required to minimize the above functional with respect to n alone. The Euler equation for this problern is

Treating J.L as before, the following scheme is obtained:

Here, the reflectance map, R.,, has been substituted back into the equation to improve the presentation.

We now assume n to be fixed and minimize the func­ tional with respect to 8. This we do, as before, by differ­ entiating with respect to 8 and equating the result to zero. Thus we have

Expanding,

Noting as before that (n.s)n = (nnT)s, the equation be­ comes

Thus we obtain the equation in 8 given by

This we may write in discrete form and combine with the iterative scheme for n derived previously to give

Note that a is assumed to be known. Interestingly, the computation of 8 proceeds as before, except that the con­ tribution of the sky is subtracted from E. This does not render the scheme trivial, however, as the calculation of shape does not follow suit.

B. Recovery of 8hape for the general reßectance map

Recall that the iterative scheme for shape and source direction is composed of two parts. Either component of the scheme may stand alone in the event that shape is required from source, or source is required from shape. Indeed, if used in this way, the shape-recovery component may be generalized to incorporate any reflectance map, R(n) . In minimizing the functional

we obtain the Euler equation

from which we derive the scheme

This is perhaps the most appealing of the current shape­ from-shading schemes that deal with a general reflectance map. It is simply derived, and is expressed elegantly in terms of unit-normals, rather than a two-parameter system such as the stereographic fg space of Ikeuchi and Horn.

V. Summary

Most current methods for obtaining shape from shading assume complete knowledge of the reflectance map. Here, we considered the situation in which a Lambertian surface is illuminated by a point light source from an unknown di­ rection. Thus we dealt with a parameterized, rather than a fixed, reflectance map. The local approach to recovering shape from shading, which is also intended to deal with unknown source direction, was found to have several draw­ backs, notably its restrictive assumption that surfaces are locally umbilical.

The adoption of unit normal vectors for describing sur­ face orientation was important to the development of our method. It led to simple derivations and elegant presen­ tations. The problern was cast as one of minimizing a positive-definite functional incorporating a brightness er­ ror term, a regularizing term, and a Lagrangian multiplier to enforce the condition that the normal be of unit length. The Euler equation for this calculus of variations problern was shown to be a second-order partial differential equation in the unknown surface-normal function. A convergent it­ erative scheme solved it in the discrete domain.

The direction of the light source can be determined in closed form if the surface shape is known. At any iter­ ation, a source-direction estimate can be obtained using the current estimate of the surface shape. The iterations for obtaining increasingly accurate estimates of the surface shape can be interlaced with estimation of the light-source direction.

We implemented and tested this method for recovering shape and source direction. We also discussed a two-way ambiguity that can appear in the solution. Further, we showed how to extend the shape-from-shading component of the iterative scheme to more general reflectance maps.

Acknowledgements

We thank Alan Yuille for many useful discussions. Com­ ments on a draft by Steve Bagley, Eric Grimson, Tomas Lozano-Perez, Bruce Nelson and, especially, Demetri Ter­ zopoulos were much appreciated. Mike Brooks wishes to express gratitude to the MIT AI laboratory for the oppor­ tunity to spend some months there, and to Flinders Univer­ sity for granting leave. Support for the laboratory's artifi­ cial intelligence research is provided in part by the System Development Foundation.