The general problem of surface matching is taken up 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 [12]. It also generalizes a 2D matching process developed in [4]. An 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 [17], 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.