Shape from Shading

Shape from shading recovers the shape of a surface from the way its brightness varies across a single image. A matte surface looks brighter where it turns toward the light and darker where it turns away, so the shading carries information about surface orientation at every pixel, even on a surface with no texture or markings. Horn and Brooks (1989) describe the task as an “inverse graphics” problem: rendering shading from a known shape is simple, while recovering the shape from the shading is much harder: each pixel gives one brightness measurement for an orientation with two degrees of freedom, so additional constraint, such as the assumption that the surface is continuous and smooth, must be brought in (Horn and Brooks, 1989). In planetary science the same problem, recovering slopes from the brightness of single images, is called photoclinometry; Horn (1990) notes that astrogeologists and machine vision researchers developed it independently and had become aware of each other’s work only a few years before he wrote.

Shape from shading inverts the shading part of image formation from one image; photometric stereo adds lights until each pixel is determined on its own, while shape from shading replaces the missing measurements with assumptions about the surface, the boundary of the region and the lighting. It is one of the cues for depth estimation.

Problem Definition

Assumptions. The classical formulation, as Horn (1990, Section 2) reviews it, assumes that the depth range of the scene is small compared with its distance from the viewer, so the projection is approximately orthographic; that the light sources are far away, so every surface element receives the same irradiance from the same direction; and that the surface is uniform in its reflecting properties. Under these assumptions the radiance of a surface patch depends only on its orientation, not on its position in space. Image irradiance is proportional to scene radiance (for a narrow field of view, up to a factor that depends on the f-number), and the constant, together with the source strength and the albedo, is removed by normalizing the image, for example by the brightness of a point whose orientation is known (Horn, 1990, Section 2.4).

Surface orientation. The surface is a height function $z(x, y)$ above a reference plane perpendicular to the optical axis, with $z$ measured toward the viewer. Its gradient components are $p = \partial z / \partial x$ and $q = \partial z / \partial y$, and the unit normal is

$$\mathbf{n} = \frac{(-p,\; -q,\; 1)^\top}{\sqrt{1 + p^2 + q^2}},$$

with the direction toward the viewer $\mathbf{v} = (0, 0, 1)^\top$ (Horn, 1990, Eq. 6; Horn and Brooks, 1989), the convention of the example in photometric stereo. Szeliski (2022, Eq. 13.2) writes $(p, q, 1)$ with the same definition of $(p, q)$; the formulas below hold in Horn’s convention. Other parametrizations are in use: the unit normal itself, slant and tilt, and the stereographic coordinates $(f, g)$ of Ikeuchi and Horn (1981), described in variational shape from shading.

The reflectance map. Because radiance depends only on orientation, it can be written as a function $R(p, q)$ on gradient space, the reflectance map, usually drawn as contours of constant brightness in the $pq$-plane (Horn, 1990, Section 2.3). It can be measured with a sample on a goniometer stage, read off the image of a calibration object of known shape such as a sphere, or derived from a reflectance model by integrating the bidirectional reflectance distribution function over the distribution of source brightness. For an ideal Lambertian surface lit by a distant point source in the direction $\mathbf{s} \propto (-p_s, -q_s, 1)^\top$, brightness is proportional to the cosine of the incident angle, $\mathbf{n}^\top\mathbf{s}$, so that

$$R(p, q) = \frac{1 + p_s\,p + q_s\,q}{\sqrt{1 + p^2 + q^2}\,\sqrt{1 + p_s^2 + q_s^2}}$$

where the numerator is positive, and $R = 0$ otherwise (Horn, 1990, Eqs. 8–9).

The image irradiance equation. Shape from shading asks for a surface whose predicted brightness matches the image at every point:

$$E(x, y) = R\big(p(x, y),\, q(x, y)\big).$$

Since $p$ and $q$ are the first partial derivatives of the unknown $z$, this is a first-order partial differential equation, and nonlinear for most reflectance maps (Horn, 1990, Section 2.4). The task is to find $z(x, y)$ given $E(x, y)$ and $R(p, q)$; additional constraints, such as boundary conditions and singular points, are needed for the solution to be unique (Horn, 1990, Section 4.4). Durou, Falcone and Sagona (2008) write the Lambertian case with a unit light vector $\boldsymbol{\omega}$ as a Hamilton–Jacobi equation,

