Variational Shape from Shading

Variational shape from shading solves the shape from shading problem on a regular grid: it minimizes a functional made of the brightness error and of terms that express smoothness or integrability, and computes the minimum with iterative schemes derived from the functional’s Euler equations or by direct minimization. The approach was introduced by Ikeuchi and Horn (1981) and developed by Horn and Brooks (1986) and Horn (1990); Durou, Falcone and Sagona (2008) classify such methods as methods using optimization.

Problem

Given an image $E(x, y)$ and a reflectance map $R(p, q)$, find the surface whose gradient $(p, q) = (z_x, z_y)$ satisfies the image irradiance equation $E(x, y) = R(p, q)$. One brightness value per pixel constrains two gradient components, the gradient field must be integrable, $p_y = q_x$, and the solution can be ambiguous; shape from shading sets out the equation, its ambiguities and the boundary data that make it determinate. A variational method replaces the equation by the minimization, over a grid, of the brightness error

$$\iint \big(E(x, y) - R(p, q)\big)^2\, dx\, dy$$

plus penalty terms, which select one of the many gradient fields that fit the image (Horn, 1990, Section 4.4).

Inputs and Outputs

• Input: the image, the reflectance map, and boundary data: orientations on the occluding boundary in Ikeuchi and Horn’s (1981) scheme; in Horn’s (1990), heights and gradients on the border of the image region where they are known and natural boundary conditions where they are not; and initial values.
• Output: a field of surface orientations on the grid, from which heights follow by a least-squares fit; or, in Horn’s (1990) scheme, heights and gradients together.

Intuition

The search for alternatives, according to Horn and Brooks (1989), began because characteristic strips suggested neither efficient, robust computation nor plausible biological mechanisms. 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. Strat (1979) minimized a weighted sum of brightness errors and of the errors of height integrals around small loops, which vanish for an integrable field; his method needed the orientation on a closed curve (Ikeuchi and Horn, 1981).

Algorithm

Smoothness and stereographic coordinates (Ikeuchi and Horn, 1981)

Ikeuchi and Horn (1981) were the first to derive iterative shape from shading schemes with the calculus of variations, inspired, as Horn and Brooks (1989) put it, by the variational estimation of optical flow in the Horn–Schunck method. To use the occluding boundary, they parametrized orientation by stereographic coordinates

$$f = \frac{2p}{p^2 + q^2}\left(\sqrt{1 + p^2 + q^2} - 1\right), \qquad g = \frac{2q}{p^2 + q^2}\left(\sqrt{1 + p^2 + q^2} - 1\right),$$

a projection of the Gaussian sphere from its south pole that maps the equator, the orientations on the occluding boundary, to a circle of radius two instead of to infinity. They minimized the sum over the grid of a departure-from-smoothness term (the squared differences of $f$ and $g$ between neighbours) and $\lambda$ times the squared brightness error, which gives the iteration

$$f_{ij}^{n+1} = f_{ij}^{*n} + \lambda\left[E_{ij} - R_s\!\left(f_{ij}^{*n}, g_{ij}^{*n}\right)\right]\frac{\partial R_s}{\partial f},$$

and the same for $g$, where $f^{*}$ and $g^{*}$ are averages over the four neighbours. Evaluating the reflectance map at the local average rather than at the cell’s own value damps a checkerboard instability. The orientations on the occluding boundary supply the boundary conditions, and the values at singular points supply initial values. On a synthetic Lambertian sphere 30 grid nodes across, lit from near the viewer, the relative root-mean-square orientation error fell below 0.01% after 30 iterations, about the number needed for information to cross the object. They also applied the method to scanning electron microscope images.

Integrability (Horn and Brooks, 1986)

