Statistical Image Models Eero Simoncelli Howard Hughes Medical - - PowerPoint PPT Presentation

Statistical Image Models

Eero Simoncelli

Howard Hughes Medical Institute, Center for Neural Science, and Courant Institute of Mathematical Sciences New York University

Photographic Images

Diverse specialized structures:

• edges/lines/contours
• shadows/highlights
• smooth regions
• textured regions

Photographic Images

Diverse specialized structures:

• edges/lines/contours
• shadows/highlights
• smooth regions
• textured regions

Occupy a small region of the full space

space of all images t y p i c a l i m a g e s

One could describe this set as a deterministic manifold....

• Step edges are rare (lighting, junctions, texture, noise)

• Step edges are rare (lighting, junctions, texture, noise)
• One scale’s texture is another scale’s edge

• Step edges are rare (lighting, junctions, texture, noise)
• One scale’s texture is another scale’s edge
• Need seamless transitions from isolated features to

dense textures

One could describe this set as a deterministic manifold....

space of all images t y p i c a l i m a g e s

One could describe this set as a deterministic manifold.... But seems more natural to use probability

space of all images t y p i c a l i m a g e s

One could describe this set as a deterministic manifold.... But seems more natural to use probability

space of all images t y p i c a l i m a g e s

P(x)

“Applications”

• Engineering: compression, denoising, restoration,

enhancement/modification, synthesis, manipulation

[Hubel ‘95]

“Applications”

• Engineering: compression, denoising, restoration,

enhancement/modification, synthesis, manipulation

• Science: optimality principles for neurobiology (evolution,

development, learning, adaptation)

[Hubel ‘95]

Density models

nonparametric parametric/ constrained

Density models

nonparametric parametric/ constrained build a histogram from lots of

• bservations...

Density models

nonparametric parametric/ constrained build a histogram from lots of

• bservations...

use “natural constraints” (geometry/photometry

• f image formation,

computation, maxEnt)

Density models

nonparametric parametric/ constrained build a histogram from lots of

• bservations...

use “natural constraints” (geometry/photometry

• f image formation,

computation, maxEnt) historical trend (technology driven)

50 100 150 200 250

histogram Original image

Range: [0, 237] Dims: [256, 256]

50 100 150 200 250

histogram Original image

Range: [0, 237] Dims: [256, 256]

50 100 150 200 250

histogram Equalized image

Range: [1.99, 238] Dims: [256, 256]

50 100 150 200 250

histogram Original image

Range: [0, 237] Dims: [256, 256]

50 100 150 200 250

histogram Equalized image

Range: [1.99, 238] Dims: [256, 256]

General methodology

Observe “interesting” Joint Statistics Transform to Optimal Representation

General methodology

Observe “interesting” Joint Statistics Transform to Optimal Representation

General methodology

Observe “interesting” Joint Statistics Transform to Optimal Representation

“Onion peeling”

Evolution of image models

• I. (1950’s): Fourier + Gaussian
• II. (mid 80’s - late 90’s): Wavelets + kurtotic marginals
• III. (mid 90’s - present): Wavelets + local context
• local amplitude (contrast)
• local orientation
• IV. (last 5 years): Hierarchical models

I(x,y) I(x+1,y) I(x,y) I(x+2,y) I(x,y) I(x+4,y)

10 1 Spatial separation (pixels)

Correlation

a. b.

Pixel correlation

I(x,y) I(x+4,y)

10 20 30 40 1 Spatial separation (pixels)

Correlation

b.

I(x,y) I(x+1,y) I(x,y) I(x+2,y) I(x,y) I(x+4,y)

10 1 Spatial separation (pixels)

Correlation

a. b.

Pixel correlation

Translation invariance

Assuming translation invariance,

Translation invariance

Assuming translation invariance, => covariance matrix is Toeplitz (convolutional)

Translation invariance

Assuming translation invariance, => covariance matrix is Toeplitz (convolutional) => eigenvectors are sinusoids

Translation invariance

Assuming translation invariance, => covariance matrix is Toeplitz (convolutional) => eigenvectors are sinusoids => can diagonalize (decorrelate) with F.T.

Translation invariance

Assuming translation invariance, => covariance matrix is Toeplitz (convolutional) => eigenvectors are sinusoids => can diagonalize (decorrelate) with F.T. Power spectrum captures full covariance structure

Spectral power

Structural:

F(sω) = spF(ω) F(ω) ∝ 1 ωp

[Ritterman 52; DeRiugin 56; Field 87; Tolhurst 92; Ruderman/Bialek 94; ...]

Assume scale-invariance: then:

Spectral power

Structural:

