Medial surfaces are popular representations of 3D objects in vision, graphics and geometric modeling. They capture relevant symmetries and part hierarchies and also allow for detailed differential geometric information to be recovered. However, exact algorithms for their computation from meshes must solve high-order polynomial equations, while approximation algorithms rarely guarantee soundness and completeness. In this article we develop a technique for computing the medial surface of an object with a polyhedral boundary, which is based on an analysis of the average outward flux of the gradient of its Euclidean distance function. This analysis leads to a coarse-to-fine algorithm implemented on a cubic lattice that reveals at each iteration the salient manifolds of the medial surface. We provide comparative results against a state-of-the-art method in the literature.