An Algebraic Approach to Surface Reconstruction from Gradient Fields 
Abstract: Several important problems in computer vision such as Shape from Shading (SFS) and Photometric Stereo (PS) require reconstructing a surface from estimated gradient field, which is usually nonintegrable, i.e. have nonzero curl. We propose a purely algebraic approach to enforce integrability in discrete domain. We first show that enforcing integrability can be formulated as solving a single linear system $Ax=b$ over the image. In general, this system is underdetermined. We show conditions under which the system can be solved and a method to get to those conditions based on graph theory. The proposed approach is noniterative, has the important property of local error confinement and can be applied to several problems. Results on SFS and PS demonstrate the applicability of our method. 
Techniques
such
as Shape from Shading (SFS) and Photometric Stereo (PS) give a
nonintegrable gradient field as output, which needs to be integrated
to get the final shape. Previous approaches such as solving Poisson
equation and FrankotChellappa algorithm are based on least squares.
These techniques distribute error across the entire surface and does
not have the property of local error confinement. Thus errors in a part
of the gradient field are not restricted leading to smoothing of
features in the final shape. Our algorithm is based on two key observations. First, for an N*N image, there are N^2 pixels but 2*N^2 gradients. The gradient field is thus overdetermined and all gradients are not required for integration. This means that if one can identify erroneous gradients using some measure, we can just use the rest of gradients for integration, or we can recover the value of erroneous gradients. Second, the curl of the gradient field gives information about the errors in the gradient field. Least squares based approaches uses the divergence of the given gradient field for reconstruction but completely ignores the curl. Thus the basic idea is to identify erroneous gradients using curl and recover their values using integrability constraints. 
Let A1 be
the true
gradient
field and the obtained gradient field using Photometric Stereo or Shape
from Shading be A2. The gradient field A2 is nonintegrable (has
nonzero curl). Least square approaches such as solving Poisson
equation will give A3 as the solution, which has the same divergence as
the A2 and zero curl. We consider A2 as the sum of the true gradient
field and a residual gradient field. Using the curl, we try to estimate
the residual gradient field.

We
consider loops of 4 square connected pixels for calculating the curl.
If the curl is nonzero for the loop, the given gradients are
considered to be erroneous in that loop. Thus the residual gradients
are nonzero for that loop. We write the loop equation for all such
loops for which the curl is nonzero, giving a linear system of
equation Ax=b. To solve for the residual gradients x, we need to solve Ax=b. However, the matrix A should be full rank. To find whether the matrix is full rank or not, consider a 2D graph where edges represents the gradients and nodes represents the surface values. Suppose, we break all the edges corresponding to erroneous gradient values. The linear system Ax=b can be solved if the graph remains connected. If it breaks into n pieces, we need n1 edges to connect it in a single piece and rank deficit is n1. To connect the graph, we assign a weight to each edge and find that set of n1 edges which has the minimum weight. After connecting the graph, we again form Ax=b corresponding to the broken edges and solve for x. The estimated residual gradient field is then subtracted from the given gradient field. The final gradient field is integrated to obtain the surface. 
(a) A
sample graph
on a 8 by 8 grid. Blue lines denotes edges and red
dots denotes nodes. (b) A configuration such as this cannot happen as
boundary edges are known (c) A particular configuration of broken edges
where the graph is broken into n = 4 parts. total number of broken
links = 34, rank(A) = 31. rank deficit = 3431 = 3 = n1 (d) If 3
more links are known to connect the graph (for e.g. dashed, green
edges), number of unknowns reduces to 31 = rank(A)
Also see my ECCV 2006 paper for improved algorithms. 