F(sω) = spF(ω) F(ω) ∝ 1 ωp

[Ritterman 52; DeRiugin 56; Field 87; Tolhurst 92; Ruderman/Bialek 94; ...]

Assume scale-invariance: then:

1 2 3 1 2 3 4 5 6

Log10 spatialfrequency (cycles/image) Log10 power

Empirical:

Principal Components Analysis (PCA) + whitening

• 20

20

• 20

20 4

• 4

4

• 4

20

• 20

20

• 20

a. b. c.

PCA basis for image blocks

PCA basis for image blocks

PCA is not unique

Maximum entropy (maxEnt)

E (f(x)) = c f(x) = x2 f(x) = |x| The density with maximal entropy satisfying pME(x) ∝ exp (−λf(x)) is of the form Examples: where λ depends on c

Model I (Fourier/Gaussian)

Basis set: Image: Coefficient density:

F -1

1/f2 P(c)

Gaussian model is weak

P(x)

F−1

ω−2

F -1

1/f2 P(c)

Gaussian model is weak

a. b.

ω2 F F−1

P(x)

F−1

ω−2

F -1

1/f2 P(c)

Gaussian model is weak

a. b.

ω2 F F−1

• 20

20

• 20

20 4

• 4

4

• 4

20

• 20

20

• 20

a. b. c.

P(x)

F−1

ω−2

Bandpass Filter Responses

500 500 10

• 4

10

• 2

10

Filter Response Probability

Response histogram Gaussian density

[Burt&Adelson 82; Field 87; Mallat 89; Daugman 89, ...]

“Independent” Components Analysis (ICA)

For Linearly Transformed Factorial (LTF) sources: guaranteed independence

(with some minor caveats)

• 4

4

• 4

4 20

• 20

20

• 20

20

• 20

20

• 20
• 4

4

• 4

4

a. b. c. d.

[Comon 94; Cardoso 96; Bell/Sejnowski 97; ...]

ICA on image blocks

[Olshausen/Field ’96; Bell/Sejnowski ’97] [example obtained with FastICA, Hyvarinen]

Marginal densities

P(x) ∝ exp −|x/s|p

[Mallat 89; Simoncelli&Adelson 96; Moulin&Liu 99; ...]

Well-fit by a generalized Gaussian:

Wavelet coefficient value log(Probability) p = 0.46 H/H = 0.0031 Wavelet coefficient value log(Probability) p = 0.58 H/H = 0.0011 Wavelet coefficient value log(Probability) p = 0.48 H/H = 0.0014

Kurtosis vs. bandwidth

0.5 1 1.5 2 2.5 3 4 6 8 10 12 14 16 Filter Bandwidth (octaves) Sample Kurtosis

[after Field 87]

Note: Bandwidth matters much more than orientation

[see Bethge 06]

Octave-bandwidth representations

Filter: Spatial Frequency Selectivity:

Model II (LTF)

Basis set: Image: Coefficient density:

LTF also a weak model...

Sample Gaussianized

Sample ICA-transformed and Gaussianized

Trouble in paradise

Trouble in paradise

• Biology: Visual system uses a cascade
• Where’s the retina? The LGN?
• What happens after V1? Why don’t responses get

sparser? [Baddeley etal 97; Chechik etal 06]

Trouble in paradise

• Biology: Visual system uses a cascade
• Where’s the retina? The LGN?
• What happens after V1? Why don’t responses get

sparser? [Baddeley etal 97; Chechik etal 06]

• Statistics: Images don’t obey ICA source model
• Any bandpass filter gives sparse marginals [Baddeley 96]

=> Shallow optimum [Bethge 06; Lyu & Simoncelli 08]

• The responses of ICA filters are highly dependent

[Wegmann & Zetzsche 90, Simoncelli 97]

Conditional densities

• 40

• 40

• 40

40

• 40

40

[Simoncelli 97; Schwartz&Simoncelli 01]

[Schwartz&Simoncelli 01]

• Large-magnitude subband coefficients are found at

neighboring positions, orientations, and scales.

Method 1: Conditional Gaussian

[Simoncelli 97; Buccigrossi&Simoncelli 99; see also ARCH models in econometrics!]

Modeling heteroscedasticity

(i.e., variable variance)

P (xn|{xk}) ∼ N

• 0;
• k

wnk |xk|2 + σ2

Joint densities

adjacent near far

• ther scale
• ther ori

100 100 150 100 50 50 100 150 100 100 150 100 50 50 100 150 100 100 150 100 50 50 100 150 500 500 150 100 50 50 100 150 100 100 150 100 50 50 100 150 100 100 150 100 50 50 100 150 100 100 150 100 50 50 100 150 100 100 150 100 50 50 100 150 500 500 150 100 50 50 100 150 100 100 150 100 50 50 100 150