$$I(x)\sqrt{1 + |\nabla u(x)|^2} + (\omega_1, \omega_2)\cdot\nabla u(x) - \omega_3 = 0,$$

where $u$ is the height and $I$ the normalized grey level. For a frontal light, $\boldsymbol{\omega} = (0, 0, 1)$, it becomes the eikonal equation $|\nabla u(x)| = \sqrt{1/I(x)^2 - 1}$, the form that most papers use.

A reflectance map linear in the gradient. Viewed from a great distance, the material of the lunar maria has a brightness that depends only on the ratio of the cosines of the incident and emittance angles, which is linear in $p$ and $q$ (Horn, 1990, Section 2.5; Horn and Brooks, 1989). With $R(p, q) = f(c\,p + s\,q)$ and an invertible $f$, the slope in the image direction $(c, s)$ is known from the brightness at each point, and profiles of the surface follow by integrating it along straight lines. This was the first shape from shading problem solved other than heuristically, by Rindfleisch (1966) for lunar topography, in a formulation with full perspective projection (Horn, 1990, Section 2.5). Durou, Falcone and Sagona (2008) trace the first mention of 3D reconstruction from photometric cues to the astronomer Van Diggelen (1951), and credit Horn with the name “shape from shading” and with showing that Rindfleisch’s particular case could be generalized by characteristic strip expansion.

Input and Output

• Input: one image whose values are proportional to scene radiance (raw or radiometrically linearized, for the reasons given in image formation); the reflectance map, which means the light direction and the reflectance of the surface, given or estimated; the region of the image to reconstruct; and, where available, boundary information: heights or orientations on the region’s boundary, the occluding boundary of the object, or the heights at singular points.
• Output: a field of surface orientations (Horn calls it a needle diagram), or heights up to an additive constant, which form a depth map of the visible surface. Methods that estimate $(p, q)$ first need a second step to find a surface: the least-squares fit of $z$ to a gradient field, minimizing $\iint (z_x - p)^2 + (z_y - q)^2\,dx\,dy$, leads to the Poisson equation $\nabla^2 z = p_x + q_y$ with the natural boundary condition that the normal derivative of $z$ matches the slope given by $(p, q)$ (Horn, 1990, Section 4.5; Horn and Brooks, 1986). Frankot and Chellappa instead project the estimated gradient onto the nearest integrable one in an orthogonal (Fourier) basis (Horn and Brooks, 1989).

A gradient field is integrable if it is the gradient of some surface, which requires the mixed second derivatives to agree, $p_y = q_x$ (Horn and Brooks, 1989). Shape from shading methods differ in whether they enforce this, and the difference matters for both accuracy and well-posedness.

Variants

• Unknown light source. Most methods assume a known source direction. Brooks and Horn estimated shape and source direction alternately, and Pentland (1984) proposed estimating a single source direction from the statistics of brightness derivatives (both as Horn and Brooks, 1989, describe them). On a real satellite image, Horn (1990, Section 7.3) found that Pentland’s estimator and Lee and Rosenfeld’s refinement of it failed to give the source azimuth, and that an azimuth found from the distribution of brightness gradients was ambiguous between two directions, corresponding to the convex and concave interpretations of the surface.
• Unknown albedo. A system that assumes uniform reflecting properties where they vary recovers the wrong shape; Horn and Brooks (1989) give countershading on animals and make-up as examples. With both lighting and albedo unknown, the generalized bas-relief ambiguity of photometric stereo applies.
• Perspective projection and nearby lights. Horn’s original work (1970, 1975) did not assume a distant viewer or distant sources, which, as Horn and Brooks (1989) remark, belies the impression that the problem can only be handled orthographically. Durou, Falcone and Sagona (2008) report that the equation of perspective shape from shading was established almost simultaneously by three papers around 2003–2004, and that some models put the light source at the centre of projection, as with an endoscope. Prados and Faugeras (2005) showed that such a model becomes well-posed when the $1/r^2$ fall-off of the light is kept (Approaches, below).
• Other reflectance and lighting. Non-Lambertian reflectance, extended or diffuse lighting and interreflections each change the image irradiance equation (Durou, Falcone and Sagona, 2008, Section 2.1). Collimated and extended sources are accommodated in the reflectance map by integrating the reflectance function against the source distribution (Horn, 1990, Section 3.1).
• The shape-from-X family. Szeliski (2022, Section 13.1) groups shape from shading with shape from texture and shape from focus (see 3D reconstruction for the other cues). Photoclinometry developed profile methods that assume a locally cylindrical surface or a zero cross-track slope, and radarclinometry applies the same ideas to synthetic aperture radar images (Horn, 1990, Section 3; Horn and Brooks, 1989).

