The integration of surface normals for the purpose of computing the shape of a surface in 3D space is a classic problem in computer vision.However,even nowadays it is still a challenging task to devise a method that is flexible enough to work on non-trivial computational domains with high accuracy,robustness,and computational efficiency.By uniting a classic approach for surface normal integration with mod computational techniques,we construct a solver that fulfils these requirements.Building upon the Poisson integration model,we use an iterative Krylov subspace solver as a core step in tackling the task.While such a method can be very efficient,it may only show its full potential when combined with suitable numerical preconditioning and problem-specific initialisation.We perform a thorough numerical study in order to identify an appropriate preconditioner for this purpose.To provide suitable initialisation,we compute this initial state using a recently developed fast marching integrator.Detailed numerical experiments illustrate the benefits of this novel combination.In addition,we show on real-world photometric stereo datasets that the developed numerical framework is flexible enough to tackle mod computer vision applications.