Horn and Brooks (1986) unified these schemes as minimizations of functionals and pointed out the cost of the smoothness term. Minimizing the brightness error alone is ill-posed in terms of orientation fields, since infinitely many fields fit; regularization picks the smoothest, but it does not guarantee integrability. Enforcing integrability strictly led to no convergent iterative scheme they could find; a penalty on departure from integrability, $\iint (E - R)^2 + \lambda\,(p_y - q_x)^2$, did. They derived further schemes in unit normals that can use occluding boundaries. Horn (1990, Section 4.4) adds that a solution of the regularized problem is not a solution of the original one: an iteration with a smoothness penalty moves away from the correct surface even when started on it, because it can trade a little brightness error for more smoothness, and its results are too smooth, so it does well only on simple smooth shapes.

Coupled height and gradient (Horn, 1990)

Horn (1990) combined height and gradient in one scheme. It minimizes

$$\iint \big(E - R(p, q)\big)^2 + \lambda\left(p_x^2 + p_y^2 + q_x^2 + q_y^2\right) + \mu\left((z_x - p)^2 + (z_y - q)^2\right) dx\, dy,$$

whose Euler equations are three coupled second-order equations for $p$, $q$ and $z$, satisfied exactly by a true solution with $p = z_x$, $q = z_y$ and $E = R$. Three ingredients make it work. The reflectance map is linearized about the current estimate of the gradient at each pixel, which Horn found speeds convergence and lets the scheme converge where simpler ones diverge or stick in local minima. The smoothness weight $\lambda$ is needed only to stabilize the iteration far from the solution and is reduced to zero as the solution is approached, so it does not distort the result. And the discrete estimators are made consistent: heights live on a grid offset by half a pixel from the gradients and the image, the first-derivative estimators average two differences over a $2 \times 2$ block, and the Laplacian estimator matches their repeated application, which leaves a checkerboard pattern of heights invisible to the derivative estimators that boundary conditions must handle with care. With all details attended to and boundary conditions for $p$ and $q$ given, the algorithm finds the solution of exact synthetic data to machine precision. On $64 \times 64$ synthetic images of Gaussians, sums of Gaussian blobs and sinusoidal gratings under linear, oblique Lambertian and rotationally symmetric reflectance maps, the normals were typically within a degree or two after a few hundred iterations; on a $231 \times 178$ image rendered from a steep, gullied digital terrain model, the simpler scheme without linearization often stuck in a local minimum after a few iterations, while the linearized one reached the solution to machine precision after a few thousand.

Complexity

In Horn’s (1990) experiments the number of iterations grew almost quadratically with image size, so total computation grew almost with its fourth power, which Horn proposed to address with multigrid methods.

Implementation

The script renders a $64 \times 64$ Lambertian image (albedo 1) of a surface with a raised-cosine bump (height 8 pixels, radius 14) and a dimple of the same size on a flat plane, so the border of the image has height zero. Heights live on the $65 \times 65$ corners of the pixels and the gradient at each pixel comes from the $2 \times 2$ block of heights around it (Horn’s staggered grid and Eq. 94), with the conventions of shape from shading: $z$ toward the viewer, $\mathbf{n} \propto (-p, -q, 1)$, $x$ along columns, $y$ along rows, and the light $\mathbf{l}$ a unit vector toward the source. The oblique light is tilted 30° from the viewing axis, at an azimuth of 30° from $+x$; the maximum slope is 0.89 and $\mathbf{n}^\top\mathbf{l} \ge 0.31$ everywhere, so there are no shadows. The array row index increases downward, so the light’s azimuth turns from $+x$ toward the lower part of the image. Rendering with the solvers’ own discrete gradient guarantees an exact discrete solution (Horn, 1990, Section 5.1).

