Photometric Stereo

Photometric stereo recovers the orientation of a surface at every pixel from several images taken by a fixed camera while the direction of the illumination changes between images. Because the camera does not move, every pixel sees the same surface point in every image, so there is no correspondence problem; what changes is the brightness, and the way the brightness of a pixel varies with the light direction reveals the surface normal there, along with the surface’s reflectance. The idea and the name come from Woodham (1980), who contrasted it with ordinary stereo vision: the method is photometric because it uses the values recorded at one image location in successive images, not the relative positions of displaced features. Szeliski (2022, Section 13.1.1) puts the analogy the other way round: the light sources play the role that cameras at different positions play in ordinary stereo. The output is a normal map, and, by integration, a depth map of the visible surface; Woodham’s conclusion summarized the complementarity: triangulation-based stereo measures distance and suits rough, textured surfaces with discontinuities, photometric stereo measures surface gradient and suits smooth surfaces with uniform reflectance.

Image formation describes the forward problem: how the reflectance of a surface, its orientation and the illumination combine into a pixel value. Photometric stereo inverts the shading part of that chain. Shape from shading attempts the same inversion from a single image under one light, which leaves one equation for two unknown orientation parameters per pixel and needs extra assumptions; photometric stereo adds images until the problem becomes a well-posed per-pixel estimation.

Problem Definition

The reflectance map. Woodham (1980) set the problem in terms of the reflectance map. Under orthographic projection, with the viewing direction along the negative $z$ axis, a surface $z = f(x, y)$ has gradient $(p, q) = (\partial f / \partial x, \partial f / \partial y)$, and its normal is proportional to $(p, q, -1)$ in Woodham’s coordinates (with $z$ measured toward the viewer, as in the example below, the normal is written $(-p, -q, 1)$). If a single distant source lights the scene and the reflected fraction of light depends only on the surface orientation, the image intensity is a function of the gradient alone, $I(x, y) = R(p, q)$. One image gives one equation in two unknowns per pixel, so it does not determine the orientation. Changing the light direction changes the reflectance map, so $k$ images taken under $k$ lights give $k$ equations $I_j(x, y) = R_j(p, q)$ for the same two unknowns.

The Lambertian case. For a perfectly diffuse, or Lambertian, surface, the reflected radiance is the same in all viewing directions and proportional to the cosine of the angle between the normal and the light (the matte term of the reflection models in image formation). With a unit normal $\mathbf{n}$, a light direction $\mathbf{l}_j$ whose length is the source’s strength, and an albedo (reflectance factor) $\rho$,

$$I_j = \rho \, \max\!\left(0,\; \mathbf{n}^\top \mathbf{l}_j\right).$$

Where the light is on the visible side of the surface, the equation is linear in the scaled normal $\mathbf{b} = \rho\,\mathbf{n}$. Woodham (1980) stacked three light directions as the rows of a matrix $N$ and the three intensities into a vector $\mathbf{I}$, so that $\mathbf{I} = \rho N \mathbf{n}$ and

$$\rho\,\mathbf{n} = N^{-1}\mathbf{I}, \qquad \rho = \lVert N^{-1}\mathbf{I} \rVert, \qquad \mathbf{n} = N^{-1}\mathbf{I} / \rho .$$

The inverse exists if and only if the three light directions do not lie in a plane. Three images therefore determine both the orientation and the albedo at each pixel, which no single image can. Woodham credits the three-source solution to a 1978 personal communication from Horn, and notes that for a reflectance function similar to that of the lunar maria, two images suffice if the albedo is known, provided the light directions are not collinear in azimuth. The paper itself is a revision of a 1978 SPIE seminar paper.

More than three lights. With $k > 3$ images, the stacked equations $\mathbf{I} = L\,\mathbf{b}$ (with $L$ the $k \times 3$ matrix of light vectors) are overdetermined and solved by least squares, $\hat{\mathbf{b}} = (L^\top L)^{-1} L^\top \mathbf{I}$, followed by $\rho = \lVert \hat{\mathbf{b}} \rVert$ and $\mathbf{n} = \hat{\mathbf{b}} / \rho$. Shi et al. (2016) call this classical method BASELINE in their benchmark, and Szeliski (2022, Section 13.1.1) states that the equations are well conditioned as long as three or more of the light vectors are linearly independent.

Assumptions. The classic formulation, as Shi et al. (2016) list it, assumes a Lambertian surface, a fixed orthographic camera with a linear radiometric response, and known directional lights. Each assumption defines a family of variants when it is relaxed.