Challenges

One equation, two unknowns. At a general image point, brightness constrains the gradient to a contour of the reflectance map; it does not fix it. Horn and Brooks (1989) point out that the brightness pattern of an arbitrary image region could in general come from infinitely many different surfaces.

Ambiguity along the dark direction. The linear reflectance map shows the ambiguity most clearly. If $z(x, y)$ is a solution of $f(c\,p + s\,q) = E$, so is

$$\tilde z(x, y) = z(x, y) + g(s\,x - c\,y)$$

for an arbitrary function $g$, because adding $g$ does not change the slope in the direction $(c, s)$ (Horn, 1990, Eqs. 18–20). Each profile has its own constant of integration, and only an initial curve from which the profiles start removes the ambiguity. Horn adds that the same effect governs iterative methods under oblique light: errors in the computed gradient are likely to be small in the direction toward or away from the source and large at right angles to it.

The convex–concave ambiguity. Under a frontal light the eikonal equation depends on $|\nabla u|$ only, so a surface and its reflection in depth, $z$ and $-z$, produce the same image. Belhumeur, Kriegman and Yuille (1999) state the general case for Lambertian surfaces: the inverted surface, lit by a correspondingly transformed source, produces the same illuminated regions, which they call the well-known up/down or convex/concave ambiguity, although the cast shadows of the two surfaces may differ. For a single distant light it follows from the normal formula (our derivation) that $\mathbf{n}(-z)^\top \mathbf{l}' = \mathbf{n}(z)^\top \mathbf{l}$ when $\mathbf{l}' = (-l_x, -l_y, l_z)$: inverting the relief and rotating the light by 180° about the viewing axis leaves every pixel unchanged. Prados and Faugeras (2005, Fig. 1) illustrate it with the crater illusion, a photograph of two ash cones that look like craters until one imagines the light coming from the other side.

Boundary conditions do not settle it. Durou, Falcone and Sagona (2008, Section 3.1.4) stress that adding boundary conditions does not resolve the concave/convex ambiguity, and that in practice boundary conditions are seldom known. The Dirichlet problem for the eikonal equation can have several solutions in the weak (viscosity) sense as well as several classical ones: a smooth surface and the non-smooth surfaces obtained from it by reflecting parts of it all fit the same image and the same boundary (their Fig. 2). A result of Ishii and Ramaswamy, as they summarize it, guarantees a unique maximal viscosity solution when the image is continuous and the number of singular points is finite, but not when the singular points cover a region of positive area, such as a flat roof; and the maximal solution need not be the true surface.

Singular points and occluding boundaries. Two kinds of image points do fix the orientation directly. At a singular point, where the reflectance map has its unique maximum (for a Lambertian surface and a point source, where the normal points at the source), the gradient is known; singular points are sources and sinks of the characteristic curves and strongly constrain the solution (Horn, 1990, Section 4.2; Horn and Brooks, 1989). On the occluding boundary (the limb), the surface is seen edge-on: its normal is perpendicular to the viewing direction and parallel to the normal of the silhouette in the image (Horn and Brooks, 1989; Ikeuchi and Horn, 1981). There the gradient is infinite, which is why Ikeuchi and Horn (1981) moved to stereographic coordinates. Horn (1990) notes that it had recently become questionable how much constraint occluding boundaries add, given that singular points already constrain the solution strongly.

Is the problem ill-posed? Horn and Brooks (1989) argue that the problem is ill-posed for an image patch that contains no singular point and no limb, but not necessarily when they are included, much as a linear partial differential equation has a unique solution only with appropriate boundary conditions. In their account the general belief that shape from shading is ill-posed comes from methods that treat $p$ and $q$ as unrelated and so need regularization. Uniqueness has been proved only in special cases. As Horn and Brooks (1989) summarize Bruss’s result, for a rotationally symmetric reflectance map the solution is unique apart from a reversal of depth if the image has a single singular point, the height vanishes to second order there, and the smooth closed silhouette is the limb of the object.