[Simoncelli, ‘97; Wainwright&Simoncelli, ‘99]

• Nearby: densities are approximately circular/elliptical
• Distant: densities are approximately factorial

ICA-transformed joint densities

d=2 d=16 d=32

kurtosis

• rientation

data (ICA’d): factorialized: sphericalized:

4 6 8 10 12 !/4 !/2 3!/4 ! 4 6 8 10 12 !/4 !/2 3!/4 ! 4 6 8 10 12 !/4 !/2 3!/4 !

ICA-transformed joint densities

d=2 d=16 d=32

kurtosis

• rientation

data (ICA’d): factorialized: sphericalized:

4 6 8 10 12 !/4 !/2 3!/4 ! 4 6 8 10 12 !/4 !/2 3!/4 ! 4 6 8 10 12 !/4 !/2 3!/4 !

• Local densities are elliptical (but non-Gaussian)
• Distant densities are factorial

[Wegmann&Zetzsche ‘90; Simoncelli ’97; + many recent models]

3 6 9 12 15 18 20 0.05 0.1 0.15 0.2

kurtosis

blk size = 3x3

blk spherical factorial

Spherical vs LTF

3x3 7x7 15x15

3 6 9 12 15 18 20 0.05 0.1 0.15 0.2

kurtosis

blk size = 7x7

blk spherical factorial 3 6 9 12 15 18 20 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

kurtosis

blk size = 11x11

blk spherical factorial

data (ICA’d): factorialized: sphericalized:

• Histograms, kurtosis of projections of image blocks onto random

unit-norm basis functions.

• These imply data are closer to spherical than factorial

[Lyu & Simoncelli 08]

• Zetzsche & Krieger, 1999;
• Huang & Mumford, 1999;
• Wainwright & Simoncelli, 2000;
• Hyvärinen and Hoyer, 2000;
• Parra et al., 2001;
• Srivastava et al., 2002;
• Sendur & Selesnick, 2002;
• Teh et al., 2003;
• Gehler and Welling, 2006
• Lyu & Simoncelli, 2008
• etc.

non-Gaussian elliptical observations and models of natural images:

• is Gaussian,
• and are independent
• is elliptically symmetric, with covariance
• marginals of are leptokurtotic

[Wainwright&Simoncelli 99]

Modeling heteroscedasticity

Method 2: Hidden scaling variable for each patch

Gaussian scale mixture (GSM)

[Andrews & Mallows 74]:

• u
• x = √z

u z

• x
• u

z > 0

• x

∝ Cu

• Empirically, z is approximately lognormal

[Portilla etal, icip-01]

• Alternatively, can use Jeffrey’s noninformative prior

[Figueiredo&Nowak, ‘01; Portilla etal, ‘03]

GSM - prior on z

pz(z) = exp (−(log z − µl)2/(2σ2

l ))

z(2πσ2

l )1/2

pz(z) ∝ 1/z

GSM simulation

!!" " !" #"

"

#"

!

!!" " !" #"

"

#"

!

Image data GSM simulation

[Wainwright & Simoncelli, NIPS*99]

Model III (GSM)

Basis set: Image: Coefficient density:

Original coefficients Normalized by √z marginal

500 500 10 8 6 4 2 Log probability 5 5 10 9 8 7 6 5 4 Log probability

joint

100 50 50 100 100 50 50 100 2 4 6 8 2 4 6 8

subband

[Schwartz&Simoncelli 01] [Ruderman&Bialek 94]

First Order Ideal Conditional Model 1 2 3 4 5 6 1 2 3 4 5 6 Empirical Conditional Entropy Model Encoding cost (bits/coeff) Gaussian Model Generalized Laplacian 3 4 5 2.5 3 3.5 4 4.5 5 5.5 Empirical First Order Entropy (bits/coeff) Model Encoding Cost (bits/coeff)

[Buccigrossi & Simoncelli 99]

Bayesian denoising

• Additive Gaussian noise:
• Bayes’ least squares solution is conditional mean:

y = x + w

P(y|x) ∝ exp[−(y − x)2/2σ2

w]

ˆ x(y) = I E(x|y)

=

• dxP(y|x)P(x)x/P(y)

• I. Classical

If signal is Gaussian, BLS estimator is linear:

denoised (ˆ x) noisy (y)

ˆ x(y) = σ2

x

σ2

x + σ2 n

· y

=> suppress fine scales, retain coarse scales

Non-Gaussian coefficients

