An Algebraic Approach to Surface Reconstruction from Gradient Fields
Amit Agrawal, Rama Chellappa & Ramesh Raskar

ICCV 2005
pdf Matlab Code


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 non-integrable, i.e. have non-zero 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 under-determined. 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 non-iterative, 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 non-integrable gradient field as output, which needs to be integrated to get the final shape. Previous approaches such as solving Poisson equation and Frankot-Chellappa 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 over-determined 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 non-integrable (has non-zero 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 non-zero for the loop, the given gradients are considered to be erroneous in that loop. Thus the residual gradients are non-zero for that loop. We write the loop equation for all such loops for which the curl is non-zero, 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 n-1 edges to connect it in a single piece and rank deficit is n-1. To connect the graph, we assign a weight to each edge and find that set of n-1 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 = 34-31  = 3 = n-1 (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.