Input and Output

• Input: $k \ge 3$ images from a fixed camera, each under a different light; for calibrated methods, the direction and strength of each light; usually a mask of the object pixels. The pixel values must be proportional to the light, which means raw or radiometrically linearized data rather than gamma-encoded images (image formation explains why). DiLiGenT, for example, provides per-light RGB intensities and asks users to divide each image by them before running photometric stereo.
• Output: a unit normal per pixel (a normal map) and, for Lambertian methods, an albedo map. Integrating the normals gives a depth map up to an additive constant; Szeliski (2022, Section 13.1.1) describes integration as a variant of regularized surface fitting, and names Nehab et al. (2005) and Harker and O’Leary (2008) for more sophisticated techniques.

Variants

Shi et al. (2016) organize the field by four properties, which they use as labels for every surveyed method: the reflectance model (Lambertian; Lambertian plus outliers handled by robust methods; analytic non-Lambertian models such as Ward or Cook-Torrance; general properties of the bidirectional reflectance distribution function, or BRDF, such as isotropy and monotonicity), whether the lighting is calibrated or uncalibrated, the lighting model, and the number of input images. In their notation Woodham’s method is LCDS: Lambertian, calibrated, directional lighting, small number of images.

• Calibrated and uncalibrated. Calibrated methods know the light directions, typically from a mirror sphere placed in the scene; Shi et al. note that this is tedious and that the sphere can cause interreflections on the object. Uncalibrated methods estimate the lighting as well, which brings the ambiguity described under Approaches.
• Lighting model. Directional lighting gives every pixel the same light direction. Point lighting models nearby sources whose direction and fall-off vary across the object, and has been used for handheld devices; general (environment) lighting treats the illumination as a spherical function, for example in spherical harmonics (Shi et al., 2016).
• Camera model. A perspective camera is more accurate than the orthographic model but makes even the Lambertian problem nonlinear; cameras with a nonlinear response either need radiometric calibration or methods that are insensitive to it (Shi et al., 2016).
• Number of images. Shi et al. group methods into small (at least three, typically 10–20 images), medium (typically 50–100) and large (typically 500–1,000) input sets; the more general the reflectance or lighting assumption, the more images the method tends to need.
• Multiple views and motion. Multi-view photometric stereo reconstructs a complete model rather than one normal map, and a coarse multi-view stereo shape can calibrate the lighting (Shi et al., 2016). Szeliski (2022, Section 13.1.1) describes the combination of multi-view stereo for coarse shape and photometric stereo for fine detail as an active research area.
• Colour-multiplexed lights. Woodham (1980) observed that the equivalent of photometric stereo can be obtained from a single view by separating several illuminations by colour.

Challenges

• Shadows. A pixel is in an attached shadow when its normal faces away from the light ($\mathbf{n}^\top\mathbf{l} \le 0$) and in a cast shadow when another part of the surface blocks the light. In both cases the observation is zero rather than the negative or positive value the linear model predicts, and fitting it as if it were linear biases the normal.
• Specular reflection. Glossy surfaces add a highlight around the mirror direction. Highlights are bright, depend on the viewing direction, and are not explained by the Lambertian model. Woodham (1980) already noted that a specular component crowds the reflectance-map contours into one region of gradient space at the expense of others.
• Interreflections. In concave regions light bounces from one part of the surface to another. Almost no photometric stereo method models this; DiLiGenT excludes the concave interior of one object from evaluation for that reason (Shi et al., 2016).
• Lighting geometry. Woodham (1980) showed how accuracy depends on the light configuration: with two lights 90° apart in azimuth, a measurement error of one grey level in ten constrains the orientation to within ±6.8° in one region of gradient space but only ±25.8° in another. A large angle between light and viewing direction makes the contours denser but puts more of gradient space in shadow, so the choice is a trade-off. He also pointed out that the sun’s path across the sky is very nearly planar, so outdoor images from a single day do not give three independent light directions.
• Calibration and linearity. Errors in light directions and intensities, nonlinear camera response and exposure changes enter the solution directly, which is why benchmark data are captured with linear cameras and calibrated lights.

Approaches

Lambertian least squares

Woodham’s method, with three lights or more, remains the reference point. It is closed-form, per-pixel, and exact when its assumptions hold. Every later method can be read as a way to deal with the observations that break the linear model $\mathbf{I} = L\mathbf{b}$: shadows, highlights, unknown $L$, or reflectance that is not Lambertian at all.

Uncalibrated photometric stereo and the bas-relief ambiguity