[Burt&Adelson ‘81; Field ‘87; Mallat ‘89; Daugman ‘89; etc]

!!"" " !"" #"

!$

#"

!%

#"

"

&'()*+,-*./01.* 2+0343'(')5

• *./01.*,6'.)07+48

94:..'41,;*1.')5,,

• II. BLS for non-Gaussian prior
• Assume marginal distribution [Mallat ‘89]:
• Then Bayes estimator is generally nonlinear:

P(x) ∝ exp −|x/s|p

p = 2.0 p = 1.0 p = 0.5

[Simoncelli & Adelson, ‘96]

MAP shrinkage

p=2.0 p=1.0 p=0.5

[Simoncelli 99]

Denoising: Joint

I E(x| y) =

dz P(z|

y) I E(x| y, z) =

dz P(z|

y)

 zCu(zCu + Cw)−1

y

 

ctr

where P(z| y) = P( y|z) P(z) P y , P( y|z) = exp(− yT(zCu + Cw)−1 y/2)

• (2π)N|zCu + Cw|

Numerical computation of solution is reasonably efficient if

• ne jointly diagonalizes Cu and Cw ...

[Portilla, Strela, Wainwright, Simoncelli, ’03]

IPAM, 9/04 20

Example estimators

Estimators for the scalar and single-neighbor cases

NOISY COEFF. ESTIMATED COEFF. !w

!!" " !" !#" " #" !!" " !" $%&'()*%+,,- $%&'()./0+$1 +'1&2/1+3)*%+,,-

[Portilla etal 03]

Comparison to other methods

Results averaged over 3 images

!" #" $" %" &" !$'& !$ !#'& !# !!'& !! !"'& " "'& ()*+,-*.)/-0123 /456,(74-.)/-0123 )89!:1:;<=>? .=9@A?-9?=@BC8D !" #" $" %" &" !$'& !$ !#'& !# !!'& !! !"'& " "'& ()*+,-*.)/-0123 ,456748+91:;<= :>6965#8*>?6<

[Portilla etal 03]

Original Noisy (22.1 dB) Matlab’s wiener2 (28 dB) BLS-GSM (30.5 dB)

Original Noisy (8.1 dB) UndWvlt Thresh (19.0 dB) BLS-GSM (21.2 dB)

Real sensor noise

400 ISO denoised

GSM summary

• GSM captures local variance
• Underlying Gaussian leads to simple computation
• Excellent denoising results
• What’s missing?
• Global model of z variables [Wainwright etal 99;

Romberg etal ‘99; Hyvarinen/Hoyer ‘02; Karklin/ Lewicki ‘02; Lyu/Simoncelli 08]

• Explicit geometry: phase and orientation

Global models for z

• Non-overlapping neighborhoods, tree-structured z

[Wainwright etal 99; Romberg etal ’99]

• Field of GSMs: z is an exponentiated GMRF, is

a GMRF, subband is the product

[Lyu&Simoncelli 08]

z

Coarse scale Fine scale

• u
• u

Lena Boats

! "#"! $! !# %! "## !& !' !$ !" # "

σ

∆()*+,

! "#"! $! !# %! "## !& !' !$ !" # "

σ

∆()*+,

FoE kSVD GSM BM3D FoGSM

[Lyu&Simoncelli, PAMI 08]

State-of-the-art denoising

2-band steerable pyramid: Image decomposition in terms of multi-scale gradient measurements

Measuring Orientation

[Simoncelli et.al., 1992; Simoncelli & Freeman 1995]

Multi-scale gradient basis

Multi-scale gradient basis

• Multi-scale bases: efficient representation

Multi-scale gradient basis

• Multi-scale bases: efficient representation
• Derivatives: good for analysis
• Local Taylor expansion of image structures
• Explicit geometry (orientation)

Multi-scale gradient basis

• Multi-scale bases: efficient representation
• Derivatives: good for analysis
• Local Taylor expansion of image structures
• Explicit geometry (orientation)
• Combination:
• Explicit incorporation of geometry in basis
• Bridge between PDE / harmonic analysis

approaches

• rientation
• rientation

magnitude

[Hammond&Simoncelli 06; cf. Oppenheim and Lim 81]

Importance of local orientation

Randomized orientation Randomized magnitude

[Hammond&Simoncelli 05]

Reconstruction from orientation

• Reconstruction by projections onto convex sets
• Resilient to quantization

Quantized to 2 bits

[Hammond&Simoncelli 06]

Original

Image patches related by rotation

two-band steerable pyramid coefficients

[Hammond&Simoncelli 06]

raw patches rotated patches

• -- Raw Patches

Rotated Patches

PCA of normalized gradient patches

