2 of 139

Motivation

Number of possible 200 x 200 images (of 256, i.e. 8 bit intensity levels) = 256^40000 = 2^320000 = 10^110000.

This is several trillion times the number of atoms in the universe (10^90).

Only a tiny subset of these are plausible as natural images.

Ajit Rajwade

3 of 139

Ajit Rajwade

4 of 139

Why study statistics of natural images?

Useful for many computer vision or image processing applications
Image denoising, deblurring, filling of missing pixels in images (inpainting)
Image compression
Classification into image categories

Ajit Rajwade

5 of 139

Ajit Rajwade

Sample State of the art result: Gaussian Noise sigma = 15

6 of 139

Ajit Rajwade

http://www.cse.cuhk.edu.hk/leojia/projects/motion_deblurring/

Motion deblurring

7 of 139

Ajit Rajwade

Inpainting

8 of 139

Why study statistics of natural images?

Natural image statistics also help us understand the human visual system better.

Ajit Rajwade

Real-world signal (scene – static or dynamic)

Convolution with Blur Kernel + Noise

Eye (Retina)

Visual Cortex

Brain is solving an inverse problem!

9 of 139

What are these amazing statistical properties?

Power law
Distribution of the DCT coefficients or wavelet coefficients of images or image patches
Relationships between these coefficients
Many more!

Ajit Rajwade

10 of 139

Some background

We know that an image is a 2D array of intensity values stored at pixels.

But images can be conveniently represented in the frequency domain as well.

Ajit Rajwade

11 of 139

Discrete Fourier transform

Given a 2D discrete signal (image) f(x,y) of size W₁ by W₂, its DFT is given as:

Here the image is being represented as a linear combination of complex exponentials of different frequencies.

Ajit Rajwade

12 of 139

2D Discrete Cosine Transform

Ajit Rajwade

Here the image is being represented as a linear combination of the cosine functions of different frequencies.

13 of 139

How do the DCT bases look like? (2D-case)

Ajit Rajwade

The DCT transforms an 8×8 block of input values to a linear combination of these 64 patterns. The patterns are referred to as the two-dimensional DCT basis vectors, and the output values are referred to as transform coefficients. Here each basis vector is reshaped to form an image.

Note: An image patch (size 8 x 8) can be represented as the linear combination of these 64 patterns.

http://en.wikipedia.org/wiki/JPEG

14 of 139

(1) Power law for natural images

The squared-amplitudes of the frequency components of an average natural image undergo rapid decay w.r.t. frequency pair magnitude f

where is a small number between 0 and 1 (usually around 0.19 for natural images), A is a constant.

This power law holds true across a range of scales (resolution) of the image.

Ajit Rajwade

Expected value of squared magnitude of Fourier transform coefficient at frequency (u,v). Expectation is over all natural images.

15 of 139

Ajit Rajwade

Spectrum of Barbara