Two solvers run from the same starts. The first minimizes the brightness error alone with the heights as unknowns, by Levenberg–Marquardt with a sparse Jacobian (the mechanics are those of nonlinear least squares; at most 300 iterations). The second is Horn’s (1990) scheme as described above (Eqs. 75 and 84–90, natural boundary conditions for $p$ and $q$), with the over-relaxation factor 1.7 and the 15-stage schedule printed in his Fig. 5 for a $64 \times 64$ image of a crater: the smoothness weight falls from 1 to 0 and $\mu$ from 0.1 to 0.001 over 12,288 iterations. Four details are our choices where the paper leaves them open: we read the Fig. 5 fields “Cor-pq 0.25” and “Cor-z 1.0” as the corner weights of the averaging stencils, read two faint $\mu$ values in the figure as 0.08 and 0.008, order the over-relaxation by four colours of pixels so that no cell is updated together with a neighbour it uses (Horn’s Section 6.4 describes such colourings), and impose the natural boundary condition by copying the border values of $p$ and $q$ outward for the local averages, whereas Horn’s Eqs. 99–100 overwrite the border cells from the interior. Cases:

• A. Oblique light, border heights given (zero), flat start.
• B. Oblique light, free border: every height is unknown except one corner, which fixes the additive constant. Only Levenberg–Marquardt runs here; Horn’s scheme needs his separate boundary scheme for a free height border, and he reports that it is not stable without the smoothness term when the boundary is completely free (Horn, 1990, Section 7.2).
• C. Frontal light, border heights given. The four surfaces obtained by flipping the bump, the dimple or both are rendered; the Jacobian of the brightness error at the flat start is printed; and both solvers start from a flat surface with Gaussian noise of standard deviation 0.01 pixels for seeds 0–4. Each result is compared with the nearest of the four sign patterns, written as multipliers of the true bump and dimple.
• D. Oblique light, but the solver is told an azimuth 180° away from the true one; flat start.
• E. Oblique light, starting from the inverted surface $-z$.

Measures: the RMS brightness error of the result under the light given to the solver, the RMS height error over all corners, and the mean angle between the result’s normals and the true ones.

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve

N = 64                                   # 64 x 64 pixels; heights on the 65 x 65 corners
Y, X = np.mgrid[0:N + 1, 0:N + 1].astype(float)

def bump(cx, cy, h, R):                  # raised cosine, exactly zero outside radius R
    r = np.hypot(X - cx, Y - cy)
    return np.where(r < R, h * np.cos(np.pi * r / (2 * R)) ** 2, 0.0)

B1, B2 = bump(20, 22, 8.0, 14), bump(44, 42, -8.0, 14)
Z = B1 + B2                              # a bump and a dimple on a plane, border = 0

def grad(z):                             # Horn (1990) eq. 94: p, q at pixel centres
    p = ((z[:-1, 1:] - z[:-1, :-1]) + (z[1:, 1:] - z[1:, :-1])) / 2
    q = ((z[1:, :-1] - z[:-1, :-1]) + (z[1:, 1:] - z[:-1, 1:])) / 2
    return p, q

def refl(p, q, l):                       # Lambertian R and its derivatives, n ~ (-p, -q, 1)
    s = np.sqrt(1 + p**2 + q**2); num = -p * l[0] - q * l[1] + l[2]
    return num / s, -l[0] / s - num * p / s**3, -l[1] / s - num * q / s**3

def shade(z, l):
    return np.maximum(refl(*grad(z), l)[0], 0)

t = np.radians(30)
l_obl = np.array([np.sin(t) * np.cos(t), np.sin(t) * np.sin(t), np.cos(t)])
l_front = np.array([0.0, 0.0, 1.0])
l_flip = np.array([-l_obl[0], -l_obl[1], l_obl[2]])     # azimuth off by 180 degrees

# ---- Solver 1: Levenberg-Marquardt on the brightness error, heights as unknowns ----
k, lq = np.mgrid[0:N, 0:N]
res_idx = (k * N + lq).ravel()
dp = {(0, 0): -0.5, (0, 1): 0.5, (1, 0): -0.5, (1, 1): 0.5}
dq = {(0, 0): -0.5, (0, 1): -0.5, (1, 0): 0.5, (1, 1): 0.5}

