I implemented the Ikeuchi and Horn shape-from-shading algorithm for the final
project in one of my computer vision courses at Berkeley.
This is an iterative method for computing
surface normals from the image of a shaded 3D surface. The algorithm
assumes the surface reflectance properties, the scene lighting, and the viewpoint
are all known. The algorithm finds a solution by minimizing an energy function based
on a brightness constraint and a smoothness constraint. The brightness
constraint ensures that the resulting surface orientations give rise to
image luminance values which closely match the original input image. A
reflectance map based on the known lighting and surface properties relates
surface orientation to image luminance. The smoothness constraint ensures
that the surface orientations of neighboring points will vary smoothly across
the surface. In order to guarantee
that the algorithm will converge to a solution, the
known orientations at the occluding boundary serve as a final
constraint on the system. My implementation of this algorithm takes
grayscale images as input and produces a set of surface normals as
output.
I then used a simple least-squares method to approximate a 3D surface
which fits the resulting normals.

Ikeuchi and Horn represent surface orientations
using a stereographic projection of the Gaussian sphere.
Every surface orientation visible to a viewer can
be represented on a Gaussian hemisphere. The 2D projection of this hemisphere allows
us to refer to surface orientation using the (f,g) coordinates of this orientation
space. The following image illustrates this projection (adapted from Ikeuchi
and Horn 1981).

The Ikeuchi and Horn algorithm assumes the
surface reflectance properties, scene lighting, and viewpoint are all known.
From this
information, we construct a reflectance map which relates images intensites to surface
orientation. For a given surface orientation (f,g), this reflectance map reveals what
the corresponding brightness value should be in the image. Here is an example
reflectance map for a lambertian sphere:

Although the
reflectance map relates image intensities to surface orientations, it does
not allow us to uniquely recover the surface orientation for any arbitrary
point in the image. Most brightness values in the image will be consistent with a
set of surface orientations, any of which could give rise to the brightness
value in question.
Here is an example reflectance map for a lambertian surface illuminated by a light source directed
from (0, 0, 1). The highlighted contour (in white) represents the surface
orientations consistent with an image intensity I(x,y) = 0.6.

Although
this constrains the solution significantly, additional constraints are
required before a final solution can be found.

The image irradiance constraint (also known as the brightness constraint)
refers to the requirement that the corresponding luminance
value R(f, g) for any estimated surface orientation (f, g) at position
(i, j) in the image should match the luminance value I(i, j)
as closely as possible. This constraint can be formally written as
follows:

The smoothness constraint minimizes total variation in surface orientation.
Given that an arbitrary image intensity value I(x, y) is
consistent with a set of possible surface orientations as described by
the reflectance function R(f, g), the smoothness constraint guarantees
that the surface orientations of the solution will smoothly vary
across surface points. This constraint can be formally written as
follows:

The algorithm iteratively updates the surface normals until the sum of the
error terms is minimized.
In order to guarantee the algorithm converges to a solution, the known
orientation of the surface at the occluding boundary is used to
further constrain the solution.

The following figures depict some example results using this algorithm.

References

Ikeuchi K and Horn BKP. Numerical shape from shading and occluding boundaries.
Artificial Intelligence. 17, Aug. 1981, 141-184.