Existence. Some shaded images correspond to no surface at all. For a Lambertian surface lit from near the viewer, no surface with continuous first derivatives produces a simply connected dark region in the middle of a bright planar region (Horn, Szeliski and Yuille, 1989, as Horn, 1990, Section 7.2 reports); on such an image Horn’s algorithm returned a “solution” with a conical singularity. A real image need not have an exact solution either, because of aliasing introduced in sampling or because no homogeneous surface with the assumed reflectance and lighting could have produced it (Horn, 1990, Section 5.1).

Low spatial frequencies. Shading depends on the first derivatives of height, so slowly varying components of the shape have little effect on the image, and the largest errors of a reconstruction are expected at low spatial frequencies (Horn and Brooks, 1989; Horn, 1990, Section 6.6). Horn and Brooks note that people, too, appear to recover rapid undulations well and to miss or misread slow ones.

What makes an instance hard. From experiments with synthetic images, Horn (1990, Section 7.4) found that surfaces with only low slopes (bas-relievo) are easy compared with steep ones (alto-relievo); that oblique illumination is easier than a light at the viewer; that problems where the reflectance map is nearly linear in the gradient are easier than those where its contours are strongly curved; that smooth, especially convex, surfaces are easier than wrinkled ones; and that discontinuities in orientation are harder still. When some patches are turned so far that their brightness falls again, a patch can be dark either because it is not turned enough toward the light or because it is turned too far, a local two-way ambiguity that traps simple schemes in local minima.

Approaches

Durou, Falcone and Sagona (2008) sort numerical methods into three classes, following a very similar classification by Zhang et al. (1999): methods that solve the partial differential equation, methods that minimize a functional, and methods that approximate the image irradiance equation, which Zhang et al. divide into local and linear methods. The history runs roughly in that order.

Characteristic strips

Horn’s original solution (1970, 1975) used the method of characteristic strip expansion for first-order partial differential equations. The image irradiance equation is replaced by five ordinary differential equations along a curve with parameter $s$ (Ikeuchi and Horn, 1981, Section 1.2):

$$\dot x = R_p, \quad \dot y = R_q, \quad \dot z = p\,R_p + q\,R_q, \quad \dot p = E_x, \quad \dot q = E_y .$$

Each solution traces a characteristic curve on the surface together with its orientation, a characteristic strip, and many closely spaced strips define the surface. The step in the image follows the gradient of the reflectance map, and the step in gradient space follows the brightness gradient of the image. For a rotationally symmetric reflectance map, as with a light at the viewer or in scanning electron microscope images, the characteristics are curves of steepest ascent or descent, and the strips grow outward from the singular points (Horn, 1990, Section 4.2; Horn and Brooks, 1989, Fig. 1).

The method needs an initial curve on which height and orientation are known, two boundary conditions at once in Durou, Falcone and Sagona’s reading (2008, Section 3.1.1). Brightness gradients measured from noisy images steer the strips off their ideal paths, errors accumulate along them, and in bad cases neighbouring strips cross. Horn’s remedy was to expand neighbouring strips in parallel and keep them consistent with a sharpening process, which also allowed new strips to be interpolated where existing ones spread apart (Horn, 1990, Section 4.1). Because the solution advances in one direction, information at the ends of the strips, such as the occluding boundary, cannot influence the result (Ikeuchi and Horn, 1981). Durou, Falcone and Sagona note that the method has been used mainly for theoretical uniqueness results.

Iterative methods on a grid: the variational approach

Iterative methods store an orientation at every grid point and repeatedly adjust it until the predicted shading matches the image, using local computations that can run in parallel; they are derived by minimizing a functional of the brightness error together with terms for smoothness or integrability. Ikeuchi and Horn (1981) introduced the approach, with stereographic coordinates that let the occluding boundary supply boundary conditions; Horn and Brooks (1986) pointed out that a smoothness penalty does not guarantee an integrable gradient field and derived a convergent scheme with a penalty on departure from integrability instead; and Horn (1990) coupled height and gradient in one scheme, linearized the reflectance map at each pixel and reduced the smoothness weight to zero near the solution, which recovers exact synthetic data to machine precision when boundary conditions for $p$ and $q$ are given. Variational shape from shading describes these solvers, later optimization methods, and an example in which Horn’s scheme and a plain least-squares solver meet the convex–concave ambiguity.

Local methods

