  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-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.
