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. The sixth and final part will compare these methods to standard deep learning techniques for denoising and filling in.

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 there are different levels. The top level (first row of the left figure) responds to high frequency information. The second row responds to medium frequency information while the bottom row responds to low frequency information.

I am going to follow the implementation details of this method based on Professor Rob Fergus' class assignment [7]. Now lets only focus on the first orientation (first column, all three rows). This first orientation will detect vertical edges. In the figure below we see three columns. The first column shows the filter again for reference. The second column shows the filters' outputs for the original uncorrupted image. The third column shows the filters' outputs for the image corrupted with gaussian noise of $\sigma =0.075$. The important thing to notice is that most of the noise shows up in the high frequency filtered output. In the two bottom rows for the original and noisy version, the filtered outputs look roughly the same. However, for the top row– the high frequency filtered output– the noisy version clearly looks a lot more corrupted than the original version.

Now let's focus only on the highest frequency filter, the first row from the above figure. We can plot out the distribution of the filtered outputs for both the original and noisy version of the image as you can see in the figure below. The top row is the same as the top row from figure 5. The bottom row shows the distributions as histograms. Notably the original filtered output's histogram has a far more sharply peaked distribution at the center than the noisy filtered output's histogram. The main intuition for this method of denoising is to try to transform the noisy filtered output's histogram to match the original filtered output's histogram.

Gradient descent + Langevin sampling

Coming next: interpreting denoising as sampling from the posterior with injected noise.

Filling in images with missing pixels

Coming next: the same Bayesian framework for inpainting when many pixels are missing.

Comparison to modern deep learning approaches

Coming next: qualitative + quantitative comparison against standard deep denoisers and inpainters.

Discussion

Takeaways, limitations, and next steps.

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