[Hammond&Simoncelli 06]

Orientation-Adaptive GSM model

patch rotation operator hidden magnitude/orientation variables Model a vectorized patch of wavelet coefficients as:

[Hammond&Simoncelli 06]

Orientation-Adaptive GSM model

patch rotation operator hidden magnitude/orientation variables Model a vectorized patch of wavelet coefficients as: Conditioned on ; is zero mean gaussian with covariance

[Hammond&Simoncelli 06]

[Hammond&Simoncelli 06]

Estimation of C(θ) from noisy data

noisy patch unknown, approximate by measured from noisy data. Assuming independent and noise rotationally invariant (assuming w.l.o.g. E[z] =1 )

Bayesian MMSE Estimator

[Hammond&Simoncelli 06]

Bayesian MMSE Estimator

condition on and integrate

• ver hidden variables

[Hammond&Simoncelli 06]

Bayesian MMSE Estimator

condition on and integrate

• ver hidden variables

[Hammond&Simoncelli 06]

Bayesian MMSE Estimator

condition on and integrate

• ver hidden variables

Wiener estimate

[Hammond&Simoncelli 06]

Bayesian MMSE Estimator

condition on and integrate

• ver hidden variables

has covariance separable prior for hidden variables Wiener estimate

[Hammond&Simoncelli 06]

Bayesian MMSE Estimator

condition on and integrate

• ver hidden variables

has covariance separable prior for hidden variables Wiener estimate

[Hammond&Simoncelli 06]

• agsm

13.1 dB gsm2 12.4 dB σ = 40 noisy 2.81 dB

Locally adaptive covariance

• Karklin & Lewicki 08: Each patch is Gaussian,

with covariance constructed from a weighted outer- product of fixed vectors:

• Guerrero-Colon, Simoncelli & Portilla 08: Each

patch is a mixture of GSMs (MGSMs):

p( x) =

• k

Pk

• p(zk) G(

x; zkCk) dzk p( y) =

• n

exp(−|yn|) Bn =

• k

wnk bk bk

T

log C( y) =

• n

ynBn p( x) = G ( x; C( y))

MGSMs generative model

• x

Patch chosen from with probabilities Parameters:

• Covariances
• Scale densities
• Component probabilities
• Number of components

Ck {√z1 u1, √z2 u2, ...√zK uK} {P1, P2, ..., PK} pk(zk) Pk K Parameters can be fit to data of one or more images by maximizing likelihood (EM-like)

[Guerrero-Colon, Simoncelli, Portilla 08]

MGSM “segmentation”

First six eigenvectors

• f GSM

covariance matrices

image 1 2 4

[Guerrero-Colon, Simoncelli, Portilla 08]

MGSM “segmentation”

Eigenvectors of GSM components represent invariant subspaces: “generalized complex cells”

Potential of local homogeneous models?

• marginal statistics [var,skew,kurtosis]
• local raw correlations
• local variance correlations
• local phase correlations

Consider an implicit model: maxEnt subject to constraints on subband coefficients:

[Portilla & Simoncelli 00;

• cf. Zhu, Wu & Mumford 97]

Visual texture

Visual texture

Homogeneous, with repeated structures

Visual texture

Homogeneous, with repeated structures “You know it when you see it”

All Images

Texture Images Equivalence class (visually indistinguishable)

Iterative synthesis algorithm

Synthesis Analysis

Transform

Measure Statistics

Example Texture Random Seed Synthesized Texture

Transform

Measure Statistics

Adjust Inverse Transform

[Portilla&Simoncelli 00; cf. Heeger&Bergen ‘95]

Examples: Artificial

Photographic, quasi-periodic

Photographic, aperiodic

Photographic, structured

Photographic, color

Non-textures?

Texture mixtures

Texture mixtures

Convex combinations in parameter space

Texture mixtures

Convex combinations in parameter space => Parameter space includes non-textures

Summary

• Fusion of empirical data with structural principles
• Statistical models have led to state-of-the-art image

processing, and are relevant for biological vision

• Local adaptation to {variance, orientation,

phase, ...} gives improvement, but makes learning harder

• Cascaded representations emerge naturally
• There’s still much room for improvement!

Cast

• Local GSM model: Martin Wainwright, Javier Portilla
• GSM Denoising: Javier Portilla, Martin Wainwright,

Vasily Strela

• Variance-adaptive compression: Robert Buccigrossi
• Local orientation and OAGSM: David Hammond
• Field of GSMs: Siwei Lyu
• Mixture of GSMs: Jose-Antonio Guerrero-Colón,

Javier Portilla

• Texture representation/synthesis: Javier Portilla