There are various methods for denoising image signals. With a preface about image sparsity I'm going to provide some matlab code.
Images are sparse in nature, not in spatial domain, but in some domains as Fourier, DCT or Wavelets domain. If you take a look at FFT of an image you'll find out only a few coefficients are non-zero, as figure below depicts fft2 of your image:

imshow(log10(abs(fftshift(fft2(Image)))),[])
and in fact image compression algorithms like JPEG take advantage of this property of images. This Sparsity property also can be exploited for denoising images. If you transform an image into wavelet domains, the random noise like the above mentioned noise tends to spread out its coefficient, however coefficients related to the image tend to concentrate to only a few bins. Therefore, the most basic, simplest technique would be sorting FFT, DCT of Wavelet coefficients of image in order of magnitude, making zero all coefficient which are smaller than a threshold and finally transforming the image back into the spatial domain. Based on this idea, some soft/hard thresholding algorithms appeared. A good Matlab implementation on this regard can be find in following page, where they provide wavelet based denosing tools:
http://eeweb.poly.edu/iselesni/WaveletSoftware/denoise.html
Another technique to take advantage of sparsity is $l_1$ regularization. Look, it is proven under some assumptions, $l_1$ norm can be used as representation for number of nonzero coefficients ($l_0$ norm, you can find the proof in Compressive Sensing literate).
Note again, random noise tend to form a non-sparse representation in transform domain (e.g white noise has flat representation in frequency). So if we define a convex optimization problem in which we try to build an image which preserve energy of noisy image and also tries to minimize number of nonzero coefficients in transform domain we would remove noise. The problem usually defined as below:
$$D = \arg \min_d ||I-d||_2 + \lambda|| f(d) ||_1 $$
Where, $D$ is denoised image, $I$ is noisy image, $f$ is transformation operator (e.g. DFT matrix) and $\lambda$ is regularization factor. The first term tries to make $D$ closer and closer to original image, and the second term tried to keep $D$ as sparse as possible. Through solving the above problem, we would remove noise from the image.
For solving this optimization problem CVX toolbox might be used.
I am not sure that the code below in CVX toolbox would work, just try to give you the idea:
I=imread('J13Wn.jpg');
I=double(I);
Rn=I(:,:,1);
Rn=imresize(Rn,[256 256]); % for out of memory error in my computer
B=dctmtx(256);
cvx_begin
variable R(256,256);
minimize(norm(Rn-R,2)+0.1*norm(B*R,1));
cvx_end
if it did not work, try this toolbox:
http://www.cs.tut.fi/~foi/GCF-BM3D/