def residual_jac(z, E, l, free):
    p, q = grad(z)
    R, Rp, Rq = refl(p, q, l)
    rows, cols, vals = [], [], []
    for (a, b) in dp:
        node = ((k + a) * (N + 1) + (lq + b)).ravel()
        v = (Rp * dp[(a, b)] + Rq * dq[(a, b)]).ravel()
        keep = free[node]
        rows.append(res_idx[keep]); cols.append(node[keep]); vals.append(v[keep])
    J = sp.csr_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
                      shape=(N * N, (N + 1) ** 2))
    return (R - E).ravel(), J[:, np.flatnonzero(free)]

def solve_lm(E, l, z0, free, max_iter=300):
    z = z0.copy().ravel(); idx = np.flatnonzero(free); mu = 1e-3
    r, J = residual_jac(z.reshape(N + 1, N + 1), E, l, free)
    cost = r @ r
    for it in range(1, max_iter + 1):
        A = (J.T @ J + mu * sp.identity(J.shape[1])).tocsc()
        step = spsolve(A, -(J.T @ r))
        zt = z.copy(); zt[idx] += step
        rt, Jt = residual_jac(zt.reshape(N + 1, N + 1), E, l, free)
        if rt @ rt < cost:
            z, r, J, cost = zt, rt, Jt, rt @ rt; mu /= 3
        else:
            mu *= 4
        if np.sqrt(cost / r.size) < 1e-13 or np.linalg.norm(step) < 1e-14:
            break
    return z.reshape(N + 1, N + 1)

# ---- Solver 2: Horn (1990), coupled height and gradient, Fig. 5 schedule ----
A_PQ, A_Z, ALPHA = 0.25, 1.0, 1.7       # stencil corner weights and SOR factor (Fig. 5)
K_PQ, K_Z = 4 / (A_PQ + 1), 4 / (A_Z + 1)
SCHEDULE = [(1.0, 0.1, 8), (0.5, 0.1, 32), (0.2, 0.1, 64), (0.1, 0.1, 128), (0.05, 0.1, 256),
            (0.02, 0.08, 512), (0.01, 0.05, 1024), (0.005, 0.03, 2048), (0.002, 0.02, 3072),
            (0.001, 0.01, 4096), (0.0005, 0.008, 5120), (0.0003, 0.007, 6144),
            (0.0002, 0.005, 8192), (0.0001, 0.002, 10240), (0.0, 0.001, 12288)]
ki, li = np.mgrid[0:N, 0:N]
COL_PQ = [(ki % 2 == a) & (li % 2 == b) for a in (0, 1) for b in (0, 1)]
kz, lz = np.mgrid[1:N, 1:N]
COL_Z = [(kz % 2 == a) & (lz % 2 == b) for a in (0, 1) for b in (0, 1)]

def avg(f, a):                           # local average of the eq. 45 stencil, edge copy
    g = np.pad(f, 1, mode="edge")
    edge = (g[:-2, 1:-1] + g[2:, 1:-1] + g[1:-1, :-2] + g[1:-1, 2:]) / 4
    corn = (g[:-2, :-2] + g[:-2, 2:] + g[2:, :-2] + g[2:, 2:]) / 4
    return a * corn + (1 - a) * edge

def pq_update(p, q, z, E, l, lam, mu):   # eqs. 86-90, linearized about the old p, q
    R, Rp, Rq = refl(p, q, l)
    zx, zy = grad(z)
    l2 = K_PQ * lam + mu
    A = K_PQ * lam * (avg(p, A_PQ) - p) + mu * (zx - p) + (E - R) * Rp
    B = K_PQ * lam * (avg(q, A_PQ) - q) + mu * (zy - q) + (E - R) * Rq
    D = l2 * (l2 + Rp**2 + Rq**2)
    return p + ((l2 + Rq**2) * A - Rp * Rq * B) / D, q + ((l2 + Rp**2) * B - Rq * Rp * A) / D