Local methods compute the normal at each pixel independently of the others, from the first and second derivatives of brightness, by assuming a local shape. Pentland (1984) assumed the surface to be locally spherical and matched the brightness derivatives to those of an image of a sphere; Lee and Rosenfeld used the same assumption in a coordinate system aligned with the light, which makes the reflectance map rotationally symmetric and leads to estimates of slant and tilt from first derivatives only; Wildey’s radarclinometry assumed a locally cylindrical surface (all as Horn and Brooks, 1989, describe them). Durou, Falcone and Sagona (2008, Section 3.3.1) judge the assumption hard to justify: the normals obtained are usually far from integrable, and the shapes are poor unless the surface really is part of a sphere or a cylinder.

Linear methods

Linear methods approximate the reflectance map by a linear function of the gradient, which makes the image irradiance equation linear. For small gradients Horn (1990, Section 2.6) expands $R(p, q) \approx R(0, 0) + p\,R_p(0, 0) + q\,R_q(0, 0)$, which fails when the reflectance map is rotationally symmetric, since the first-order terms vanish, and which still leaves the ambiguity at right angles to the profiles. As Durou, Falcone and Sagona (2008, Section 3.3.2) describe the global linear methods, Pentland’s (1990) uses one linear approximation for the whole image and solves it with the Fourier transform, and later ones with finite differences; they note that a single approximation is reasonable only when the normal varies little over the surface, and that Pentland showed it improves as the angle between the viewing and lighting directions grows.

Tsai and Shah (1994) linearize locally rather than globally (Durou, Falcone and Sagona, 2008, Section 3.3.2). In the scheme as Durou, Falcone and Sagona (2008, Section 4.2.3) describe it, the backward differences $p_{ij} = u_{ij} - u_{i-1,j}$ and $q_{ij} = u_{ij} - u_{i,j-1}$ turn the image irradiance equation into one equation per pixel in the height, solved by a Newton-like update

$$u^{\text{new}}_{ij} = u^{\text{cur}}_{ij} - K_{ij}\left(r^{\text{cur}}_{ij} - I_{ij}\right),$$

where the gain $K_{ij}$ approximates the inverse of $\partial r/\partial p + \partial r/\partial q$ and is set to zero where that derivative vanishes. Durou, Falcone and Sagona point out that this is a variant of Newton–Raphson with no general proof of convergence; that the method imposes Neumann conditions at the boundary, so the height is determined only up to a constant and drifts from one iteration to the next; and that it is stopped after a fixed number of iterations (Tsai and Shah recommended two, Durou, Falcone and Sagona used five). They also found that even for a frontal light, the method worked only when the reflectance function was modelled with an oblique source, for which they had no explanation.

Viscosity solutions and propagation from the boundary

From the 1990s, methods for the partial differential equation used the theory of viscosity solutions, weak solutions that are typically Lipschitz continuous and so can represent surfaces with creases (Durou, Falcone and Sagona, 2008, Section 3.1.2). Starting with Rouy and Tourin, finite-difference schemes, a Markov chain approximation by Oliensis and Dupuis, fast-marching schemes for the eikonal equation, and semi-Lagrangian schemes such as Falcone and Sagona’s, which provably converge to the maximal solution, compute the solution on the grid from boundary data; with the state-constraint boundary condition the computed solutions grow inward from the boundary. Because the Dirichlet problem has several weak solutions, these methods either add information, such as the height at each singular point, or compute the maximal solution. Kimmel and Bruckstein instead propagate equal-height contours from the singular points with level sets (level set methods) and merge the local solutions.

A well-posed perspective model

Prados and Faugeras (2005) changed the image formation model rather than the solver. With a pinhole camera, a point light source at the optical centre (approximately the situation of an endoscope or a camera with a flash) and the $1/r^2$ fall-off of irradiance with distance kept in the brightness equation, they derived a new Hamilton–Jacobi equation and proved, for images that are differentiable and bounded away from zero, that it has a unique viscosity solution under “state constraints” on the image boundary, a condition that needs no boundary data and roughly requires the height to increase toward the image border, as for an object in front of a background. In this model the notion of singular points disappears, and with it the concave/convex ambiguity they cause. Their provably convergent scheme reconstructed a synthetic field of bumps with several local minima and a synthetic face without any additional data, gave what the authors call relatively satisfying results under strong noise, gamma, albedo and focal-length errors, and was applied to real face images after the eyes were painted over by hand, since the eyes violate the constant-albedo assumption.

