Suppose the local surface normal is
For a light source at infinity, we can similarly write the light source direction as
=(-pl,-ql,1)T. If the surface is Lambertian the reflectance map is given by
The image irradiance equation  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
minimized through the iterative adjustment of the needle map.
are updated with a step-size dictated by Euler's equation. Here we consider the formulation of Brooks and
Horn  which is couched in terms of unit
surface normals. Their error functional is defined to be
The terms and 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 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 is a Lagrange multiplier. For numerical stability, 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  and solving the resulting Euler equation. The general Euler equation for a functional of the form in Equation 3 is
Here, In is the result of taking the partial derivative of I with respect to the surface normal, n. Similarly, Inx and Iny 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
We discretize Equation 5 using the following approximation to the Laplacian
is the spacing of pixel-sites on the lattice and
is the local 4-neighbourhood mean of the surface normals around pixel position i,j
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
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.