def z_update(z, p, q):                   # eq. 75, diagonal stencil, interior nodes
    zbar = (z[:-2, :-2] + z[:-2, 2:] + z[2:, :-2] + z[2:, 2:]) / 4
    px = ((p[:-1, 1:] - p[:-1, :-1]) + (p[1:, 1:] - p[1:, :-1])) / 2
    qy = ((q[1:, :-1] - q[:-1, :-1]) + (q[1:, 1:] - q[:-1, 1:])) / 2
    return zbar - (px + qy) / K_Z

def solve_horn(E, l, z0):                # border heights stay fixed (Dirichlet)
    z = z0.copy(); p, q = grad(z); it = 0
    for lam, mu, end in SCHEDULE:
        while it < end:
            for c in COL_PQ:             # SOR over four colours
                pn, qn = pq_update(p, q, z, E, l, lam, mu)
                p = np.where(c, p + ALPHA * (pn - p), p); q = np.where(c, q + ALPHA * (qn - q), q)
            for c in COL_Z:
                inner = z[1:-1, 1:-1]
                z[1:-1, 1:-1] = np.where(c, inner + ALPHA * (z_update(z, p, q) - inner), inner)
            it += 1
    return z

# ---- Measures ----
def normals(z):
    p, q = grad(z); n = np.stack([-p, -q, np.ones_like(p)], -1)
    return n / np.linalg.norm(n, axis=-1, keepdims=True)

def report(name, z, E, l, ref=Z):
    ang = np.degrees(np.arccos(np.clip((normals(z) * normals(ref)).sum(-1), -1, 1))).mean()
    print("  %-26s brightness RMS %.1e  depth RMSE %.2e px  normal error %.2e deg"
          % (name, np.sqrt(((shade(z, l) - E) ** 2).mean()), np.sqrt(((z - ref) ** 2).mean()), ang))

border = np.ones((N + 1, N + 1), bool)
border[0, :] = border[-1, :] = border[:, 0] = border[:, -1] = False
flat = np.zeros_like(Z)
p, q = grad(Z)
print("max slope %.2f, min n.l under oblique light %.2f" % (np.hypot(p, q).max(), refl(p, q, l_obl)[0].min()))

E_obl = shade(Z, l_obl)
print("A. oblique light, border heights given, flat start")
report("LM", solve_lm(E_obl, l_obl, flat, border.ravel()), E_obl, l_obl)
report("Horn (1990)", solve_horn(E_obl, l_obl, flat), E_obl, l_obl)

print("B. oblique light, free border (one corner fixed), flat start")
corner = np.ones((N + 1, N + 1), bool); corner[0, 0] = False
report("LM", solve_lm(E_obl, l_obl, flat, corner.ravel()), E_obl, l_obl)

print("C. frontal light, border heights given")
E_fr = shade(Z, l_front)
for s1, s2 in [(1, 1), (-1, 1), (1, -1), (-1, -1)]:
    print("  bump x %+d, dimple x %+d: max |image difference| %.1e"
          % (s1, s2, np.abs(shade(s1 * B1 + s2 * B2, l_front) - E_fr).max()))
r0, J0 = residual_jac(flat, E_fr, l_front, border.ravel())
print("  max |Jacobian| at the flat start: %.1e" % np.abs(J0.data).max())
combos = {(a, b): a * B1 + b * B2 for a in (1, -1) for b in (1, -1)}
for seed in range(5):
    z0 = flat + np.random.default_rng(seed).normal(0, 0.01, Z.shape) * border
    for name, solver in [("LM", solve_lm), ("Horn", solve_horn)]:
        z = solver(E_fr, l_front, z0, border.ravel()) if name == "LM" else solver(E_fr, l_front, z0)
        best = min(combos, key=lambda c: np.sqrt(((z - combos[c]) ** 2).mean()))
        print("  seed %d %-4s brightness RMS %.1e  nearest pattern (bump x %+d, dimple x %+d), depth RMSE %.2e px"
              % (seed, name, np.sqrt(((shade(z, l_front) - E_fr) ** 2).mean()), *best,
                 np.sqrt(((z - combos[best]) ** 2).mean())))