If the lights are unknown, the images still have a low-rank structure. Write the intensities of $P$ pixels under $Q$ lights as a $P \times Q$ matrix. Without shadows, it is the product of a $P \times 3$ matrix whose rows are the scaled normals $\mathbf{b}^\top$ and a $3 \times Q$ matrix whose columns are the light vectors, so its rank is at most three. Most uncalibrated methods are based on the factorization technique that Shi et al. (2016) attribute to Hayakawa (1994): a singular value decomposition gives a rank-3 factorization, but only up to an invertible $3 \times 3$ matrix $A$, because

$$\mathbf{b}^\top \mathbf{s} = (A\mathbf{b})^\top (A^{-\top}\mathbf{s})$$

for every normal $\mathbf{b}$ and light $\mathbf{s}$ (Belhumeur, Kriegman and Yuille, 1999, Eq. 5). According to Shi et al. (2016), Hayakawa resolved $A$ up to a rotation by identifying six pixels with the same albedo or six lights of the same intensity.

Belhumeur, Kriegman and Yuille (1999) showed how much of $A$ the images can determine. Not every transformed normal field corresponds to a surface: a field of normals comes from a surface $z = f(x, y)$ only if it is integrable, $f_{xy} = f_{yx}$. They proved that the only invertible linear transformations of a surface’s normal field that preserve integrability are those of the generalized bas-relief (GBR) group,

$$G = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \mu & \nu & \lambda \end{pmatrix}, \qquad \bar f(x, y) = \lambda f(x, y) + \mu x + \nu y, \qquad \bar{\mathbf{b}} = G^{-\top}\mathbf{b}, \qquad \bar{\mathbf{s}} = G\,\mathbf{s},$$

with a corresponding change of albedo (their Eq. 3). With $\mu = \nu = 0$ and $0 < \lambda < 1$ this is the flattening of a classical relief sculpture. The paper’s central result is that for a Lambertian surface viewed orthographically and lit by distant sources, every image of the surface is also an image of each of its GBR transforms under transformed lights, including the attached and cast shadows: the shadowed regions are the same because $\bar{\mathbf{b}}^\top\bar{\mathbf{s}} = \mathbf{b}^\top\mathbf{s}$ and because a GBR transformation maps a grazing light ray to a grazing light ray. For convex objects they also proved that the GBR transformation is the only transformation with this property. The consequences for photometric stereo, as the paper lists them:

• With unknown lights, at least three images and the integrability constraint determine the surface only up to a GBR transformation, a three-parameter family. Shadows give no further information.
• If the albedo is constant or known, or all lights have the same or known strength, the ambiguity shrinks to $\lambda = \pm 1$, $\mu = \nu = 0$: the surface is known up to a reflection in depth, the convex–concave ambiguity familiar from shape from shading. Cast shadows can distinguish the two.
• Known albedo alone does not recover the surface, contrary to Hayakawa (1994) as the paper reads it: it restricts $A$ to rotations and reflections, and a rotated normal field is not integrable, so it does not correspond to a surface.
• Specular reflection and interreflections depend on the relief, so they are not preserved by GBR transformations.

Later uncalibrated methods resolve the remaining three parameters with extra assumptions. Shi et al. (2016) summarize several: a natural image has a few dominant albedos, so the GBR transformation that minimizes the entropy of the albedo distribution is chosen (Alldrin et al., 2007, AM07); pixels with equal albedo are found by chromaticity clustering (Shi et al., 2010, SM10); local diffuse maxima, where the normal coincides with a light direction, give a closed-form solution (Papadhimitri and Favaro, 2014, PF14); or specular reflection, which is not GBR-invariant, constrains the transformation (Wu and Tan, 2013, WT13). Others use specular spikes, interreflections or a perspective camera.

Robust methods: shadows and highlights as outliers

The first family of non-Lambertian methods keeps the Lambertian model and discards the observations that violate it. Shi et al. (2016) trace it from early four-light methods that select the three most Lambertian measurements, through methods that select Lambertian subsets with graph cuts, RANSAC, maximum-likelihood estimation or the median, to methods that use robust statistics on the whole image matrix. As Shi et al. describe them, Wu et al. (2011, ACCV 2010) and Ikehata et al. (2012) model the observation matrix as a low-rank (Lambertian) matrix plus a sparse matrix of outliers from shadows and highlights; the first minimizes the rank, the second enforces rank three with sparse Bayesian regression. The advantage is that only the sparsity of the outliers is assumed, not their size or distribution. The limitation is that broad, soft specular lobes are not sparse, and such methods usually need more images for their statistics.

