Recover Normal from Three Lights

Photometric stereo with three lights is about solving a small linear system at each pixel to recover the surface normal from measured intensities under known lighting directions.

Under a Lambertian model, the image intensity from a directional light is proportional to the dot product between the light direction and the surface normal, scaled by the albedo (reflectance). With three lights, you get three linear equations in the three unknown components of the vector $\rho \mathbf{n}$. If the three light directions are linearly independent (their direction vectors are not coplanar), this system has a unique solution. Since a normal is by definition a unit vector, you then normalize this solution to obtain $\mathbf{n}$ (assuming unit albedo, or equivalently absorbing albedo into the scale before normalization).

---

1. Background Knowledge

• Lambertian reflectance model For a Lambertian (matte) surface, the intensity under a directional light $\mathbf{l}_i$ is

$$I_i = \rho \max(0, \mathbf{l}_i^\top \mathbf{n})$$

where $\rho$ is albedo, $\mathbf{n}$ is the unit surface normal, and $\mathbf{l}_i$ is usually taken as a known unit vector in lighting direction.

• Linear system formulation With three lights, stack the equations:

$$\begin{bmatrix} I_1 \\ I_2 \\ I_3 \end{bmatrix} = \begin{bmatrix} - \\ \mathbf{l}_1^\top \\ - \\ - \\ \mathbf{l}_2^\top \\ - \\ - \\ \mathbf{l}_3^\top \\ - \end{bmatrix} \cdot (\rho \mathbf{n}) \quad\Longrightarrow\quad I = L g,\; g = \rho \mathbf{n}$$

where $L$ is $3 \times 3$. If $L$ is invertible, $g = L^{-1} I$. This $g$ is a 3D vector pointing in the normal direction but scaled by albedo. Normalizing it removes the unknown scale.

---

2. Algorithm / Approach

At a high level, the algorithm pattern is:

1. Precompute light matrix inverse: from known light directions, build matrix $L$ (each row is a light direction) and compute its inverse $L^{-1}$.
2. Per-pixel solving: for each pixel, collect its observed intensities under the three lights into a vector $I$.
3. Solve for scaled normal: compute $g = L^{-1} I$, which is proportional to $\mathbf{n}$.
4. Normalize: set $\mathbf{n} = g / \|g\|$ to obtain a unit normal (ignoring absolute albedo).

This is a direct linear algebra solution; no iteration or optimization is needed in the ideal noiseless three-light case.

---

3. Step-by-Step Strategy

1. Form the light direction matrix $L$

• You are given three light directions $\mathbf{l}_1,$\mathbf{l}$_2,$\mathbf{l}$_3 \in \mathbb{R}^3$.
• Ensure each is a row of $L$:

import numpy as np

L = np.array([
l1, # shape (3,)
l2,
l3
]) # shape (3, 3)

2. Check invertibility of $L$

• Compute determinant or rank:

if np.linalg.matrix_rank(L) < 3:
raise ValueError("Lights are coplanar / not linearly independent")

• Optionally precompute inverse:

L_inv = np.linalg.inv(L)

3. For each pixel, gather intensity vector $I$

• Suppose you have three grayscale images I1, I2, I3 (same size H×W).
• At pixel $(y, x)$:

I = np.array([I1[y, x], I2[y, x], I3[y, x]]) # shape (3,)

4. Solve for $g = \rho \mathbf{n}$

• Multiply with $L^{-1}$:

g = L_inv @ I # shape (3,)

5. Normalize to get the unit normal $\mathbf{n}$

• Compute norm and normalize:

norm = np.linalg.norm(g)
if norm > 1e-8: # avoid division by zero
n = g / norm
else:
# handle degenerate / very dark pixel
n = np.array([0.0, 0.0, 0.0]) # or any convention

6. Vectorize for efficiency (optional)

• Instead of looping pixels, stack intensities into shape (3, N) and multiply once with L_inv, then normalize each column.

---

4. Common Pitfalls

• Non-invertible or ill-conditioned light matrix

• If lights are nearly coplanar or have very similar directions, $L$ becomes ill-conditioned, amplifying noise. Even if the determinant is non-zero, numerical errors may be large.

• In a contest setting, ensure test lights are well-separated; in general code, you could check the condition number of $L$.

• Forgetting normalization

• $L^{-1} I$ gives $\rho \mathbf{n}$, not a unit normal. Forgetting to normalize yields incorrect normals and angles.

• Zero or near-zero intensities

• If all intensities are very small (e.g., point is in shadow or background), the norm of $g$ will be tiny, leading to unstable normalization. Handle this explicitly (e.g., output a sentinel normal or skip).

• Incorrect matrix layout (rows vs columns)

• The formula assumes each row of $L$ is a light direction. If you treat them as columns, you must adjust the equations accordingly. Inconsistent layout leads to wrong results or dimension mismatch.

• Non-unit light directions

• If your light directions are not normalized, the proportionality between intensity and dot product changes. Either normalize $\mathbf{l}_i$ beforehand or be consistent with your model; for this simple problem, lights are typically given as unit vectors.

---

5. Time & Space Complexity

Let there be $N$ pixels (e.g., $N = H \times W$):

• Precomputation

• Inverting $L$ (a $3 \times 3$ matrix) is $O(1)$ with a small constant.

• Per-pixel computation

• Each pixel: one matrix-vector multiplication (3×3 by 3×1) and one normalization.

• This is $O(1)$ per pixel, so total time is $O(N)$.

• Space complexity

• Storing three input images: $O(N)$.

• Storing the normal map (3 channels per pixel): $O(N)$.

• The matrices $L$ and $L^{-1}$ are constant size.

• Overall space is $O(N)$.