print("D. light azimuth assumed 180 degrees off, flat start")
for name, z in [("LM", solve_lm(E_obl, l_flip, flat, border.ravel())), ("Horn (1990)", solve_horn(E_obl, l_flip, flat))]:
    report(name + " vs truth", z, E_obl, l_flip)
    report(name + " vs inverted truth", z, E_obl, l_flip, ref=-Z)

print("E. oblique light, start from the inverted surface")
report("LM", solve_lm(E_obl, l_obl, -Z, border.ravel()), E_obl, l_obl)
report("Horn (1990)", solve_horn(E_obl, l_obl, -Z.copy()), E_obl, l_obl)

Output (NumPy 2.3.5, SciPy 1.16.3; the script takes about three minutes and gives identical output on repeated runs):

max slope 0.89, min n.l under oblique light 0.31
A. oblique light, border heights given, flat start
  LM                         brightness RMS 4.4e-03  depth RMSE 9.53e-01 px  normal error 6.71e+00 deg
  Horn (1990)                brightness RMS 4.7e-06  depth RMSE 1.34e-03 px  normal error 2.26e-02 deg
B. oblique light, free border (one corner fixed), flat start
  LM                         brightness RMS 5.1e-04  depth RMSE 1.12e+00 px  normal error 7.45e+00 deg
C. frontal light, border heights given
  bump x +1, dimple x +1: max |image difference| 0.0e+00
  bump x -1, dimple x +1: max |image difference| 0.0e+00
  bump x +1, dimple x -1: max |image difference| 0.0e+00
  bump x -1, dimple x -1: max |image difference| 0.0e+00
  max |Jacobian| at the flat start: 0.0e+00
  seed 0 LM   brightness RMS 4.1e-04  nearest pattern (bump x -1, dimple x -1), depth RMSE 1.78e+00 px
  seed 0 Horn brightness RMS 3.1e-06  nearest pattern (bump x -1, dimple x +1), depth RMSE 2.15e-02 px
  seed 1 LM   brightness RMS 3.7e-04  nearest pattern (bump x +1, dimple x -1), depth RMSE 1.80e+00 px
  seed 1 Horn brightness RMS 3.0e-06  nearest pattern (bump x +1, dimple x -1), depth RMSE 2.10e-02 px
  seed 2 LM   brightness RMS 3.5e-04  nearest pattern (bump x -1, dimple x +1), depth RMSE 1.76e+00 px
  seed 2 Horn brightness RMS 3.1e-06  nearest pattern (bump x -1, dimple x +1), depth RMSE 2.17e-02 px
  seed 3 LM   brightness RMS 2.5e-04  nearest pattern (bump x -1, dimple x +1), depth RMSE 1.79e+00 px
  seed 3 Horn brightness RMS 2.0e-06  nearest pattern (bump x +1, dimple x +1), depth RMSE 7.39e-03 px
  seed 4 LM   brightness RMS 3.3e-04  nearest pattern (bump x +1, dimple x +1), depth RMSE 1.80e+00 px
  seed 4 Horn brightness RMS 2.0e-06  nearest pattern (bump x -1, dimple x -1), depth RMSE 7.54e-03 px
D. light azimuth assumed 180 degrees off, flat start
  LM vs truth                brightness RMS 4.4e-03  depth RMSE 2.90e+00 px  normal error 1.77e+01 deg
  LM vs inverted truth       brightness RMS 4.4e-03  depth RMSE 9.53e-01 px  normal error 6.71e+00 deg
  Horn (1990) vs truth       brightness RMS 4.7e-06  depth RMSE 3.59e+00 px  normal error 1.71e+01 deg
  Horn (1990) vs inverted truth brightness RMS 4.7e-06  depth RMSE 1.34e-03 px  normal error 2.26e-02 deg
E. oblique light, start from the inverted surface
  LM                         brightness RMS 2.7e-02  depth RMSE 2.06e+00 px  normal error 1.04e+01 deg
  Horn (1990)                brightness RMS 4.7e-06  depth RMSE 1.34e-03 px  normal error 2.26e-02 deg

