Previously we presented a shape reconstruction method from photometric stereo, which applies the Jacobi iterative method to reflectance map equations for M images and linearly combines the resulting iterative relations, to directly estimate the depth map of the object. For the case of two images, however; the method gives rise to noticeable distortions for certain lighting directions. In this paper, four approximations of the surface normal are introduced and the resulting 4M iterative relations are linearly combined as constraints, to effectively realize a symmetric discretization and achieve robust estimation free from such distortions. The method is investigated numerically using both synthetic and real images.