In this paper we demonstrate how to recover surface shape from single images of faces using shape-from-shading when shadows are present. We make use of a statistical representation of the distribution of surface normal directions based on the equidistant azimuthal projection. This is allows us to develop a statistical model of the variations in facial shape in the surface normal domain. We show how ideas from robust statistics can be used to fit the model to facial images in which there is significant self-shadowing. The method is evaluated on both synthetic and real-world images. It is demonstrated to effectively fill-in the facial surface when more than 30% of the area is subject to self-shadowing.