A simpler member of the family is the position threshold that Shi et al. (2016) evaluated: sort each pixel’s observations, discard the lowest fraction $T_\text{low}$ as shadows and the highest fraction above $T_\text{high}$ as highlights, and run least squares on the rest. On DiLiGenT they found that narrower intervals gave more accurate normals on most objects, and that the tightest setting they tried ($T_\text{low} = 40\%$, $T_\text{high} = 60\%$) was comparable to ST14 and IA14, the methods that performed best on most objects in their evaluation (below).

Analytic and general reflectance models

Instead of discarding highlights, the second family fits a reflectance model that explains them. Methods have used the Torrance–Sparrow model, the Ward model and mixtures of Ward models with spatially varying weights (Goldman et al., 2010, GC10), solving a nonlinear system for the material parameters together with the normals. Shi et al. (2016) note the drawbacks: each analytic model covers a narrow class of materials, and the fits are nonlinear and expensive.

The third family uses properties shared by broad classes of materials. Isotropic BRDFs give equal reflectance for normals symmetric about the plane spanned by the viewing and lighting directions, from which the azimuth of the normal can be computed; bivariate approximations of isotropic BRDFs, monotonicity of reflectance, and visibility constraints then give the elevation. Shi et al. (2016) list AZ08, ST12, HM10, ST14 and IA14 among these, and found ST14 and IA14 the most accurate of the calibrated methods they evaluated on most objects. Methods of this family often need far more images than the Lambertian baseline.

Learned methods

From 2017, networks trained on rendered data took over the role of the reflectance model. Santo et al. (2017), as Chen, Han and Wong (2018) and Ikehata (2018) describe their deep photometric stereo network (DPSN), concatenated a pixel’s observations under a fixed, predefined set of lights (Chen, Han and Wong give 96 as an example) into a vector and regressed the normal with a fully connected network, so the same light set was needed in training and testing.

PS-FCN (Chen, Han and Wong, 2018) removed that restriction. Each image is concatenated with its light direction, repeated as a three-channel image, and passed through a shared-weight convolutional feature extractor; the per-image feature maps are fused by max-pooling, which makes the network independent of the number and order of the inputs, and a regression network produces a normal map, trained with a cosine-similarity loss. The training images were rendered from two synthetic shape datasets (blobby and sculpture shapes) with the 100 measured BRDFs of the MERL database. Removing the light directions from the input gave an uncalibrated variant, UPS-FCN.

CNN-PS (Ikehata, 2018) works per pixel. It projects each light direction onto the image plane and writes the normalized intensity observed under that light into a $32 \times 32$ observation map, so that any number of unordered images becomes a fixed-size input for a small DenseNet-style convolutional neural network. For isotropic materials, rotating the lights and the normal together about the viewing axis leaves the observations unchanged, which Ikehata calls rotational pseudo-invariance and uses both to augment the training data and to average predictions over $K$ rotations at test time. The training set, CyclesPS, was rendered with the Cycles path tracer from fifteen 3D models with spatially varying parameters of the Disney principled BSDF, so that cast shadows and interreflections appear in the training data. Ikehata reports that CNN-PS handled the concave HARVEST object of DiLiGenT better than the other methods compared, and attributes its weaker result on READING to interreflected narrow highlights that were rare in the training data.

Example: Normals and Albedo from a Rendered Shape

The script renders a $128 \times 128$ orthographic height field with two bumps and a pit, textured with the scikit-image camera photograph as albedo (scaled to 0.25–1.0), under 24 known unit lights: 8 each at 15°, 35° and 55° from the viewing axis. The ground-truth normals are those of the height map, so the error is measured exactly, as the angle between estimated and true normals in degrees, the DiLiGenT measure. Five renders add one effect at a time: (1) the linear Lambertian model without clipping, (2) attached shadows, (3) cast shadows, found by marching from each pixel toward the light over the height map, (4) a Blinn–Phong highlight of strength 0.5 and exponent 50, and (5) Gaussian noise with standard deviation 0.01. Five estimators run on each: least squares on all images; least squares on observations above 0.02; Shi et al.’s position thresholds at (20%, 80%) and (40%, 60%), since the window is a free choice; and an $L_1$ regression by iteratively reweighted least squares. Part B factorizes the images without the lights, applies a GBR transformation, and checks integrability.

import numpy as np
from scipy import ndimage
from skimage import data, transform

