We address the problem of estimating the three-dimensional shape and radiance of a surface in space from images obtained with different focal settings. We pose the problem as an infinite-dimensional optimization and seek for the global shape of the surface by numerically solving a partial differential equation (PDE). Our method has the advantage of being global (so that regularization can be imposed explicitly), efficient (we use level set methods to solve the PDE), and geometrically correct (we do not assume a shift-invariant imaging model, and therefore are not restricted to equifocal surfaces).