log(abs(fftshift(fft2(im)+1)

Spectrum of Noise

log(abs(fftshift(fft2(im)+1)

16 of 139

Ajit Rajwade

(1) Power Law

17 of 139

(1) Power Law

Power law is tremendously useful for (lossy) image compression.
Image energy concentrated in just first few Fourier coefficients.
Remaining coefficients are small and can be ignored (i.e. considered to be 0) without much loss of information.
This principle is used by the JPEG algorithm (DCT coefficients), usually at a patch (block) level (patch-size of 8 x 8, usually).

Ajit Rajwade

18 of 139

(2) Joint statistics:�pixel and an immediate neighbor

Ajit Rajwade

What do you think the corresponding plot for a pure noise image would look like?

19 of 139

(2) Joint statistics:�pixel and an immediate neighbor

Ajit Rajwade

What do you think the corresponding plot for a pure noise image would look like?

20 of 139

(2) Joint statistics of pixel and an immediate neighbor

The histogram shows a high degree of correlation between the values of a pixel and one of its immediate neighbors.

Ajit Rajwade

21 of 139

(2a) Joint statistics of a pixel and nearby pixels: Experiment

Suppose you extract M ~ 100,000 small-sized (8 x 8) patches from a set of images.
Compute the column-column and row-row correlation matrices.

The correlation values in the matrix decrease almost in geometric progression with respect to distance – for short distances.

Ajit Rajwade

Code: http://www.cse.iitb.ac.in/~ajitvr/CS663/DCT_PCA/

22 of 139

Ajit Rajwade

1.0000 0.9902 0.9795 0.9733 0.9682 0.9639 0.9604 0.9570

0.9902 1.0005 0.9908 0.9795 0.9734 0.9684 0.9643 0.9605

0.9795 0.9908 1.0010 0.9908 0.9796 0.9735 0.9689 0.9646

0.9733 0.9795 0.9908 1.0005 0.9904 0.9793 0.9735 0.9686

0.9682 0.9734 0.9796 0.9904 1.0004 0.9903 0.9794 0.9734

0.9639 0.9684 0.9735 0.9793 0.9903 1.0001 0.9903 0.9793

0.9604 0.9643 0.9689 0.9735 0.9794 0.9903 1.0004 0.9904

0.9570 0.9605 0.9646 0.9686 0.9734 0.9793 0.9904 1.0002

1.0000 0.9888 0.9770 0.9704 0.9648 0.9599 0.9554 0.9510

0.9888 1.0004 0.9891 0.9768 0.9703 0.9646 0.9596 0.9548

0.9770 0.9891 1.0004 0.9886 0.9764 0.9698 0.9640 0.9587

0.9704 0.9768 0.9886 0.9994 0.9878 0.9755 0.9687 0.9627

0.9648 0.9703 0.9764 0.9878 0.9986 0.9870 0.9746 0.9676

0.9599 0.9646 0.9698 0.9755 0.9870 0.9978 0.9861 0.9734

0.9554 0.9596 0.9640 0.9687 0.9746 0.9861 0.9967 0.9847

0.9510 0.9548 0.9587 0.9627 0.9676 0.9734 0.9847 0.9951

CR/CR(1,1) -

Notice it can be approximated by the form shown two slides before, with ρ~0.99

CC/CC(1,1) -

Notice it can be approximated by the form shown two slides before, with ρ~0.9888

23 of 139

(3) Distribution of DCT coefficients

Due to the JPEG standard, the DCT is widely used in image processing.

The DCT is performed on 8 x 8 blocks.

The distribution of the DCT coefficients of a fixed frequency value (u,v) across different blocks has a peculiar shape – see next slide.

Ajit Rajwade

24 of 139

Ajit Rajwade

DCT coefficients computed for small patches. The distribution is represented by means of a histogram per coefficient. The samples to build the distribution for each coefficient come from the image patches – for each of which the DCT was computed.

Image source: Lam and Goodman, “A Mathematical Analysis of the DCT coefficient distributions for images”, IEEE Transactions on Image Processing, 2000.

https://www.researchgate.net/publication/3327260_A_mathematical_analysis_of_the_DCT_coefficient_distributions_for_images

25 of 139

(3) Distribution of DCT coefficients

The shapes of these histograms can be approximated to a high degree of accuracy by means of some distributions with a precise parametric form.

In particular, a Laplacian distribution has been experimentally shown to be a very good fit for all except the DC coefficient (u = v = 0).

Ajit Rajwade

26 of 139

Segway: Generalized Gaussian Distribution

Ajit Rajwade

Shape

parameter

Scale

parameter

Gaussian

Laplacian

Uniform

Gaussian

Generalized Gaussian

27 of 139

Ajit Rajwade

Generalized Gaussian Distributions: The Laplacian is a GGD with beta = 1.

https://en.wikipedia.org/wiki/Generalized_normal_distribution

28 of 139

Laplacian versus Gaussian distribution

Ajit Rajwade

The probability density of a sample decreases as the sample value moves farther away from the mean. But for the Gaussian distribution, this decrease is much quicker (why?).

29 of 139

(3) Distribution of DCT coefficients

Is there an explanation as to why the DCT coefficients have a Laplacian distribution?

Consider the formula for a DCT coefficient:

Ajit Rajwade

By Lindeberg’s central limit theorem, the distribution of F(u,v) should be Gaussian – the weighted summation of identically distributed random variables can be approximated as Gaussian.

30 of 139

(3) Distribution of DCT coefficients: Lindeberg’s central limit theorem

Consider independent random variables X_k, 1≤ k ≤ n – with mean μ_k and std. dev. σ_k respectively. Here n = N x N.
Define .

If the following condition holds, then the distribution of the random variable Z_n converges to N(0,1) as n→∞:

Ajit Rajwade

https://en.wikipedia.org/wiki/Lindeberg%27s_condition

31 of 139

(3) Distribution of DCT coefficients: Lindeberg’s central limit theorem

If we consider all pixels in a N x N (n=N²) patch to be random variables with variance σ², then we have:

Clearly the sufficient condition for Lindeberg’s CLT is satisfied in our scenario! This is because the variance of no one r.v. dominates over the sum s_n.

Ajit Rajwade

Spatial index k corresponding to location (x,y)

Random variable X_k

32 of 139

(3) Distribution of DCT coefficients

Note that the distribution is Gaussian even if the pixel values in a patch are spatially correlated – as long as the correlation coefficient is less than 1.
A small number of elements in the weighted summation is enough for the Gaussian distribution via the central limit theorem.
But note that this is under the assumption that the pixel values are identically distributed.
This assumption is true within a patch but not true across patches – as the variance of the intensity values in different patches can be very different!

Ajit Rajwade

33 of 139

(3) Distribution of DCT coefficients

In terms of mathematical equations:

What can we say about p(σ²)?

Ajit Rajwade

If σ² = variance of pixel intensity values, then

Var(F(u,v)) = σ²

34 of 139

Ajit Rajwade

Histogram of patch variance values – collected over all non-overlapping patches of size 8 x 8 from a grayscale image

35 of 139

(3) Distribution of DCT coefficients

What can we say about p(σ²)?
It can be modelled quite closely by an exponential distribution.

Ajit Rajwade

https://en.wikipedia.org/wiki/Exponential_distribution

36 of 139

(3) Distribution of DCT coefficients

So let us complete the math!
We have

Ajit Rajwade

This is a Laplacian distribution!

37 of 139

(3) Distribution of DCT coefficients

The widths of the Laplacian distribution of the coefficients however decrease with increase in frequency – go back a few slides.
Why so?
Remember that each DCT coefficient is a weighted summation of the intensity values.
The weights are of the following form:

Ajit Rajwade

38 of 139

(3) Distribution of DCT coefficients

The weights are of the following form:

At lower frequencies, the weights applied to adjacent pixels are close in value, and this increases the magnitude of the summation.

At higher frequencies, the weights applied to adjacent pixels are of similar magnitude but opposite sign, due to which spatially correlated values nullify each other giving smaller-valued summations. Remember: most images are piecewise smooth!

Hence the average value of the high frequency coefficients is often small in magnitude, and the average value of low frequency coefficients is large in magnitude.

This also affects the variances and second order uncentralized moments.

Ajit Rajwade

39 of 139

(3) Distribution of DCT coefficients

The DCT coefficients of small patches of images are quite sparse – very few of them have significant value and the rest have value close to 0.

The distribution is heavy-tailed, i.e. the probability of values significantly larger than the mean does not go to zero as quickly as a Gaussian distribution would predict!

Ajit Rajwade

40 of 139

(3) Distribution of DCT coefficients

If we assumed a different prior for p(σ²), note that the final distribution will be different from a Laplacian.

In general, the DCT coefficients of natural images can be modelled as a Generalized Gaussian distribution with shape parameter less than or equal to 1.

Ajit Rajwade

41 of 139

(4) Distribution of Wavelet coefficients

We will study first a very simple form of the wavelet transform – called the Haar wavelet.

The Haar wavelet is a sequence of rescaled square-shaped functions which together form an orthonormal basis.

The wavelet transform basically involves expressing the image as a linear combination of Haar wavelet basis functions.

Ajit Rajwade

42 of 139

(4) Distribution of Wavelet coefficients

The basic Haar wavelet functions are shown below:

Ajit Rajwade

Corresponding orthonormal basis (each column is a basis vector)

Image source: Huang & Mumford, Statistics of Natural Images and Models, CVPR 1999

http://www.dam.brown.edu/ptg/MDbook/Huangthesis.pdf

43 of 139

(4) Distribution of Wavelet coefficients

The wavelet transform is computed in a multi-scale manner.

Imagine you had a 64 x 64 image.

Take each 2 x 2 patch (non-overlapping) and compute the four wavelet coefficients – low-pass, vertical, horizontal and diagonal.

Ajit Rajwade

44 of 139

(4) Distribution of Wavelet coefficients

Each set of coefficients can be organized to form a 32 x 32 image – called a sub-band image.

This is called as a level 1 wavelet decomposition.

The low-pass sub-band can then be subjected to a second-level wavelet decomposition to generate 16 x 16 sub-band images.

The level two 16 x 16 low-pass sub-band can then be subjected to a level 3 wavelet decomposition, and so on till a maximum level 6 (yielding a single coefficient).

Ajit Rajwade

45 of 139

(4) Distribution of Wavelet coefficients

A sample level-two Haar wavelet decomposition is shown below:

Ajit Rajwade

46 of 139

(4) Distribution of Wavelet coefficients

Ajit Rajwade

Smoother regions: larger coefficient values

Textured regions/edges: smaller values

Image source: Simoncelli, Bayesian denoising of visual images in the wavelet domain, 1999, http://www.cns.nyu.edu/pub/lcv/simoncelli98e.pdf

47 of 139

(4) Distribution of Wavelet coefficients

Follow exponential decay rule (like Fourier coefficients) – small coefficients can be ignored, the remaining can be coded.

Ajit Rajwade

48 of 139

Ajit Rajwade

512 x 512 Barbara image

Image reconstructed from top 80,000 largest DCT coefficients

49 of 139

(4) Distribution of Wavelet coefficients

The significant wavelet coefficients can be (say) Huffman encoded using their histogram.

This can be used in image compression algorithms.

But you can do even better!

Ajit Rajwade

50 of 139

�(4a) Joint Statistics of Haar wavelet coefficients

Ajit Rajwade

Concept of parent, child, sibling and cousin coefficients (all are called wavelet sub-bands). Sibling = adjacent spatial locations in a sub-band, cousins = same spatial location at adjacent orientations.

Coefficients computed at multiples scales of the Haar wavelet pyramid

Image source: Huang & Mumford, Statistics of Natural Images and Models, CVPR 1999

http://www.dam.brown.edu/ptg/MDbook/Huangthesis.pdf

51 of 139

Ajit Rajwade

HL1

HH1

LL1

LH1

HH2

LH2

LL2

HL2

Three-level wavelet decomposition of an image

LL3

HL3

LH3

HH3

52 of 139

Ajit Rajwade

These joint statistics reveal that wavelet coefficients are NOT independent

Image source: Huang & Mumford, Statistics of Natural Images and Models, CVPR 1999

http://www.dam.brown.edu/ptg/MDbook/Huangthesis.pdf

53 of 139

Ajit Rajwade

Large magnitude coefficients tend to occur at neighboring spatial locations within a sub-band, or at the same locations in sub-bands of adjacent scale/orientation, or at related spatial locations in parent and child sub-bands.

Image source: Buccigrossi et al, Image Compression via Joint Statistical

Characterization in the Wavelet Domain.

http://www.cns.nyu.edu/pub/lcv/buccigrossi97a.pdf

54 of 139

Ajit Rajwade

Joint histogram of logarithms of absolute value of child and parent coefficients (also true for coefficients at adjacent orientations) – features like prominent edges or textures have large magnitude coefficients at multiple scales

Joint histogram of logarithm of squared value of child coefficient, and logarithm of square of linear combination of squared values of coefficients from neighboring sub-bands

55 of 139

Ajit Rajwade

HL1

HH1

LL1

LH1

HH2

LH2

LL2

HL2

LL3

HL3

LH3

HH3

Wavelet coefficient from neighboring sub-band (at corresponding location)

56 of 139

Ajit Rajwade

Conditional distributions of a wavelet coefficient (log absolute) given linear combinations of its neighbors (log absolute) – the shapes are robust across images! The plots are mean and variance normalized.

This redundancy means that we need not store all the wavelet coefficients – we can just predict some coefficients directly given their neighbors using the linear model that was fit. Useful for lossy image compression!

57 of 139

Why study Natural Image Statistics: Bayesian Framework

Ajit Rajwade

Observation

Unknown (to be determined) signal

noise

Prior Model on signal

Likelihood

Bayes rule

Known

operator

Posterior

probability

58 of 139

Why study Natural Image Statistics: Bayesian Framework

Ajit Rajwade

Assuming zero-mean i.i.d. (additive) Gaussian noise model,

the likelihood is as follows:

NOTE: Likelihood derived from the assumed noise model

59 of 139

Bayesian Framework: To Estimate x

Maximum a posteriori (MAP) estimate:

Ajit Rajwade

As y does not affect maximization w.r.t. x

The MAP estimate asks the following question: Given the observation y, what x is the most likely, taking into account that we have prior information on x in the form of p(x)?

If p(x) were a uniform distribution (or effectively we had no prior information about x), then MAP reduces to maximizing p(y|x) – which is called the maximum likelihood estimate.

60 of 139

Bayesian Framework: To Estimate x

Minimum mean square error (MMSE) estimate:

Ajit Rajwade

Prior is important!

Integrate to 1

61 of 139

Simple Example: 1

Ajit Rajwade

Observed Value of y = 14. Determine x given y and the knowledge of the noise model (likelihood) and prior on x.

Prior

Likelihood

62 of 139

Simple Example :1

Ajit Rajwade

63 of 139

Simple Example: 2

Ajit Rajwade

Observed Value of y = 2. Determine x given y and the knowledge of the noise model (likelihood) and prior on x.

64 of 139

Simple Example: 2

Ajit Rajwade

65 of 139

Simple Example: 2

Ajit Rajwade

When the likelihood and prior are both Gaussian, the MAP and MMSE estimates are equal.

66 of 139

Maximum Likelihood Estimation

Ajit Rajwade

If there is only one observation sample available (assume Gaussian noise), what is the maximum likelihood estimate of x?

If there are some N observation samples available (under Gaussian noise), what is the maximum likelihood estimate of x?

67 of 139

Ajit Rajwade

These two examples were very simple and involved scalar quantities. In future lectures, we will use more complex examples, where the unknown quantity x will be multivariate – in fact, it will be an image.

We will study the applications of natural image statistics in the following applications:

Image denoising and deblurring
Scene categorization
Image denoising (another flavour)
Reflection Removal

68 of 139

Application in Denoising or Deblurring

Ajit Rajwade

69 of 139

Application in Denoising

Consider the following noise model:

Given y, and knowing σ, determine the underlying image x.

Exploit the prior (fact) that the image x has DCT coefficients which are Laplacian distributed.

Ajit Rajwade

70 of 139

Application in Denoising

Let the DCT coefficients be given as follows:

So the estimation problem is

Ajit Rajwade

71 of 139

Application in Deblurring

Consider the following noise model:

Given y, and knowing H and σ, determine the underlying image x.

Exploit the prior (fact) that the image x has DCT coefficients which are Laplacian distributed.

Ajit Rajwade

72 of 139

Application in Deblurring

Let the DCT coefficients be given as follows:

So, the estimation problem is

Ajit Rajwade

73 of 139

Application in deblurring

The convolution kernel h is assumed to be known here – we will not immediately deal with the case where it is unknown!

The convolution of x with h can be equivalently represented by the product of a circulant matrix H derived from h, with the signal vector x.

Ajit Rajwade

74 of 139

Application in deblurring

A circulant matrix is a matrix where each row is a right circular shift of its preceding row in the following form:

Ajit Rajwade

75 of 139

Gaussian instead of Laplacian prior

What would happen if you imposed a Gaussian prior on the DCT coefficients?

Ajit Rajwade

76 of 139

Gaussian instead of Laplacian prior

Taking derivative w.r.t. θ, we get:

However for natural images or image patches, the Laplacian prior on the DCT or wavelet coefficients, is better suited – and yields better results in denoising.

Ajit Rajwade

This is the Wiener filter which we have seen last semester! The Wiener filter is the optimal linear filter regardless of the signal prior, which is what we proved in CS 663. For Gaussian likelihood and Gaussian prior, the Wiener filter is the optimal filter (among linear as well as non-linear) in a MAP or MMSE sense.

77 of 139

Ajit Rajwade

i.e. solution with Laplacian prior

78 of 139

Limitation of this model

For some images, a GGD with shape parameter less than 1 is more suitable to model the DCT coefficients than a Laplacian.
In such cases, the optimization problem becomes:

The problem however is that this is a non-convex optimization problem – and hence the local minima are different from the global minimum.
With Laplacian or Gaussian priors, the problems were convex. Many convex problems have efficient solutions.

Ajit Rajwade

79 of 139

Limitation of this model

As pointed out in class, the DC coefficient is often not well modeled by the Laplacian distribution.

The model also assumes that all DCT coefficients (or wavelet coefficients) are statistically independent – which is not true always as we shall see later.

Ajit Rajwade

80 of 139

Statistical Compressed Sensing

This is another view of compressed sensing based on Bayesian statistics.
Consider compressive measurements of the form:

Suppose .

Ajit Rajwade

81 of 139

Statistical Compressed Sensing

Consider the MAP solution for x given y and Φ:

Ajit Rajwade

82 of 139

Statistical Compressed Sensing

Consider the MAP solution for x given y and Φ:

Ajit Rajwade

The latter expression follows by the Woodbury matrix identity.

https://en.wikipedia.org/wiki/Woodbury_matrix_identity

83 of 139

Statistical Compressed Sensing

Since the likelihood and prior are both Gaussian, this is also the MMSE solution.

Unlike traditional compressed sensing, this solution is in closed form!

Hence better reconstruction speed!

Ajit Rajwade

84 of 139

Results

Ajit Rajwade

k = m = # of measurements

Assumption:

Eigen-values in the covariance matrix for the signal (i.e. Σ) are of the form: i^-αwhere 1 ≤ i ≤ n. The larger the value of α, the lower the reconstruction error.

The decay of the eigenvalues of the covariance matrix is the equivalent of signal sparsity or compressibility in an appropriate orthonormal basis.

For example, a signal with identity covariance matrix would not yield itself to good reconstruction via this method.

https://arxiv.org/pdf/1101.5785.pdf

Best-k = best k-term approximation obtained by keeping the k largest (absolute value) entries of x and setting the rest t0 0. This is for x ~ N(0,C) where C is a covariance matrix with decaying eigenvalues.

85 of 139

Gaussian assumption on signals?

Statistical compressed sensing has been applied to patch-based compressed sensing (where have we encountered this?)

Is the Gaussian assumption valid on patches?

Not a single Gaussian, but a mixture of Gaussians – also called a Gaussian mixture model (GMM).

Ajit Rajwade

86 of 139

Gaussian assumption on signals?

The probability density function of an image patch is expressed as a GMM, i.e. as follows:

In fact, GMMs are known to be universal PDF estimators.

That is any smooth PDF can be approximated to an arbitrary degree of accuracy using an appropriate number of Gaussians.

Ajit Rajwade

87 of 139

Gaussian assumption on signals?

Choice of K is a bit of an art, though there are techniques for it (eg: cross-validation).

Too large a K leads to overfitting, too small a value of K leads to underfitting.

A GMM could be fit to a bag of some N small-sized patches using an algorithm called Expectation Maximization (EM).

You will study EM in a machine learning or non-parametric statistics course.

Ajit Rajwade

88 of 139

GMMs for compressed sensing

Perform the GMM fitting for a large bag of patches.

Consider compressive measurements y of a patch x as follows:

Ajit Rajwade

89 of 139

GMMs for compressed sensing

From y, we find component-wise MAP estimates (for each of the k different mixture components):

Which of these estimates do we select? The j-th one as per the following:

Ajit Rajwade

90 of 139

GMMs for compressed sensing

This is called the MAP step.

In some variants, the K mixture components (means, covariance matrices and mixing values) are re-estimated from the signal reconstructions in the previous step.

This step is called the ML step (maximum likelihood).

These steps are repeated till convergence, i.e. till when the signal estimates do not change much.

This is called the MAP-EM or the MAP-ML or the Max-Max algorithm.

Ajit Rajwade

91 of 139

Ajit Rajwade

Yu and Sapiro, Statistical Compressed Sensing of Gaussian Mixture Models, IEEE TSP 2011

92 of 139

Ajit Rajwade

Yu and Sapiro, Statistical Compressed Sensing of Gaussian Mixture Models, IEEE TSP 2011

93 of 139

Application in Scene Categorization

Ajit Rajwade

94 of 139

Problem statement

Consider a set of different scene categories.

Given an image belonging to one of those categories, classify to which one it belongs.

Next slide shows average images of different categories – with the main object (if any) chosen at a fixed scale and images translated such that the main object is in the center.

Ajit Rajwade

95 of 139

Ajit Rajwade

Image source: Torralba and Olvia, Statistics of Natural Image Categories, Network, 2003

96 of 139

Man-made versus natural scenes

We know the power law for an average image:

But across image categories, the plots of the Fourier spectrum vary considerably even though they obey the power law.

In man-made scenes, the horizontal and vertical edges are typically dominant – in natural scenes, they aren’t.

Ajit Rajwade

97 of 139

Ajit Rajwade

50% contour

80% contour

X% contour means that the energy inside the contour is X%. Energy refers to sum of the squares of the amplitudes of the Fourier components inside the contour. How is the contour computed? The Fourier coefficients are sorted in descending order of squared magnitude and the first set of coefficients amounting to X% of energy are collected together to create this contour.

Image source: Torralba and Olvia, Statistics of Natural Image Categories, Network, 2003

Natural scene (more “isotropic” than the signature for man-made scenes)

Man-made scene (notice the larger energy concentrated near the axes)

98 of 139

Ajit Rajwade

Domination by the horizon (strong horizontal edges)

More isotropic spectral signatures!

99 of 139

Ajit Rajwade

Average spectral signatures for scenes at different scales.

Close-ups of man-made scenes are dominated by smooth surfaces, as opposed to close-ups of natural scenes. Smooth surfaces = dominance of low frequency information! Hence the spectra of man-made scenes have a dominance of low-frequency content – as is evident in the plots above.

Image source: Torralba and Olvia, Statistics of Natural Image Categories, Network, 2003

100 of 139

Ajit Rajwade

100

Scene scale: average scene depth

Image scale: related to its resolution (different scales generated by upsampling the image by a factor of 2 successively)

Image source: Torralba and Olvia, Statistics of Natural Image Categories, Network, 2003

101 of 139

Scene classification

So how can we apply the observations about scene spectra to scene classification?

We can apply PCA or LDA on such signatures.

Ajit Rajwade

101

102 of 139

Application of Dependencies between wavelet coefficients: in denoising

Ajit Rajwade

102

103 of 139

Recall 1

Wavelet coefficients of images are Laplacian distributed!

The various wavelet coefficients are not statistically independent.

Ajit Rajwade

103

Large wavelet coefficients tend to occur near each within the same sub-band.

And at the same relative spatial locations in sub-bands at adjacent scales or orientations

104 of 139

Wavelet coefficient dependency

Ajit Rajwade

104

Image source: Buccigrossi et al, Image Compression via Joint Statistical

Characterization in the Wavelet Domain, IEEE Transactions on Image Processing, 1997

105 of 139

Wavelet coefficient dependency

The conditional density of the child wavelet coefficient (c) given the parent (p) (figure 6A,B two slides before) reveals:
E(c|p) = 0 for all values of p.
They are not independent statistically – because the variance of c depends on the value of p.
The right side of the conditional density of the log of the squared coefficient is unimodal and concentrated on a unit slope line. In fact, E(c²|p²) is proportional to p².
Left side shows c being constant (not dependent on p).
This pattern is also observed for siblings (adjacent spatial locations), cousins (same spatial location, adjacent orientations).
This pattern is robust across a wide range of images.

Ajit Rajwade

105

106 of 139

Wavelet coefficient dependency

So how do we model this mathematically?

Here is one statistical model:

Ajit Rajwade

106

Neighbors of coefficient c

(some cousins & siblings, parent)

Can be obtained by least squares method

107 of 139

Ajit Rajwade

107

Can be obtained by least squares method

Least squares estimate of w

108 of 139

Wavelet coefficient dependency

Ajit Rajwade

108

109 of 139

Application to denoising

We have already seen the formula:

The same formula is applicable here with the following modification:

Ajit Rajwade

109

110 of 139

Application to denoising

But this is a chicken and egg problem – because we do not know the values of {w_k} or α beforehand!
But we can estimate these values by minimizing:

The values of {p_k} are obtained from a denoising algorithm that ignores wavelet coefficient dependency, e.g. using

and then used for estimating {w_k} or α.

Ajit Rajwade

110

111 of 139

Sample results

Ajit Rajwade

111

(Left) Original and (Right) noisy image

(Left) Marginal MAP with independent Gaussian prior and (Right) new model

MMSE estimator using GGD prior on wavelet coefficients

112 of 139

User-assisted reflection removal

Ajit Rajwade

112

113 of 139

Problem statement

A picture of the scene outside taken through a glass window contains undesirable reflections – of the scene on the side of the photographer.

The reflection can be weakened using a polarizer in front of the camera lens but that is unavailable in all cameras and it is unwieldy to use.

Ajit Rajwade

113

114 of 139

Ajit Rajwade

114

115 of 139

Ajit Rajwade

115

This is an ill-posed problem, as there are infinitely many I₁and I₂ that could sum up to I.

116 of 139

Method for reflection removal: user-assisted

Allow a user to mark a set of points S₁ which the user thinks belong to I₁.

And another set of points S₂which the user thinks belong to I₂.

Intuitively, these points will belong to strong edges in the two images.

Ajit Rajwade

116

117 of 139

Ajit Rajwade

117

118 of 139

Method for reflection removal: user-assisted

Find I₁ and I₂ such that:
I₁ and I₂ sum up to I
The gradient of I₁ at points in S₁ should match the gradient of I at those points.
The gradient of I₂ at points in S₂ should match the gradient of I at those points.

Ajit Rajwade

118

119 of 139

Method for reflection removal: statistical model

Exploit a statistical property of a natural image.

The gradients are sparse!

Ajit Rajwade

119

120 of 139

Ajit Rajwade

120

The gradients (in this case, the y derivatives of the intensity) are sparse!

121 of 139

Method for reflection removal: statistical model

One possible statistical model for the gradients of the image is the following:

But the gradients can be computed in x and y directions, and we can have first as well as second order derivatives.
Each choice of derivative will be called as a gradient filter – of which we will have many.

Ajit Rajwade

121

z = gradient value,

q = shape parameter

122 of 139

Method for reflection removal: statistical model

Note that each gradient filter is applied to every pixel location of image I.
The output of filter k at pixel location i on image I₁ will be denoted as f_ik,1.
We will assume statistical independence of all these filter outputs.
So the overall model is:

Ajit Rajwade

122

123 of 139

Method for reflection removal: statistical model

Recall that we have I = I₁+ I₂ where I₁ and I₂ are both unknown.
Gradients of both images are sparse and independent of each other.
So we have to maximize wr.t. I₁ and I₂:

Equivalent to minimizing w.r.t. I₁ and I₂:

Ajit Rajwade

123

124 of 139

Method for reflection removal: statistical model

But we also have an important constraint that I₁ and I₂ sum up to I.
So the objective function becomes:

But that’s not enough. We need to enforce that the gradients of the given image I and the estimated I₁ agree at all points in set S₁.
Likewise for I₂ and S₂.

Ajit Rajwade

124

125 of 139

Method for reflection removal: statistical model

So the objective function becomes:

Ajit Rajwade

125

Because I₂ = I –I₁.

126 of 139

Optimization algorithm

Given the statistical model for the gradient filter outputs, the function ρ is non-convex if q < 1 (it is convex if q = 1).

The optimization procedure for this is not very easy.

The authors use a method called iteratively reweighted least squares (IRLS).

Ajit Rajwade

126

127 of 139

Segway: Least squares method

Consider the solution to the following problem:

This is a least squares problem, and it has a well-known pseudo-inverse based solution.

Now we will look at some flavours of least squares.

Ajit Rajwade

127

128 of 139

Segway: Weighted least squares method

Now consider the solution to the following problem:

Here W is a n x n diagonal matrix containing weights which give different levels of importance to each entry of y.
The solution of this is again in terms of a pseudo-inverse.

Ajit Rajwade

128

129 of 139

Segway: Least p-norm problem

Consider the solution to the following problem:

This has no known closed form solution!

Instead an iterative procedure has been proposed – called IRLS.

Ajit Rajwade

129

130 of 139

Segway: IRLS

The IRLS at step t+1 involves a weighted least squares problem:

At t = 1, the weights are set to 1.
The weights are updated as follows:

This is done till convergence.

Ajit Rajwade

130

Weight for point i at iteration t

Diagonal matrix of weights at iteration t (for all points)

131 of 139

Segway: IRLS

The weights are updated as follows:

Why these weights? Simply because the problem can be re-written as follows:

Ajit Rajwade

131

132 of 139

Back to the Optimization algorithm

The objective function is:

It can be expressed as:

Ajit Rajwade

132

133 of 139

Optimization algorithm

Ajit Rajwade

133

134 of 139

Method for reflection removal: actual statistical model used in the paper

The statistical model for the gradients of the image is chosen to be the following:

This is a mixture of two Laplacian distributions and it is seen to be sparser than a single Laplacian.

Ajit Rajwade

134

z = gradient value