H = W = 128
y, x = np.mgrid[0:H, 0:W].astype(float)
g = lambda cx, cy, h, s: h * np.exp(-((x - cx)**2 + (y - cy)**2) / (2 * s**2))
Z = g(45, 60, 30, 16) + g(95, 35, 14, 8) + g(90, 95, -20, 10)  # height (pixels)
alb = 0.25 + 0.75 * transform.resize(data.camera() / 255.0, (H, W),
                                     anti_aliasing=True)  # albedo map
inner = np.zeros((H, W), bool); inner[2:-2, 2:-2] = True

def normals(Z):                                  # n ~ (-z_x, -z_y, 1)
    zy, zx = np.gradient(Z)
    n = np.stack([-zx, -zy, np.ones_like(Z)], -1)
    return n / np.linalg.norm(n, axis=-1, keepdims=True)

L = np.array([[np.sin(p) * np.cos(a), np.sin(p) * np.sin(a), np.cos(p)]
              for k, p in enumerate(np.radians([15, 35, 55]))
              for a in np.radians(45 * np.arange(8) + 15 * k)])  # 24 lights

def cast_lit(Z, l, step=0.5):                    # march toward the light
    d = np.hypot(l[0], l[1]); lit = np.ones(Z.shape, bool); t = step
    while d > 0 and t < 1.5 * H:
        xs, ys = x + t * l[0] / d, y + t * l[1] / d
        inside = (xs >= 0) & (xs <= W - 1) & (ys >= 0) & (ys <= H - 1)
        zs = ndimage.map_coordinates(Z, [ys, xs], order=1, mode='nearest')
        lit &= ~(inside & (zs > Z + t * l[2] / d + 1e-9)); t += step
    return lit

def render(Z, alb, L, clip=True, cast=True, spec=0.0):
    n = normals(Z); out = []
    for l in L:
        s = np.linalg.norm(l); nl = n @ l           # |l| is the intensity
        im = alb * (np.maximum(nl, 0) if clip else nl)
        if spec:
            h = l / s + [0, 0, 1]; h /= np.linalg.norm(h)
            im = im + spec * s * np.maximum(n @ h, 0) ** 50 * (nl > 0)
        out.append(im * cast_lit(Z, l) if cast else im)
    return np.stack(out, -1).reshape(-1, len(L))

def wls(I, w):                                   # per-pixel weighted LS for b
    A = np.einsum('pq,qi,qj->pij', w, L, L)
    return np.linalg.solve(A, np.einsum('pq,pq,qi->pi', w, I, L)[..., None])[..., 0]

def methods(I):
    rank = np.argsort(np.argsort(I, 1), 1)
    out = {"LS all": wls(I, np.ones_like(I)),
           "LS I>0.02": wls(I, (I > 0.02) + 1e-12)}
    for lo, hi in [(20, 80), (40, 60)]:
        out[f"rank {lo}-{hi}%"] = wls(I, (rank >= lo * 0.24) & (rank < hi * 0.24))
    B = out["LS all"]
    for _ in range(30):                          # L1 by reweighting
        B = wls(I, 1 / np.maximum(np.abs(I - B @ L.T), 1e-4))
    out["L1 IRLS"] = B
    return out

N = normals(Z).reshape(-1, 3); m = inner.ravel()
def errors(B):
    rho = np.linalg.norm(B, axis=1)
    c = np.clip((B / rho[:, None] * N).sum(1), -1, 1)
    e = np.degrees(np.arccos(c))[m]
    return e.mean(), np.median(e), np.abs(rho / alb.ravel() - 1)[m].mean()

conds = {"1 linear, no shadows": render(Z, alb, L, clip=False, cast=False),
         "2 + attached shadows": render(Z, alb, L, cast=False),
         "3 + cast shadows": render(Z, alb, L),
         "4 + highlights": render(Z, alb, L, spec=0.5)}
rng = np.random.default_rng(0)
conds["5 + noise 0.01"] = conds["4 + highlights"] + rng.normal(0, 0.01, (H * W, 24))
print("mean / median angular error (deg), mean albedo relative error")
for name, I in conds.items():
    print(name)
    for k, B in methods(I).items():
        print("   %-12s %6.2f %6.2f  %.3f" % ((k,) + errors(B)))

# Part B: uncalibrated factorization and the generalized bas-relief ambiguity
for name in ["1 linear, no shadows", "2 + attached shadows", "4 + highlights"]:
    s = np.linalg.svd(conds[name], compute_uv=False)
    print(f"s4/s3 {name}: {s[3] / s[2]:.2e}")