Case A reproduces the documented limit of the simple formulation. Levenberg–Marquardt, which linearizes the whole brightness error at each step and solves for all heights at once, still ends at a surface whose brightness error is far from zero and whose relief is wrong by about a pixel; in our reading this is the trap of local minima that Horn (1990, Sections 5.1 and 7.4), Durou, Falcone and Sagona (2008, Section 3.2.3) and Szeliski (2022, Section 13.1.1) describe. Horn’s scheme, with the same image, boundary and start, recovers the surface to within Horn’s own description of his results, though not to machine precision. Horn (1990, Section 7) reports that the brightness error falls rapidly while integrability improves much more slowly and that the shape becomes exact only eventually, and his machine-precision result assumes boundary conditions for $p$ and $q$, which this example leaves free. Since the smoothness weight is zero in the last stage, it does not bias the final result.

Case B removes the boundary data. Levenberg–Marquardt ends with a lower brightness error than in case A but a larger height error. The problem has 128 more unknowns than pixel equations (4,224 against 4,096; our count), so many surfaces fit the image exactly, the true one among them, yet the solver still stops short of a zero brightness error.

Case C is the convex–concave ambiguity in its purest form. Under a frontal light all four sign patterns render to exactly the same image, and all four satisfy the same zero-height border, so the boundary condition does not decide between them, as Durou, Falcone and Sagona (2008) state. The Jacobian at the flat start is exactly zero, because a flat surface facing the light is a singular point everywhere, so a solver started there cannot move at all; hence the small random starts. Horn’s scheme ends within about 0.02 pixels of one of the four exact solutions for every seed, and which one depends only on the noise in the start: across the five seeds, every one of the four patterns appears. Levenberg–Marquardt stalls for every seed, more than a pixel away from all four.

Case D is the crater illusion. With the azimuth wrong by 180°, Horn’s scheme returns the inverted surface and Levenberg–Marquardt the inverse of its case A result, each with exactly its case A errors: inverting the relief and the light’s azimuth together is an exact symmetry of the discrete problem, and the flat start is unchanged by it, so each computation is the mirror image of case A (our reasoning). The image cannot tell the two interpretations apart; only the assumed light does.

Case E starts on the wrong side of the ambiguity under the correct, oblique light. Here the inverted surface is not a solution, and Horn’s scheme returns to the true surface with the same printed accuracy as from the flat start, while Levenberg–Marquardt stays far from it with a large brightness error.

The example was set up after a pilot in which only the Levenberg–Marquardt solver had been planned; it did not recover the surface in case A, so Horn’s scheme was added. A second pilot of Horn’s scheme on case A recovered the surface but not to the planned “depth RMSE and angular error near machine precision”, so the success criterion was set to Horn’s own description of his results (normals within a degree or two after a few hundred iterations, the exact shape eventually); Horn’s printed schedule was kept unchanged. After the first run, one of two printed fields in case C was removed and the other relabelled, because the removed field gave the sign of the mean height where the other gave multipliers of the true features; no computed value changed.

Variants

Later optimization methods, as Durou, Falcone and Sagona (2008, Section 3.2) survey them, differ in three choices. The unknowns may be the height itself, rarely used because of reported problems of convergence or slowness; the gradient $(p, q)$; height and gradient together; stereographic coordinates; the normal; or a parametric surface model. The functional combines the brightness error with an integrability term, a smoothness term or both, or imposes integrability as a hard constraint at each step, as in Frankot and Chellappa’s method. The minimization either solves the Euler equations, which is easier and faster but can diverge, or minimizes directly, for example by gradient descent with line search. The solutions are typically local minima; global minima need stochastic algorithms such as simulated annealing at a much higher cost, and multiresolution reduces the time. Szeliski (2022, Section 13.1.1) summarizes the same experience: shape from shading is susceptible to local minima and can converge slowly, multi-resolution techniques help with speed and more sophisticated optimization with local minima.