Denoising Images with Natural Scene Statistics

Introduction

This blog will introduce and demonstrate how to use natural image statistics in order to denoise an image. The first part will explore the 1/f² power spectrum of natural images. The second part will discuss using bayesian inference in order to denoise an image. The third part will implement Simoncelli et al.’s “NOISE REMOVAL VIA BAYESIAN WAVELET CORING” in order to do denoising. The fourth part will discuss using gradient descent with langevin sampling in order to denoise an image. The fifth part will demonstrate filling in an image where a large percentage of pixels are missing.

What do we know about natural images? (1/f² power spectrum)

Click here to expand Fourier transform explanation Click here to hide Fourier transform explanation

Fourier's theorem was first proposed in 1807 but not published until 1822 due to the skepticism of many leading mathematicians who did not think it could be true. Today, it provides the theoretical foundation for the Fourier transform which is one of the most widely used tools in all of science and engineering, from radio astronomy and signal processing to neuroscience and modern machine learning. For example, it is used for:

• analyzing the frequency content of natural signals such as sound waveforms.
• detecting rhythmic activity in brain signals - e.g., EEG, LFP and spike train recordings.
• even when there is no periodic structure - characterizing the length scale of correlations over time or space, since the power spectrum is the Fourier transform of the auto-correlation function.
• identifying the input-output relationship of linear systems.

Here we will experiment with Fourier decompositions of simple 1D waveforms and 2D images, to gain some intuition about what the Fourier transform tells you about a signal, and how it can be used to manipulate images to demonstrate perceptual phenomena.

Fourier's theorem states that any continuous function $f(x)$ can be described as a weighted sum of sine and cosine waveforms at different frequencies:

$$f(x)=\sum _{\omega =0}^{\infty }{a}_s(\omega )\sin (\omega x)+a_c(\omega )\cos (\omega x)$$

where $a_s(\omega )$ and $a_c(\omega )$ are the weights of the sine and cosine components, respectively, and $\omega$ is the frequency in $2\pi$ radians per unit of $x$ (space or time). Alternatively if we define

$$A(\omega )\equiv \sqrt{{a_c(\omega )}^2+{a_s(\omega )}^2}$$ $$\varphi (\omega )\equiv {\tan }^{-1}\frac{a_c(\omega )}{a_s(\omega )}$$

we may rewrite the equation equivalently as

$$f(x)=\sum _{\omega =0}^{\infty }A(\omega )\sin (\omega x+\varphi (\omega ))$$

where now $f(x)$ is described simply as a sum of cosine waves, each with amplitude $A(\omega )$ and phase-shifted by $\varphi (\omega )$. We refer to $A(\omega )$ as the amplitude spectrum and $\varphi (\omega )$ as the phase spectrum of $f(x)$.

To gain some intuition for the power and generality of the Fourier theorem, let us first decompose a simple step-function this way. The function $f(x)$ is