U, s, Vt = np.linalg.svd(conds["1 linear, no shadows"], full_matrices=False)
Bs = U[:, :3] * np.sqrt(s[:3])                   # pseudo-normals, b up to A
print("pseudo-normal error: %.2f deg" % errors(Bs)[0])
A = np.linalg.lstsq(Bs, N * alb.reshape(-1, 1), rcond=None)[0]
print("after best linear map to truth: %.2e deg" % errors(Bs @ A)[0])

lam, mu, nu = 0.6, 0.4, -0.3
G = np.array([[1, 0, 0], [0, 1, 0], [mu, nu, lam]])
Zg = lam * Z + mu * x + nu * y                   # GBR-transformed surface
zy, zx = np.gradient(Z)
albg = alb / lam * np.sqrt(((lam * zx + mu)**2 + (lam * zy + nu)**2 + 1)
                           / (zx**2 + zy**2 + 1))   # Belhumeur et al., Eq. (3)
Lg = L @ G.T                                     # lights s -> G s
d0 = np.abs(render(Zg, albg, Lg) - conds["3 + cast shadows"]).max()
d1 = np.abs(render(Zg, albg, Lg, spec=0.5) - conds["4 + highlights"]).max()
c = np.clip((normals(Zg).reshape(-1, 3) * N).sum(1), -1, 1)
print(f"GBR images, Lambertian with shadows: max |dI| = {d0:.1e}")
print(f"GBR images, with highlights:         max |dI| = {d1:.2f}")
print(f"GBR normals differ from truth by {np.degrees(np.arccos(c))[m].mean():.1f} deg")

def curl(n):                                     # integrability residual
    p, q = -n[..., 0] / n[..., 2], -n[..., 1] / n[..., 2]
    return np.abs(np.gradient(p, axis=0) - np.gradient(q, axis=1))[inner].mean()
a = np.radians(20); R = np.array([[1, 0, 0], [0, np.cos(a), -np.sin(a)],
                                  [0, np.sin(a), np.cos(a)]])
print("integrability residual: true %.1e, GBR %.1e, rotated %.1e"
      % (curl(normals(Z)), curl(normals(Zg)), curl(normals(Z) @ R.T)))

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2; the script takes about 15 seconds, mostly for the shadow ray marching):

mean / median angular error (deg), mean albedo relative error
1 linear, no shadows
   LS all         0.00   0.00  0.000
   LS I>0.02      0.00   0.00  0.000
   rank 20-80%    0.00   0.00  0.000
   rank 40-60%    0.00   0.00  0.000
   L1 IRLS        0.00   0.00  0.000
2 + attached shadows
   LS all         0.49   0.00  0.006
   LS I>0.02      0.00   0.00  0.000
   rank 20-80%    0.00   0.00  0.000
   rank 40-60%    0.00   0.00  0.000
   L1 IRLS        0.00   0.00  0.000
3 + cast shadows
   LS all         1.06   0.00  0.008
   LS I>0.02      0.00   0.00  0.000
   rank 20-80%    0.08   0.00  0.000
   rank 40-60%    0.00   0.00  0.000
   L1 IRLS        0.00   0.00  0.000
4 + highlights
   LS all         4.70   2.88  0.215
   LS I>0.02      4.41   2.57  0.214
   rank 20-80%    3.78   3.45  0.131
   rank 40-60%    2.47   0.91  0.062
   L1 IRLS        3.09   1.48  0.094
5 + noise 0.01
   LS all         4.79   3.01  0.215
   LS I>0.02      4.54   2.67  0.215
   rank 20-80%    3.41   2.89  0.131
   rank 40-60%    3.72   2.68  0.078
   L1 IRLS        3.61   2.30  0.105
s4/s3 1 linear, no shadows: 4.73e-15
s4/s3 2 + attached shadows: 6.90e-02
s4/s3 4 + highlights: 5.80e-01
pseudo-normal error: 96.91 deg
after best linear map to truth: 2.45e-07 deg
GBR images, Lambertian with shadows: max |dI| = 1.4e-14
GBR images, with highlights:         max |dI| = 0.50
GBR normals differ from truth by 26.0 deg
integrability residual: true 1.4e-17, GBR 2.0e-17, rotated 6.8e-03

On the linear render every estimator recovers the normals and albedo exactly, which confirms that the solvers and the error measure are right before any effect is added. Shadows then bias plain least squares, because it fits zeros as if they were the negative values of the linear model. The bias is small on average, and the median of zero shows that it is concentrated on the few pixels that some light leaves in shadow. Dropping dark observations removes it exactly, since without noise a shadowed sample is exactly zero. The (20%, 80%) window drops the five lowest of 24 observations; in our reading, its small residual error with cast shadows comes from pixels with more than five shadowed observations, which keep some zeros.

