The general problem of surface matching is considered in this study. The process described in this work hinges on a geodesic distance equation for a family of surfaces embedded in the graph of a cost function. The cost function represents the geometrical matching criterion between the two 3D surfaces. This graph is a hypersurface in 4-dimensional space, and the theory presented herein is a generalization of the geodesic curve evolution method introduced by R. Kimmel et al. It also generalizes the 2D matching process developed by I. Cohen et al. A Eulerian level-set formulation of the geodesic surface evolution is also used, leading to a numerical scheme for solving partial differential equations originating from hyperbolic conservation laws proposed by Sethian, which has proven to be very robust and stable. The method is applied on examples showing both small and large deformations, and arbitrary topological changes.