We present the first numerical implementation of a non-perturbative renormalization method for lattice operators, based on the study of correlation functions in coordinate space at short Euclidean distance. The method is applied to compute the renormalization constants of bilinear quark operators for the non-perturbative O(a)-improved Wilson action in the quenched approximation. The matching with perturbative schemes, such as (MS) over bar, is computed at the next-to-leading order in continuum perturbation theory. A feasibility study of this technique with Neuberger fermions is also presented. (C) 2004 Published by Elsevier B.V.