A new method for smoothing both gray-scale and color images is presented that relies on the heat diffusion equation on a graph. We represent the image pixel lattice using a weighted undirected graph. The edge weights of the graph are determined by the Gaussian weighted distances between local neighbouring windows. We then compute the associated Laplacian matrix (the degree matrix minus the adjacency matrix). Anisotropic diffusion across this weighted graph-structure with time is captured by the heat equation, and the solution, i.e. the heat kernel, is found by exponentiating the Laplacian eigen-system with time. Image smoothing is accomplished by convolving the heat kernel with the image, and its numerical implementation is realized by using the Krylov subspace technique. The method has the effect of smoothing within regions, but does not blur region boundaries. We also demonstrate the relationship between our method, standard diffusion-based PDEs, Fourier domain signal processing and spectral clustering. Experiments and comparisons on standard images illustrate the effectiveness of the method.