Data-driven alternatives

Monocular depth estimation predicts depth from a single image with networks trained on data, without an explicit reflectance map or lighting model; the classical shape from shading literature surveyed here inverts a stated physical model instead.

Datasets and Benchmarks

The two comparative surveys rely mainly on synthetic images rendered from known surfaces, which, as Durou, Falcone and Sagona (2008) note, cannot be impossible images and can be compared with the true shape; real images with measured shape are few, and Durou, Falcone and Sagona add two.

Zhang, Tsai, Cryer and Shah (1999). As Durou, Falcone and Sagona (2008) describe this survey, it compared six methods: those of Lee and Kuo, Lee and Rosenfeld, Pentland, Tsai and Shah, Zheng and Chellappa, and Bichsel and Pentland, so that four of the six approximate the image irradiance equation; the test images included a synthetic vase, the Mozart image and the real Pepper, Lena and vase images; and Tsai and Shah’s method had the worst overall score (a total error of 59.3 against 41.3 for the best), although Durou, Falcone and Sagona judge it qualitatively the best of the six on the synthetic images.

Durou, Falcone and Sagona (2008). The survey compares one method per class, chosen so that they need the same data: Falcone and Sagona’s semi-Lagrangian scheme (FS), Daniel and Durou’s gradient descent on a functional with integrability and smoothness terms (DD), and Tsai and Shah’s method (TS). The three synthetic images are $256 \times 256$ renderings, under a frontal light, of a vase (SV, the same image as in Zhang et al.), a “Canadian tent” with sharp ridges (CT) and a smooth multi-peaked surface (DEM, Matlab’s peaks function); the real images are a vase (RV), whose shape is recovered by assuming rotational symmetry, the moulding of an elk’s head (Elk), measured at 1,104 points by mechanical metrology, and Pepper and Lena, which have no ground truth and violate the assumptions. Boundary conditions are part of the protocol: FS needs a Dirichlet condition, TS is run without one, and DD with either; shapes reconstructed without a Dirichlet condition are shifted by the constant that minimizes their RMS height error. Selected mean absolute errors (the survey’s $L^1$ errors) on the synthetic images follow; FS has the zero-height Dirichlet condition throughout, DD has none on SV and zero height on CT and DEM, and TS never has one:

Image	Method	Height error $\lvert u\rvert_1$	Normal error $\lvert n\rvert_1$	Grey-level error $\lvert I\rvert_1$
SV	FS	0.80	0.49	0.01
SV	DD	0.29	0.21	0.05
SV	TS	0.53	0.24	0.07
CT	FS	0.03	0.03	0.01
CT	DD	0.57	0.26	0.08
CT	TS	1.19	0.83	0.36
DEM	FS	0.52	0.49	0.01
DEM	DD	0.59	0.52	0.06
DEM	TS	0.64	0.67	0.12

Heights are in the survey’s units on the domain $[-6.4, 6.4]^2$; the normal error is the Euclidean norm of the difference of unit normals (bounded by 2) and the grey-level error is on images in $[0, 1]$ (Durou, Falcone and Sagona, 2008, Tables 1, 2 and 4). The zero-height condition is wrong for the vase, whose height does not vanish at the top and bottom of its outline, and with it FS produces a wrong shape, which the authors attribute to the concave/convex ambiguity; given the true boundary heights, FS and DD both improve (to 0.23 and 0.19 in height error). On the tent, FS follows the creases and is best, DD smooths them away and TS returns a nearly flat surface. On DEM, all three reconstructions are qualitatively wrong: FS returns the maximal solution, which is not the true surface, DD a minimum-energy configuration, and TS interprets every singular point as a local maximum. TS was the fastest in every test because it stops after five iterations. The authors conclude that FS and DD always converge (DD in general to a local minimum) while TS has no convergence guarantee, that DD and TS are accurate on smooth surfaces and FS adapts to smooth and non-smooth ones provided the surface is the maximal solution, and that FS always needs a Dirichlet condition, even a wrong one. The preprint announces a companion website with the codes and an online test interface; it was not checked for this article.

Evaluation Metrics