Highlights affect every method, because they are bright, cover a neighbourhood of the mirror direction, and the threshold on dark values cannot see them. The narrow (40%, 60%) window discards most of them and gives the lowest error, in line with what Shi et al. (2016) found for position thresholds on DiLiGenT, and the $L_1$ fit helps but less: in our reading, a highlight lobe with exponent 50 touches several of the 24 observations at many pixels, so the outliers are not sparse enough for $L_1$ to ignore them. Albedo suffers most, since the extra light is absorbed into the length of $\mathbf{b}$. With noise the ordering changes: the narrow window keeps only five observations per pixel and averages the noise least, so its mean error rises above that of the (20%, 80%) window. Which window is best therefore depends on the balance between outliers and noise, which is why the threshold is printed both ways rather than chosen.

Part B shows the uncalibrated problem. The linear render has rank three to machine precision; shadows and especially highlights add significant further singular values, so the factorization’s rank-3 assumption is the first thing to fail on real data (Shi et al., 2016, name the Lambertian factorization as a major bottleneck). The rank-3 factors give pseudo-normals that are wrong by about 97° on average, yet a single $3 \times 3$ linear map turns them into the exact normals: the factorization recovers the normals up to $A$ and nothing more. The GBR test reproduces the paper’s central result: a surface flattened to 0.6 of its relief and tilted, with transformed albedo and lights, produces the same images to $10^{-14}$, attached and cast shadows included, although its normals differ by 26° on average. With highlights the images differ by up to 0.5, since specular reflection is not GBR-invariant. Finally, the GBR field is as integrable as the true one, while a 20° rotation of the normal field (which reproduces all Lambertian images with rotated lights, since $\mathbf{b}^\top\mathbf{s} = (R\mathbf{b})^\top(R\mathbf{s})$) violates integrability, which is why integrability reduces the general linear ambiguity to the GBR group.

Datasets and Benchmarks

Shi et al. (2016) found no earlier real dataset suitable for evaluating non-Lambertian photometric stereo: the public ones lacked ground-truth shapes, or had simple reflectance and few lighting conditions. DiLiGenT, named for its calibrated Directional Lightings, objects of General reflectance and “ground Truth” shapes, fills that gap for calibrated and uncalibrated methods under a fixed orthographic camera and directional lights. Both learned methods above report results on it, and Ikehata (2018) describes it as the largest real photometric stereo dataset. The protocol below follows the CVPR 2016 paper and the official dataset page as of November 2026; the TPAMI 2019 version, which the page asks users to cite, could not be read for this article.

• Objects. Ten objects, chosen to span reflectance and shape: mostly diffuse (CAT), rough (POT1), sparse specular spikes (BALL, READING), broad specular lobes on uniform (BEAR, BUDDHA) and spatially varying (POT2, GOBLET) materials, and metallic paint (COW, HARVEST); shapes range from a sphere to the concave HARVEST.
• Capture. A Point Grey Grasshopper camera with a linear response and a 50 mm lens recorded 2448 × 2048 images of objects about 20 cm across from about 1.5 m, to approximate orthographic projection, in a dark room with everything but the object covered in black cloth. Each light’s image is an HDR merge of four exposures. The 96 white LEDs sit on a 12 × 8 grid spanning 74.6° × 51.4°, with the camera at its centre.
• Calibration. Light intensities were calibrated with a Macbeth white balance chart; light directions by imaging a mirror sphere with all LEDs on and fitting the measured grid to the reflected rays with a rigid transformation. Point-light positions and camera intrinsics are also distributed for methods that use point lights or a perspective camera.
• Ground truth. Each shape was scanned with a Rexcan CS+ scanner (precision up to 0.01 mm) and stitched from 162 scans, then aligned to the images with a mutual-information method and checked by comparing the rendered normal map with one estimated by photometric stereo. The concave interior of GOBLET is excluded from evaluation because of strong interreflection.
• Released data. 16-bit PNG images at 612 × 512 for each of the 96 lights, an object mask, light directions and RGB intensities (images are divided by the intensities before use), and ground-truth normals for the main dataset. A test dataset images nine of the objects (all but BALL) from a different viewpoint and withholds their normals; as of November 2026 the page says that test results are computed manually by the authors from normal maps sent to them, and that an automatic evaluation site is still under construction.

Evaluation Metrics

