The Horn and Brooks Algorithm

Suppose the local surface normal is n=(-p,-q,1)^T, where $p = \frac{\partial z}{\partial x}$ and $q = \frac{\partial z}{\partial y}$. For a light source at infinity, we can similarly write the light source direction as s =(-p_l,-q_l,1)^T. If the surface is Lambertian the reflectance map is given by

$$R(p,q) = {\mathbf n}\cdot{\mathbf s}$$ (1)

The image irradiance equation [12] states that the measured brightness of the image is proportional to the radiance at the corresponding point on the surface; that is, just the value of R(p,q) for p,q corresponding to the orientation of the surface. Normalizing both the image intensity, E(x,y), and the reflectance map, the constant of proportionality becomes unity, and the image irradiance equation is simply

E(x,y) = R(p,q) (2)

Although the image irradiance equation succinctly describes the mapping between the x,y co-ordinate space of the image and the the p,q gradient-space of the surface, it provides insufficient constraints for the unique recovery of the needle-map. To overcome this problem, a further constraint must be applied. Usually, the needle-map is assumed to vary smoothly. Another way to look at this constraint is that the goal of computation is to recover the smoothest surface satisfying the image irradiance equation.

The process of smooth surface recovery is posed as a variational problem in which a global error-functional is minimized through the iterative adjustment of the needle map. Surface normals are updated with a step-size dictated by Euler's equation. Here we consider the formulation of Brooks and Horn [5] which is couched in terms of unit surface normals. Their error functional is defined to be

 

$$\int\int \underbrace{\Big(E(x,y)-{\mathbf n}\cdot{\mathbf s})\Big... ...ial \mathbf n}\over{\partial y}}\right\Vert^2 \right)}_{Regularizing Term} dxdy$$ I = (3)

The terms ${{\partial \mathbf n}\over{\partial x}}$ and ${{\partial \mathbf n}\over{\partial y}}$ above are the directional derivatives of the needle-map in the x and y directions respectively. The magnitudes of these quantities are used to measure the smoothness of the surface, with a large value indicating a highly-curved region. A planar surface has all its normals parallel, so ${{\partial \mathbf n}\over{\partial x}}={{\partial \mathbf n}\over{\partial y}}=\mathbf 0$ in this case.

The functional of Equation 3 possesses two distinct terms. Firstly, the brightness error encourages data-closeness of the measured image intensity and the reflectance function. The regularizing term imposes the smoothness constraint on the recovered surface normals - it penalises large local changes in surface orientation. Departures from smoothness are measured by computing the magnitudes of the partial derivatives of the surface normals in the x and ydirections. The constant $\lambda$ is a Lagrange multiplier. For numerical stability, $\lambda$ must often be large, resulting in the smoothness term dominating.

Minimization of the functional defined in Equation 3 is accomplished by applying the calculus of variations [8] and solving the resulting Euler equation. The general Euler equation for a functional of the form in Equation 3 is

$$I_{\mathbf n} - {{\partial}\over{\partial x}}I_{\mathbf n_x} - {{\partial}\over{\partial y}}I_{\mathbf n_y} = {\mathbf 0}$$ (4)

Here, I_n is the result of taking the partial derivative of I with respect to the surface normal, n. Similarly, I_nx and I_ny are partial derivatives of I with respect to the directional partial derivatives of n in the x and y directions.

Substituting Equation 3 into the general Euler equation yields

$$\Big(E-{\mathbf n}\cdot{\mathbf s}\Big){\mathbf s} + \lambda\nabla^2{\mathbf n} - \mu{\mathbf n} = {\mathbf 0}$$ (5)

We discretize Equation 5 using the following approximation to the Laplacian

$$\Big\{\nabla ^2{\mathbf n}\Big\}_{i,j} \approx \frac{4}{\epsi... ...2}\Big(\bar{\mathbf n}_{i,j} - {\mathbf n}_{i,j}\Big)\nonumber$$

where $\epsilon$ is the spacing of pixel-sites on the lattice and $\bar{\mathbf{n}}_{i,j}$ is the local 4-neighbourhood mean of the surface normals around pixel position i,j

$$\bar{\mathbf{n}}_{i,j} = \frac{1}{4}\Big({\mathbf n}_{i,j+1} ... ...bf n}_{i,j-1} + {\mathbf n}_{i+1,j} + {\mathbf n}_{i-1,j}\Big)$$ (6)

Discretizing Equation 5 and re-arranging, we obtain an update equation to estimate the needle-map at epoch k+1 using the estimate at epoch k

$${\mathbf n}_{i,j}^{\left(k+1\right)}= \bar{\mathbf{n}}_{i,j}^... ...bf n}_{i,j}^{\left(k\right)}\cdot{\mathbf s}\right){\mathbf s}$$ (7)

Considering the form of Equation 8, it is clear that, at each iteration, there is a small step in the direction of the light source, coupled with a smoothing action using the local 4-neighbour mean.