Horn (1990, Section 6.6) discusses which measures are informative on synthetic images. The best test compares the orientations of the computed surface with those of the true one, by the root-mean-square angle between normals or the RMS difference of gradients. Pointwise height comparisons are less useful, because errors at low spatial frequencies, which shading constrains least, dominate them and grow with the width of the region. Re-rendering the solution under the same light and comparing it with the input is not useful: the brightness error falls quickly in most iterative schemes, and a gradient field can reproduce the input image exactly without being integrable. To make the point, Horn built a gradient field by solving the photometric stereo equations for the input image and an arbitrary second image under two assumed lights: it reproduces the input under the original light, yet under the second light it shows the arbitrary image instead of the true surface. Rendering the solution and the true surface under a different light is a useful test, and measures of departure from integrability, especially the misfit $\iint (z_x - p)^2 + (z_y - q)^2$ between the recovered height and gradient, track convergence; they also serve as stopping criteria (Horn, 1990, Section 6.7).

Durou, Falcone and Sagona (2008, Section 5.2) report the mean absolute, root-mean-square and maximum errors of the height, of the unit normal and of the re-rendered grey level, computing normals and grey levels from the reconstructed height in the same way for every method.

Applications

• Planetary topography. The earliest quantitative use of shading, in the mid-1960s, recovered the shape of parts of the lunar surface in preparation for human exploration of the Moon, using the special reflectance of the lunar maria (Horn and Brooks, 1989). Photoclinometry in astrogeology has since concentrated on profiles across features such as ridges, grabens and craters (Horn, 1990, Section 3), and radarclinometry recovers slopes from synthetic aperture radar images, as in Wildey’s work for the Venus Radar Mapper (Horn and Brooks, 1989). Horn (1990) applied his algorithm to Viking Orbiter images of Mars, for which no ground truth was available.
• Satellite images of terrain. On a $108 \times 128$ patch of a SPOT image of a wooded hilly area near Huntsville, Alabama, with 10 m ground resolution, Horn (1990, Section 7.3) had to estimate the light direction and the atmospheric scatter from the image itself, and could not remove the smoothness term with free boundaries; the reconstruction had a little less than 200 m of relief where the actual relief is about 250 m. His conclusion was that quantitatively meaningful results need an accurate light source direction in camera coordinates, and the camera’s position and attitude if the result is to be related to the ground (remote sensing).
• Microscopy. Scanning electron microscope images have a rotationally symmetric reflectance map and are, as Horn and Brooks (1989) note, easy for people to interpret as shapes; Horn (1975) and Ikeuchi and Horn (1981) recovered shapes from them.
• Medical and close-range imaging. Models with the light at the centre of projection fit endoscopy and flash photography (Prados and Faugeras, 2005); Durou, Falcone and Sagona (2008) cite perspective methods applied to reconstruct organs from medical images and to digitize documents, and Prados and Faugeras reconstructed faces.
• Combination with other cues. Szeliski (2022, Section 13.1.1) notes that surfaces other than plaster casts rarely have a single uniform albedo, so shape from shading is usually combined with another technique, such as stereo matching or known texture: stereo and texture supply information in textured regions, while shading fills in uniformly coloured regions and adds finer detail. Horn and Brooks (1989) also cite a uniqueness result for the case where the image boundary corresponds to a known curve in space, of the kind stereo provides.

Open Problems

• Uniqueness and existence. Horn and Brooks (1989) called existence and uniqueness very difficult questions with only a few special cases settled; Durou, Falcone and Sagona (2008) still list open technical questions, such as uniqueness without continuity assumptions, and note that a global method under realistic assumptions is still lacking.
• Realistic image formation. Durou, Falcone and Sagona (2008) observe that the classical problem’s impact on applications has been rather limited, and Prados and Faugeras (2005) call the results of recent surveys quite unsatisfactory even on simple synthetic images that satisfy the model. Both point toward models closer to real images: perspective projection, nearby lights with fall-off, non-frontal lighting, non-Lambertian reflectance and interreflections.
• Lighting and albedo. Real scenes rarely have uniform albedo or a known light; with both unknown, the generalized bas-relief ambiguity applies (Belhumeur, Kriegman and Yuille, 1999).
• Optimization. Local minima and slow convergence on large images remain the practical obstacles of the minimization methods (Horn, 1990; Durou, Falcone and Sagona, 2008; Szeliski, 2022). Horn (1990) left open how much the implementation details that give exact solutions on synthetic data matter on real images, where errors from other sources may swamp them.