$$f(x)=\{\begin{array}{ll}1& \text{for}x\ge 0\\ -1& \text{for}x<0\end{array}$$

You are given $f(x)$ as a discretely sampled function. In the image below you can see that even though $f(x)$ looks nothing like a sinusoid, we can nevertheless describe it perfectly by adding sinusoids of different frequencies together, as long as they are given the correct amplitudes and phases.

The direct method of computing the amplitude and phase spectrum of a function $f(x)$ is the Fourier transform! Looking back at the Fourier theorem, we can compute the Fourier coefficients $a_s(\omega )$ and $a_c(\omega )$ from the discretely sampled function $\hat{f}\left[n\right]=f(n\Delta x),n=0:N-1$ , as follows:

$$a_s({\omega }_k)=\frac{2}{N}\sum _{n=0}^{N-1}\hat{f}\left[n\right]\sin ({\omega }_{k}n)\Delta x$$ $$a_c({\omega }_k)=\frac{2}{N}\sum _{n=0}^{N-1}\hat{f}\left[n\right]\cos ({\omega }_{k}n)\Delta x$$

where ${\omega }_k=2\pi k/N,k=0:N/2$. Intuitively, we can think of these equations as measuring the similarity between $\hat{f}\left[n\right]$ and the sine function, and the similarity between $\hat{f}\left[n\right]$ and the cosine function, for each frequency ${\omega }_k$.

Formally, we can show this gives us exactly the Fourier coefficients $a_s(\omega )$ and $a_c(\omega )$ in the Fourier theorem by plugging the expression for $f(x)$ into these equations and utilizing the orthogonality of sines and cosines of different frequencies, or sine and cosine of the same frequency.

Once you have computed $a_s(\omega )$ and $a_c(\omega )$, the amplitude spectrum and phase spectrum can be computed via $A(\omega )\equiv \sqrt{{a_c(\omega )}^2+{a_s(\omega )}^2}$ and $\varphi (\omega )\equiv {\tan }^{-1}(\frac{a_c(\omega )}{a_s(\omega )})$.

Click here to hide Fourier transform explanation

A 2D image may also be described in terms of a Fourier decomposition as follows:

where ${\omega }_x$ and ${\omega }_y$ correspond to frequencies in the $x$ and $y$ dimensions of the image, respectively. Each sinewave component now corresponds to a 2D grating of a specific spatial-frequency and orientation. You can use the demo below to build up an image from its Fourier components, one component at a time.

Note how the amplitudes fall off with frequency. At some point, the amplitudes are so small that the gratings are not even visible, yet they are clearly needed to reconstruct. This seems paradoxical - gratings which are invisible in isolation nevertheless create structure when combined. The squared amplitudes have a 1/f² [1] fall off relationship, which is shown in the dashed line. This phenomenon of 1/f power spectra is not unique to natural images — it occurs all throughout nature such as in music [2].

Bayesian denoising with a natural-image prior

If you want an incredible and easy-to-understand introduction to bayesian inference please read this write-up by Bruno Olshausen [3]. We will use bayesian inference in order to denoise images. Let’s begin with the following problem set-up (copied from the handout write-up [3]): “A classic example of where Bayesian inference is employed is in the problem of estimation, where we must guess the value of an underlying parameter from an observation that is corrupted by noise. Let’s say we have some quantity in the world, $x$, and our observation of this quantity, $y$, is corrupted by additive Gaussian noise, $n$, with zero mean:”

Our job is to make the best guess as to the value of $x$ given the observed value $y$. Let’s call our best guess $\hat{x}$. In this set-up, when $x\sim N({\mu }_x,{\sigma }_x^2)$ and $n\sim N({\mu }_n,{\sigma }_n^2)$ are both gaussian, there is a standard closed-form solution for $\hat{x}$. This is also known as the Wiener filter:

Now assume ${\mu }_x=0$ and ${\sigma }_x=1.0$. Then $\hat{x}$ simplifies to the following form:

$a$ is an attenuation factor between 0 and 1. $a$ varies with ${\sigma }_n$ as shown in the plot below:

Now let's apply the same framework to the problem of denoising an image. We are given a noisy observation $I(\overrightarrow{x})$ of an image $I_0(\overrightarrow{x})$ as follows:

where $\overrightarrow{x}$ denotes the pixel position within the 2D image. Our task is to estimate $I_0$ from $I$. We could do this pixel by pixel, but since the signal and noise variance is the same for each pixel all it will do is attenuate the entire image, which isn't very interesting. However if we transform to the frequency domain we can exploit our prior knowledge about the statistics of natural images - namely, they have 1/f² power spectra!

Letting $\tilde{I}(\overrightarrow{f})$ denote the Fourier transform of $I(\overrightarrow{x})$, we have

Since the Fourier transform is a unitary transformation (it is like multiplying by an orthonormal matrix), the noise variance in the frequency domain is the same as it was in the space domain. But the signal variance is now different for each frequency as it follows the 1/f² power law:

We just need to set the constant of proportionality $k$ so that the total signal power in the frequency domain equals the signal power in the space domain, which can be accomplished via

We can then denoise the image frequency by frequency as follows:

where $a(\overrightarrow{f})$ is the attenuation factor computed from the signal and noise variance for each frequency. In this case, then

Examine the code below which implements the above equations and be sure you understand it. Try running it for different amounts of noise to see how far you can go before the image is unrecoverable.

Wavelet Coring

In the previous section we used the Wiener filter to attenuate each frequency. In this section, instead of modifying frequency content, we will take wavelet decomposition of the image and modify the wavelet coefficients. This is a method developed by Simoncelli and Adelson in 1996 [4]. It relies on steerable pyramids [5], which is a method for decomposing an image. If you want a great explanation of steerable pyramids check out the actively maintained pyrtools documentation page about steerable pyramids [6]. Below is an example of the steerable pyramid decomposition. The left figure shows the filters and the right figure shows the activations on the Einstein image for each of the filters.

You can see from the figure above that the decomposition separates the image by both scale and orientation. The top row responds to the highest frequencies, the middle row responds to medium frequencies, and the bottom row responds to lower frequencies.

I am going to follow the implementation style from Professor Rob Fergus' computational photography assignment [7]. To keep the picture simple, focus on one orientation band: the first column below, which is most sensitive to vertical edges. The first column shows the filter, the second column shows the filter output for the clean Einstein image, and the third column shows the same output after adding Gaussian noise with $\sigma =0.075$. Notice that the highest-frequency band is much more corrupted by the noise than the lower-frequency bands.

Now focus only on the highest-frequency band. The clean image has a sharply peaked coefficient distribution: most coefficients are close to zero, but a few large coefficients correspond to real edges. Adding Gaussian noise spreads this distribution out. Wavelet coring uses this difference in distributions to decide which noisy coefficients should be suppressed and which should be kept.

The Bayesian version of coring starts with the same scalar observation model we used before, but now $x$ is a clean wavelet coefficient and $y$ is the noisy coefficient:

The least-squares optimal estimate is the posterior mean, $\hat{x}(y)$. By Bayes' rule this can be written as

Here $P_n$ is the noise distribution and $P_x$ is the prior distribution of clean coefficients. If both are Gaussian, this reduces to a linear Wiener-style shrinkage. But natural-image wavelet coefficients are not Gaussian: they are sharply peaked at zero with heavy tails. That changes the posterior mean into a nonlinear coring curve. Small observed coefficients are pulled toward zero because they are likely to be noise, while large coefficients are mostly preserved because they are more likely to be real image structure.

In the code, this posterior mean is estimated directly from histograms. The implementation cores the residual highpass band and the four finest oriented bands. For each band, the denominator is the convolution of the clean-coefficient histogram with the noise histogram, which gives the distribution of noisy observations. The numerator is the same convolution after weighting the clean coefficients by their value. Dividing the two gives a lookup table $\hat{x}(y)$ that is applied to every coefficient in that band. After coring the high-frequency bands, the steerable pyramid is reconstructed back into an image.

Gradient descent + Langevin sampling

The Wiener filter gave us the posterior mean directly. Another way to use the same posterior is to move an image estimate downhill on the negative log posterior. In the implementation, the prior term is still the 1/f² natural-image prior, so high frequencies are penalized more strongly than low frequencies.

The update is done in the Fourier domain. Let $\tilde{I}$ be the current estimate and $\tilde{Y}$ be the noisy observation. The code uses the gradient of the quadratic posterior energy, preconditions each frequency by its curvature, and optionally adds Langevin noise:

When $T=0$, this is just gradient descent toward the MAP estimate. When $T>0$, the injected noise turns the procedure into Langevin sampling, so the estimate jitters around plausible images under the posterior instead of collapsing to only one optimum.

Filling in images with missing pixels

The inpainting demo uses the same prior, but changes the likelihood. Some pixels are observed and must stay fixed; the rest are unknown. In the implementation, the mask keeps a small fraction of pixels (2% by default in the web demo) and shows the missing pixels in blue.

Since the known pixels are already constrained by the data, the update only changes the missing pixels. At each step, the code computes the gradient of the 1/f² prior energy by transforming the current image to the Fourier domain, multiplying by $\lambda (f)$, transforming back, and stepping downhill only where the mask is missing:

This is why a recognizable image can emerge from very few pixels: the observed pixels anchor the solution, while the 1/f² prior fills in the missing values with a power spectrum that looks like a natural image.

Discussion

The common thread through all of these examples is that natural images are highly structured. A simple 1/f² power-spectrum prior already gives useful denoising and inpainting behavior, and the Bayesian framing tells us exactly how to combine that prior with noisy or incomplete observations.

The Fourier-domain methods are deliberately simple: they capture the average falloff of power with frequency, but not the localized edge structure of images. Wavelet coring adds one important piece of that structure by using the sparse, heavy-tailed statistics of steerable-pyramid coefficients. This is why it can preserve edges while suppressing small noisy coefficients.

These methods are not meant to be the final word on image restoration. Their value is that they make the assumptions visible: choose a prior, write down the observation model, and then compute or sample from the posterior. That is the same basic logic behind many more powerful modern methods, even when the prior is learned rather than hand-specified.

@article{kothapalli2026denoising,
  author={Kothapalli, Tejasvi and Olshausen, Bruno},
  title={Denoising Images with Natural Scene Statistics},
  year={2026},
  url={https://tejasvikothapalli.github.io/blogs/denoising/}
}