The benchmark scores a normal map by the per-pixel angular error $\arccos(\mathbf{n}_0^\top\mathbf{n})$ in degrees between the ground-truth normal $\mathbf{n}_0$ and the estimate $\mathbf{n}$, inside the object mask (Shi et al., 2016). The mean angular error per object is the headline number in most papers; the benchmark paper also reports the minimum, maximum, median and quartiles as box plots. Papers then average the per-object means over the ten objects. Albedo and depth are not scored.

Selected mean angular errors on the DiLiGenT main dataset, averaged over the ten objects with all 96 images unless stated, are below. The averages for the non-learned methods are those of Chen, Han and Wong (2018, Table 2 for calibrated, Table 4 for uncalibrated methods), whose per-object values match the official summary page; DPSN, PS-FCN and UPS-FCN are from the same tables. Both PS-FCN rows use the model trained on both synthetic datasets with 32 inputs per sample, and the result with 16 test images is the average of 100 random trials. CNN-PS is from Ikehata (2018, Table 2), who reports one decimal place and discarded the first 20 images of BEAR as corrupted.

Method	Setting	Mean angular error (°)
BASELINE (Woodham least squares)	calibrated	15.39
WG10 (Wu et al., low rank)	calibrated	13.35
AZ08 (bivariate BRDF)	calibrated	12.61
IA14	calibrated	10.60
ST14	calibrated	10.30
DPSN	calibrated, learned	9.41
PS-FCN, 16 test images	calibrated, learned	9.37
PS-FCN, 96 test images	calibrated, learned	8.39
CNN-PS, $K = 10$	calibrated, learned	7.2
AM07	uncalibrated	37.25
PF14	uncalibrated	16.66
UPS-FCN	uncalibrated, learned	16.02

Per object the spread is large: on the official results the BASELINE error runs from 4.10° on BALL to 30.62° on HARVEST, and every calibrated method in the benchmark paper shows large errors in concave areas such as HARVEST, where cast shadows and interreflections are frequent. Uncalibrated methods are far less accurate than calibrated ones on non-Lambertian objects (Shi et al., 2016). The benchmark paper also fitted the best linear transformation (“Opt. A”) and the best GBR transformation (“Opt. G”) of the distorted normals that uncalibrated methods start from to the ground truth, and found surprisingly large residual errors; the official results give 3.37° on BALL and 27.14° on HARVEST for Opt. A, against the 2.45e-07° of the example above, where the images obey the Lambertian model exactly.

Numbers for the same method are comparable only when they come from the same run: the IW14 row on the official summary page and the IW14 row of Ikehata (2018, Table 2), which comes from the CNN-PS author’s own run of the method, differ (for BALL, 2.54° against 2.0°). The labels also need care: the CVPR 2016 paper evaluates IW12, the 2012 version of Ikehata et al.’s sparse-regression method, while the summary page has no IW12 row and labels the row for Ikehata et al.’s method IW14; Ikehata (2018) runs the 2014 journal version (TPAMI, sparse Bayesian regression).

Applications

Woodham (1980) proposed two uses: measuring surface orientation at every pixel, and the reverse problem of finding the image points with a given orientation, which he related to the industrial “bin-of-parts” problem of locating a known object for a robot to grasp; he described photometric stereo as practical wherever the illumination can be controlled, as in industrial inspection. Szeliski (2022, Section 13.1.1) describes its use for fine surface detail on top of multi-view stereo in high-quality laboratory scanning, a flash/no-flash variant with a stereo camera, and outdoor photometric stereo from webcam sequences that uses the sun’s trajectory as a moving light. In 3D reconstruction more broadly it belongs to the shape-from-X cues that need controlled lighting, and the normals it produces complement the depth from depth estimation methods, which measure distance rather than orientation.

Open Problems

Shi et al. (2016) closed their evaluation with open problems that their benchmark exposed: complicated BRDFs, spatially varying materials and concave shapes remain difficult; methods built on general reflectance properties need many more images than the Lambertian baseline and much more computation; global illumination degrades every evaluated method, and how to integrate its removal into photometric stereo is open; and uncalibrated methods are limited by the Lambertian factorization they start from, which outlier rejection fixes only when a dominant Lambertian component exists. Learned methods have since lowered the benchmark errors, but they inherit the reflectance and light transport of their synthetic training data, as Ikehata’s analysis of READING shows. Beyond the classic setting, the methods Shi et al. (2016) survey for near point lights, perspective cameras, general environment lighting, multiple views and moving objects each relax one of the assumptions of the fixed-camera, distant-light formulation.