I'd like to welcome back guest blogger Stan Reeves, professor of Electrical and Computer Engineering at Auburn University. Stan will be writing a few blogs here about image deblurring.
In my last blog, I looked at image deblurring using an inverse filter and some variations. The inverse filter does a terrible job due to the fact that it divides in the frequency domain by numbers that are very small, which amplifies any observation noise in the image. In this blog, I'll look at a better approach, based on the Wiener filter.
I will continue to assume a shift-invariant blur model with independent additive noise, which is given by
y(m,n) = h(m,n)*x(m,n) + u(m,n)
where * is 2-D convolution, h(m,n) is the point-spread function (PSF), f(m,n) is the original image, and u(m,n) is noise.
The Wiener filter can be understood better in the frequency domain. Suppose we want to design a frequency-domain filter G(k,l) so that the restored image is given by
We can choose G(k,l) so that we minimize
E[] is the expected value of the expression. We assume that both the noise and the signal are random processes and are independent of one another. The minimizer of this expression is given by
This filter gives the minimum mean-square error estimate of X(k,l). Now that sounds really great until we realize that we must supply the signal and noise power spectra (or equivalently the autocorrelation functions of the signal and noise).
The noise power spectrum is fairly easy. We usually assume white noise, which makes the noise power spectrum a constant. But how do we determine what this constant is? If
then the noise power spectrum is given by
for an MxN image. This noise variance may be known based on knowledge of the image acquisition process or may be estimated from the local variance of a smooth region of the image.
The signal power spectrum is a little more challenging in principle, since it is not flat. However, we have two factors working in our favor: 1) most images have fairly similar power spectra, and 2) the Wiener filter is insensitive to small variations in the signal power spectrum.
Consider two very different images -- the cameraman and house.
cam = im2double(imread('cameraman.tif')); house_url = 'http://blogs.mathworks.com/images/steve/180/house.tif'; house = im2double(imread(house_url)); subplot(1,2,1) imshow(cam) colormap(gray(256)) title('cameraman') subplot(1,2,2) imshow(house) title('house')
We show the actual (log) spectrum for these two images.
cf = abs(fft2(cam)).^2; hf = abs(fft2(house)).^2; subplot(1,2,1) surf([-127:128]/128,[-127:128]/128,log(fftshift(cf)+1e-6)) shading interp, colormap gray title('cameraman power spectrum') subplot(1,2,2) surf([-127:128]/128,[-127:128]/128,log(fftshift(hf)+1e-6)) shading interp, colormap gray title('house power spectrum')
Let's see what happens if we restore the cameraman with the actual power spectrum.
h = fspecial('disk',4); cam_blur = imfilter(cam,h,'circular'); % 40 dB PSNR sigma_u = 10^(-40/20)*abs(1-0); cam_blur_noise = cam_blur + sigma_u*randn(size(cam_blur)); cam_wnr = deconvwnr(cam_blur_noise,h,numel(cam)*sigma_u^2./cf); subplot(1,1,1) imshow(cam_wnr) colormap(gray(256)) title('restored image with exact spectrum')
For comparison purposes, we restore the cameraman using the power spectrum obtained from the house image.
cam_wnr2 = deconvwnr(cam_blur_noise,h,numel(cam)*sigma_u^2./hf);
imshow(cam_wnr2)
colormap(gray(256))
title('restored image with house spectrum')
Visually, the two are very similar in quality. In terms of mean-square error (MSE), the former is better (lower), as the theory predicts:
format short e mse1 = mean((cam(:)-cam_wnr(:)).^2) mse2 = mean((cam(:)-cam_wnr2(:)).^2)
mse1 = 2.9905e-003 mse2 = 3.9191e-003
In both cases, the spectrum is concentrated in the low frequencies and falls away at higher frequencies. A smoothed version of the spectra would look even more similar. Rather than using the power spectrum from a specific image, one can either average a large number of images or use a simple model of the power spectrum or autocorrelation function.
A common model for the image autocorrelation function is
We calculate these parameters for the cameraman, averaging over horizontal and vertical shifts to get a single parameter for the correlation coefficient. Then we calculate an autocorrelation function.
sigma2_x = var(cam(:)) mean_x = mean(cam(:)) cam_r = circshift(cam,[1 0]); cam_c = circshift(cam,[0 1]); rho_mat = corrcoef([cam(:); cam(:)],[cam_r(:); cam_c(:)]) rho = rho_mat(1,2); [rr,cc] = ndgrid([-128:127],[-128:127]); r_x = sigma2_x*rho.^sqrt(rr.^2+cc.^2) + mean_x^2; surf([-128:127],[-128:127],r_x) axis tight shading interp, camlight, colormap jet title('image autocorrelation model approximation')
sigma2_x = 5.9769e-002 mean_x = 4.6559e-001 rho_mat = 1.0000e+000 9.4492e-001 9.4492e-001 1.0000e+000
From this we calculate another restored image:
cam_wnr3 = deconvwnr(cam_blur_noise,h,sigma_u^2,r_x);
imshow(cam_wnr3)
colormap(gray(256))
title('restored image using correlation model')
mse3 = mean((cam(:)-cam_wnr3(:)).^2)mse3 = 3.5255e-003
The restored image from this method is better than the restored image using the house spectrum but not quite as good as the one using the exact cameraman spectrum. All in all, we can see that the exact noise-to-signal spectrum isn't necessary to yield good results.
- by Stan Reeves, Department of Electrical and Computer Engineering, Auburn University
Get
the MATLAB code
Published with MATLAB® 7.5
The very first equation is not consistent with the following text. It should read:
y(m,n) = h(m,n)*x(m,n) + u(m,n)
where x is original image and y is degraded image.
Tomash—Thanks for pointing out the problem. I’ve fixed it.
I attempted to use Wiener filter to deconv a degraded image. But the psf is imaginary, not real. For deconvwnr function, it requires psf be real. Is it possible to use this function? In my case, psf is exp(j(x^2+y^2)). If deconvwnr function is not useful, could you recommend some ways to deconv it? Thank you.
Zhang—Your point-spread function has infinite extent and does not decay. How does that make sense physically?
Yes. You are right. The exact psf is exp(j(x^2+y^2)) convolved by a Gaussian envelope, i.e., exp(-k(x^2+y^2)). Because the envelope is what I can control all the time, I can neglect it and try to deconv exp(j(x^2+y^2)). If I use the convolution to work as psf, it is still imaginary, so I can’t use deconvwnr function.
Zhang—I believe you probably modify deconvwnr to accept a complex-valued PSF. However, I doubt that you will be able to obtain meaningful results for that PSF.
Zhang–This is interesting. I just recently was puzzling over an almost identical PSF! I would be interested to know the application in which your PSF arises. Actually, I think you don’t need a Wiener filter. This PSF is an all-pass filter, so you should just take the FFT of the sampled version and divide by it in the frequency domain. Since it’s all-pass, you don’t need to worry about noise amplification, which is what the Wiener filter is intended to combat. I’d be interested to hear how it goes.
Stanley－Thanks for your interest. The psf is used in the optics, to be exact, holography. You mentioned the difficulty of this reconstruction. Noise in my case is a little different. It is the defocus noise, i.e. the noise is also a convolution between certain signal and the psf. So the imaging procedure can be modeled, y = conv(x1, psf1) + conv(x2, psf2). psf = k*exp(jz(x^2＋y^2)). Here k and z are constants. psf1 and psf2 are different in z. I am trying to reconstruct x1 and diminish x2. So noise is equal to conv(x2, psf2). Wiener filter is the first method I wanna try, because it is simple and optimal to implement for deconvolving. But now it doesn’t work well to diminish the noise.
Zhang–Thanks for the details. It sounds like a very interesting application. The Wiener filter isn’t necessary just to deal with the PSF exp(j(x^2+y^2)). If you want to first use an inverse filter with this PSF like I described in my previous comment, then you can apply the Wiener filter using the Gaussian PSF and a noise spectrum determined from conv(x2, psf2) to deblur and reduce noise. After the inverse filter with exp(j(x^2+y^2)), the noise spectrum will not be changed, since this inverse filter only affects the phase of the signal and not the magnitude. Then you can apply the Wiener filter.
You can do noise filtering with a Wiener filter even if the PSF is all-pass. It’s just that you have to define the noise spectrum so that it reflects the actual spectrum of conv(x2, psf2).
Do you have a reference that explains the image formation equation?
Stan-Sorry for my late reply. The method You mentioned is really good at a reconstruction of x1 in terms of efficiency and computation. Now through a multiplication in frequency domain with a conjugate PSF, I can obtain a reconstructed result of x1. But the result is still contaminated seriously in some cases by x2. After applying the conjugate PSF, a model can be expressed as {x1+conv(conv(x2, psf2), conj(psf1)}. Though in some sense, x2 has been suppressed owing to the convolution with conj(psf1), it still influces negatively the final result, especially in certain cases. So at present, how to suppress x2 or the last part effectively in the model is a key to solve this problem.
The second point you mentioned is what I am thinking now. If I can put up with an appropriate estimation of noise, it will be convinient to remove it. I am struggling in the direction, but haven’t got an idea. Do you have some advice?
Reference for my work is mainly in the book:
Contemporary optical image processing with MATLAB
In section 3.2, it tells how the point spread function comes into being.
Million thanks
Hi Steve
I’m trying to correlate two images aplying normxcorr2. I read that the template is typically smaller than the image f(x,y), it is correct correlate images with the same size?
Thanks
Hilda - Yes, you can use NORMXCORR2 on images having the same size. Note that the normalization is not symmetric with respect to the two input arguments, so you may get different results depending on which image you pass as the TEMPLATE and which you pass as A.
Hi Steve,
I have just started working on depth from defocus. I am trying to find difference of sigma squares for two defocused images. Thus, I firstly created two blurred images with two known gaussian filters. However, I1f./If2 does not seem to me like the original equation exp(-0.5*(sigma1^2-sigma2^2)*(fx^2+fy^2)) but very degraded version of it, although there is no noise added.
What may it be the reason for that? I avoid of dividing by zeros in I2f. How could I find the sigma difference between two images?
Thanks a lot.
Thank you Mara!!
Hakan - I’m not entirely clear on the steps you’ve taken, but I’ll take a shot at a response. First, you blurred the same image with two different Gaussians, right? Then you took the FFT of both and divided one by the other and then took the IFFT of the result. Is that right? If that’s what you did and the result doesn’t look like a Gaussian, then it may be that the denominator — though not zero — is so small that you’re getting errors due to roundoff noise. If you didn’t do the IFFT step, it still should look like a Gaussian, at least around the low frequencies. I can’t guess anything else without having a clearer picture of what you’re doing.
Hello,
i have a question. i would like to avoid the ringing effect in the deblured image using edgetaper(). however, my transfer function is the same size a the image, and this function requires half the size for the otf. shoud i downsample my H? or upsample the image? thx
I tried to do the opposite - to estimate the noise variance from the blured image.
i did like this:
1.cropped the upper left corner of the cameraman image
2.
h=imhist(cam);
[v,unv]=statmoments(h,2);
i got for unv(1)=164 for the mean and unv(2)=46 for the variance.
my question is how from this result i obtain the orignal
sigma_u = 10^(-40/20)*abs(1-0) ?
==================================================
function statmoments.m:
function [v,unv]=statmoments(p,n)
lp = length(p);
G=lp-1;
p=p/sum(p);
p=p(:);
z=0:G;
m=z*p;
z=z-m;
v=zeros(1,n);
v(1)=m;
for j=2:n
v(j)=(z.^j)*p;
end
if nargout>1
unv=zeros(1,n);
unv(1)=m.*G;
for j=2:n
unv(j)=((z*G).^j)*p;
end
end
Assaf–Regarding the estimate of the noise variance, the way you propose won’t work well. It overestimates the noise variance, because it assumes a constant mean when the signal actually varies slowly. Instead, you need to compute a pointwise local variance. One way to do this is to compute a local mean by using the filter function, then subtracting the local mean from each pixel value, and averaging the result. Something like this should work reasonably well:
local_mean_p = filter2(ones(3)/9,p,’valid’);
p_valid = p(2:end-1,2:end-1);
var = mean((p_valid(:)-local_mean_p(:)).^2);
Note that you need to exclude the boundaries of the local mean and data from the calculation, because the filter2 function averages in zero boundaries and throws off the mean calculation.
Assaf—Following up on Stan’s comment, you might want to look at the Image Processing Toolbox function stdfilt to compute your pointwise local variance. This function uses symmetric extension at the image boundaries.
Assaf–If your PSF is more than half the size of the image, then one of two things are probably true. Either you’re not going to get very good results because the image just isn’t big enough for the spreading caused by the PSF, or the PSF is mostly concentrated in the middle. If it’s the latter, then you might try windowing the PSF to half the size of the image. It may be best to window the PSF and then taper the result using a 2-D window. Something like:
x_taper = x.*(hamming(N)*hamming(N)’);
It won’t work to interpolate the image or decimate the PSF, because then you’ll mess up the spatial scaling correspondence between the image and the PSF.
If the PSF is symmetric, you can try symmetrically extending the image in both dimensions. Even if it’s not symmetric, that often works fairly well.
Thank you very much Stan and Steve. I would like to tell you the whole picture because i think i’m doing something wrong. I want to deblurr an image using the wiener filter. the image is effected by turbulence. the OTF is exp(-lamda*(u^2+v^2)^(5/6)) and the size of the picture is 576*720. I wrote this program:
0. Blurred=imread(’blurred.tif’)
1. lamda=0.001
2. [u,v]=computed_size_with_meshgrid(576,720)
3. H=exp(-lamda.*(u.^2+v.^2).^(5/6))
4. h=real(ifftshif(ifft2(H)))
5. restored= deconvwnr(Blurred,h,1)
I get result that looks better than the blured image, but i’m sure more can be achived.
My questions:
1. Is the size of the OTF should be the same as the size of the picture in line 3?
2. Is the size of the PSF should be the same as the size of the OTF in line 4?
3. Do i convert the OTF to PSF correctly in line 4? or should i use other functions like fsamp2(), or psf2otf()?
4. If i want to use edgetaper(), how do i make the PSF be half the size of the picture, if my OTF is the same size of the picture to begin with?
5. Can you please correct my program?
Thank you very much
Assaf–Let me suggest a couple of things. First, H must be symmetric around the origin (thinking of H as one period of a periodic function). The way you’re computing it, you’re only getting one side of H. Regarding #4, I think I addressed that in Comment 20 above.
Stan - thank you! i changed my H so it be symmetric , and the results are much better! i didnt know how to make it symetric, so i looked at the help of function fspecial() and it says it produces symmetric filters… so i used fspecial(’gaussian’) insted :)
Hello, regarding the estimate of noise variance at smooth regions of the frame, how can i estimate where the smooth regions are automaticlly?
is using functions rangefilt() a good idea?
thx
Assaf—I don’t think it’s a good idea to use rangefilt, because I imagine that operator is sensitive to noise. You might try something like sorting the estimated noise variances, taking some fraction of the lowest values, and then taking the median.
Hello - if i implement the wiener filter formula as is, a.k.a h(k,l)/(…) will i get the same result as deconvwnr?
Assaf—You can look inside deconvwnr to see exactly the computation it uses.
Hey, i looked but i have a problem:
the deconvwnr() usess fftn and it accepts any dimensions, whereas in simulink i am using 2d-fft and it tells me it can only work with dimensions that are power of 2..
what should i do?
Assaf—You’ll have to zero-pad your arrays.
Will it not cause me a problems is the edges and theres a chance i get ringing?
Assaf—Ringing at the image boundaries is always a potential issue with FFT-based deblurring. You can reduce the problem by applying a spatial-domain window to the image.
can you please reffer me to a metarial regarding the subject?
Assaf—This paper is useful: Stanley J. Reeves, “Fast restoration of PMMW imagery without boundary artifacts,” in SPIE Vol. 4719 — Infrared and Passive Millimeter-Wave Imaging Systems: Design, Analysis, Modeling, and Testing (R. Appleby, G. C. Holst, and D. A. Wikner, eds.), (Orlando, FL), pp. 289-295, SPIE - Int. Soc. Opt. Eng. (US), April 2002. You can probably get Dr. Reeves to send you a copy. He’s at Auburn University. A more general reference is the book Deblurring Images: Matrices, Spectra, and Filtering. It’s got examples and code based on MATLAB. I’ve only thumbed through it briefly, but it looked interesting.
Thank you very much
Hello,
Congratulations for the blog.
Can you give me a reference of model for the image autocorrelation function?
Thanks a lot.
Ernesto Moya.
Ernesto—I don’t understand your question.
Excuse my english, it is not my native languaje, where I can find the explication del model of autocorrelation used?, any book or paper, you say “…A common model for the image autocorrelation function is…”, where can I read more for understand the model.
Thansk for your time.
Ernesto—Try Fundamentals of Digital Image Processing, by Anil K. Jain.
Thank you very much.
Hi,
i want to know how PSF can be constructed for a motion blur. For a particular length and angle what will be the size of psf . What will be the function for this?? will it be symmetric function for a linear motion.
Hi Steeves ,
i have a question. if the psf is symmetric the phase response function should be zero. But in matlab i have generated the motion blur PSF using fspecial function which comes to symmetric in some cases. but phase is not zero in this case
Ratnakar—I assume you are computing the phase response by calling fft2, or something similar. The function fft2 computes the two-dimensional discrete Fourier transform (DFT). The symmetry properties of the DFT may not be what you are accustomed to. For example, for the one-dimensional sequence x[n], 0 < = n <= N-1, the DFT X[k] is real if x[n] = x*[N - n]. This is sometimes called "circularly symmetric," and sometimes written using mod: x[n] = x*[mod(-n, N)].
You should zero pad your PSF according to the size of the 2-D FFT you want to use, and then you should use circshift to circularly shift the PSF so that it’s “center” is located at the upper left corner of the matrix. Then call fft2. For example, see the implementation of psf2otf.
Ratnakar—You can use fspecial to construct a linear motion PSF.
What does the function “deconvwnr” assume about the noise to signal power (K) when it is called by the following command:
deconvwnr(blurredImage, psf)?
In other words, the NSR or noise/image correlation aren’t specified in the command line. If K=0, then the reconstructed image would be blurredImage./psf .
Thanks for the reply..
I have another question if you don’t mind. if i use fspecial to generate a motion PSF, and then use imfilter like below
f = imread(’cameraman.tif’);
h = fspecial(’motion’,7,45);
g = imfilter(f,h,’circular’);
G = fft2(g);
if i take the phase spectrum of g then it should be equal to phase spectrum of f. But i am not getting it? where is the problem.. can u help me in solve this problem of phase.
The question i have just given is for symmetrical PSF
J.—As noted in the documentation, the default value assumed for NSR is 0.
Ratnakar—There are two issues. The first is that imfilter produces a uint8 output image if the input is uint8. Quantizing the output image to uint8 destroys the phase equality you are expecting. Try this:
f = im2double(imread('cameraman.tif'));The second issue is that the frequency response of your filter has zeros. Floating-point round-off effects will prevent an exactly zero-phase response wherever the frequency response is zero or nearly zero.
Thanks for the reply ..
i have a doubt. How to calculate the SNR from a noisy image. You are adding PSNR to the blurred image . But , most of the literatures writing that 30 dB SNR , 40dB SNR etc.
please make me understand this noise addition to image and how to calculate SNR from a burred noisy image or noisy image
Hi,
I have another question. Can any one name some features that change with blurring ? The features should be such that it is applicable for all standard images
Ratnakar—A Google search on the terms “image PSNR” should get you the needed formulas.
Ratnakar—Just about any non-constant region of an image will change with blurring.
Hi Steve,
I was wondering whether anyone has witnessed the strange dependence of deconvwnr.m on the overall signal intensity range of the image.
I am using deconvwnr.m to deconvolve a known blur to the standard cameraman.tif matlab image. I do not add any noise. The initial result is good. However if I then scale the image by x100, apply the same blur, and then deconvolve using deconvwnr, the result becomes more blurred. I have seen the same behaviour with other images, including a simple checkerboard image.
I have tried changing the NSR (even though it should be 0 theoretically as I’ve added no additional noise) and it does not make any difference, just exacerbates the blurriness!
I have now written my own deconvolution routine using a Wiener Filter and do not get this behaviour. But I would love to find out why this matlab routine shows this strange dependence.
Many thanks for your help,
Joan
Joan—I have reproduced the effect you describe and will investigate.
Joan—I apologize for the problem you encountered. My tentative conclusion is that there’s a bug in this line of deconvwnr:
It should instead read:
I’ll continue to investigate. When an official patch is available on our tech support site, I’ll post again.
Thanks for taking the time to let us know about the problem.
Joan—There is now a bug report, including a work-around, available on our tech support site.
Hi,
I have 3 questions regarding the following sentence:
“Rather than using the power spectrum from a specific image, one can either average a large number of images or use a simple model…”
My questions are:
1. is the sigma2_x= the pointwise local variance for the image? (comupted with stdfilt for example)
2. is the sigma2_x=the pointwise local variance of average a large number of images?
3.is the power spectrum can be calculated like this:
PS=abs(fft2(avarage_image).^2)?
thanks..
Assaf–here are some answers:
1. Yes, if that’s the way you want to define the random process that created the image.
1. Yes, if that’s the way you want to define tathe random process that created the image. Either this or #1 is ok, just two different ways to do it.
3. Yes, that’s one way to do it. The idea of a power spectrum of an image is that the image is a sample of a random process. So, in effect you’d be computing the power spectrum from one sample. If you assume that the process is stationary, then you can average spectral estimates over blocks of the image. If you assume that all images are samples of the same random process, you can average the spectrum over many images.
Dear Sir
I want to restore a vibrated image by
applying weiner filter. I dont know how to apply the filter to the vibrated image. can you help?
thank you
david
David—Deconvolution using a Wiener filter requires that you provide the point-spread function (PSF) for the blurring process. If you have the PSF, you can use deconvwnr.
Hi
when I am doing this surf([-127:128]/128,[-127:128]/128,log(fftshift(cf)+1e-6))
It says data dimensions must agree and size of this one log(fftshift(cf)+1e-6)) is 1X2? I dont know whats the problem
is
Is there an elegant way to preserve radiometric calibration while applying a Wiener filter to very noisy images? I want to preserve the absolute radiometric intensity of faint stars in a very noisy image, where PSDNoise>PSDSignal. This results in a Wiener filter W < 1, sometimes substantially so, at all spatial frequencies. The Wiener filter indeed reveals faint stars buried in the noise, but the absolute intensities are grossly underpredicted.
Bhanuprkash—If log(fftshift(cf)+1e-6)) is 1-by-2, then cf must be 1-by-2. cf is computed from cam, which in turn is read from a 256-by-256 image cameraman.tif. So cf should be 256-by-256. Since it’s not, I’d guess that you have a bug somewhere in your code. Try clicking on the link at the bottom of the post that says “Get the MATLAB code”.
Bob–—The Wiener filter introduces bias into the solution to reduce the noise variance. Essentially, it allows for a small amount of well-defined blurring to smooth the noise while still restoring the image. Thus, impulsive signals like stars will be spread out somewhat and lose their original amplitude.
You can calculate the effective PSF of the blurring followed by Wiener filter and then infer the height of an impulse that has been blurred by this PSF just by comparing the PSF height to the height of the “bump” caused by the star. However, this still won’t recover the original radiometric intensity unless you either know or can estimate the original radius of the star. As you probably know better than I do, if the resolution is fairly limited, you won’t be able to tell the difference between a tiny but bright star and a larger but dim one—especially if the star size is smaller than a pixel in the image plane.
Hi,
I have a question regarding the calculations of image autocorrelation function, in our example:
cam_r = circshift(cam,[1 0]);
cam_c = circshift(cam,[0 1]);
rho_mat = corrcoef([cam(:); cam(:)],[cam_r(:); cam_c(:)])
rho = rho_mat(1,2);
[rr,cc] = ndgrid([-128:127],[-128:127]);
is there an easier way o calculate the rho? i have explored the corrcoef() function and it invloves many calculations… can i calculate in another, easier way?
thx
hi
what does the different between inverse filter and wiener filter ?
Hi, i have some images which are blurred due to circular motion. The different images with angular speed are given.
Can any one help me in removing the blur ??
1.The minimizer of this expression is given by
G(k,l)=H(k,l)/(H(k,l)^2+S_u(k,l)/S_x(k,l))
2.A common model for the image autocorrelation function is
sigma_x*rho.^sqrt(m^2+n^2) + mean_x^2
I am astudent,
I Did not know how it does come,
Ask how to infer this formula?
Hi.
I found a paper on fingerprint enhancement. In that paper there was Wiener filter with formula given as
w(n1,n2)= u+((s-v)/s)*(I(n1,n2)-u)
where:
v = noise variance
s = local variance
u = local mean
I = grey-level intensity
and everything is for 3×3 neighbourhood.
i’ve found everything else but don’t know how to find out v(noise variance).
Please help in this regard.
Thanks in advance
Assaf—Image autocorrelation is related to the image power spectrum, so you can probably use FFTs to speed up the calculation. This is a topic you can probably find in a book on multidimensional digital signal processing.
Ali—Wiener filtering simplifies to a pure inverse filter when you assume there is no noise. In the presence of noise, a pure inverse filter will tend to greatly amplify any noise present, often making the result unusable.
Meng-cian—1. I’ll leave it to you to work out the math on that one, or you can ask one of your instructors. 2. I’m going from a very old memory here, so you’ll need to verify this, but I believe that form of autocorrelation model comes from the assumption that the signal can be approximated by a first-order autoregressive model. You might be able to find more information in a reference on spectrum analysis.
Umar—I’m not an expert on this subject, but I’m only aware of fairly ad hoc methods for estimating noise levels in an image. For example, you could compute the pixel value variance for a subset of image pixels that are considered to be in “smooth” regions, for some suitable definition of smooth.