Since this is apparently now a canonical question, and it hasn't been mentioned here yet:
One natural extension of k-means to use distance metrics other than the standard Euclidean distance on $\mathbb R^d$ is to use the kernel trick. This refers to the idea of implicitly mapping the inputs to a high-, or infinite-, dimensional Hilbert space, where distances correspond to the distance function we want to use, and run the algorithm there. That is, letting $\varphi : \mathbb R^p \to \mathcal H$ be some feature map such that the desired metric $d$ can be written $d(x, y) = \lVert \varphi(x) - \varphi(y) \rVert_{\mathcal H}$, we run k-means on the points $\{ \varphi(x_i) \}$. In many cases, we can't compute the map $\varphi$ explicitly, but we can compute the kernel $k(x, y) = \langle \varphi(x), \varphi(y) \rangle_{\mathcal H}$. Not all distance metrics fit this model, but many do, and there are such functions defined on strings, graphs, images, probability distributions, and more....
In this situation, in the standard (Lloyd's) k-means algorithm, we can assign easily points to their clusters, but we represent the cluster centers implicitly (as linear combinations of the input points in Hilbert space). Finding the best representation in the input space would require finding a Fréchet mean, which is quite expensive. So it's easy to get cluster assignments with a kernel, harder to get the means.
The following paper discusses this algorithm, and relates it to spectral clustering:
I. Dhillon, Y. Guan, and B. Kulis. Kernel k-means, Spectral Clustering and Normalized Cuts. KDD 2005.