Implementation of Fourier Domain Denoising with Hard Threshold

Question

I just tried the Fourier denoising method with a hard threshold and my code is as follow:

F = imnoise(phantom(128),'gaussian');
B = fft(F);
s = 0.3;
for i = 1:128
for j = 1:128
    if abs(B(i,j))<s
        B(i,j)=0;
    end
end
end

C=ifft(B);
subplot(121)
imshow(F,[])
subplot(122)
imshow(C,[])

Left is the phantom after adding noise and the right is the denoised image from my code.

The Fourier denoising hard threshold method just uses threshold value to keep high frequency coefficients and the coefficients below the threshold to be 0.

In my code, I just added gaussian noise with default mean 0 and variance 0.01 to my phantom image and then I set the threshold to be 3 times of the standard deviation of the noise, which is 0.3 (this is some good practical threshold I looked at from literature). But why the result is not good? It seems it is still very noisy. Even I change my threshold the result is not good. Why? And I want to my denoised image as close to the original phantom.

Thanks in advance!

Answer 1

You added noise with Variance of 0.01 in the spatial domain.
Now you need to calculate what's the variance of the Noise per each bin in the DFT.

It's not the same, there are gains in the DFT domain, otherwise no one would use it.

By the way, there are much faster methods to apply "Hard Threshold" than the loop in your program. Use logical indexing to do it faster.

Answer 2

Thresholding in the Fourier domain is an archaic method sometimes called spectral subtraction, used for background noise removal in speech. Bad results in image processing can be due to several factors:

1. Misinterpretation: you say "The Fourier denoising hard threshold method just uses threshold value to keep high frequency coefficients". Not quite. You keep "high amplitude coefficients", no matter low or high frequency.
2. Test data: the Shepp-Logan model is known to be quite specific. Bad results on it do not imply bad results on real data, yet:
3. Transform choice: Fourier is very good for close-to-stationnary data. Standard images with smooth parts inside contours are not well represented in the Fourier domain, at least not concentrated enough for acurate data compaction and reconstruction with few (high) coeffcients. Derivatives, or multiscale directional derivatives like wavelets can help better.
4. Implementation: Fourier implementation via FFT suffers boundary artifacts coefficients, not talking about the 1D implementation you seem to use (prefer fft2 here)
5. Coefficient selection: hard and scalar thresholding typically yields more important visual artifacts than "smoother" thresholding and vector or group thresholding.
6. Formulation: you can generally get better results when you do not use a "parameter-tuned algorithm", but prefer a variational formulation based on your priors, including noise statistics.

So you have plenty of room to improve your results. Two references for multivariate image denoising with automated threshold (for Gaussian noise):

• A nonlinear Stein based estimator for multichannel image denoising, 2008
• Wavelet transform for the denoising of multivariate images, 2010

Answer 3

You should use the 2D Fourier transform function fft2 instead of fft. When provided with a 2D array arguement, fft returns 1D FFT results for the columns.

In order to remove the noise you should apply a low pass filter, as opposed to zeroing coefficients with small magnitudes. By using the fftshift and meshgrid functions you can easily do this in frequency space.

Answer 4

In literature it is shown that direct thresholding does not perform well. Consider iterative thresholding. You can find some papers on that ( see basis pursuit, matching pursuit etc.,). You can also read this paper - http://arxiv.org/abs/1504.00976