Prior Learning and Gibbs Reaction-Diffusion.
AbstractThis article addresses two important themes in early visual computation: First, it presents a novel theory for learning the universal statistics of natural imagesa prior model for typical cluttered scenes of the worldfrom a set of natural images, and, second, it proposes a general framework of designing reaction-diffusion equations for image processing. We start by studying the statistics of natural images including the scale invariant properties, then generic prior models were learned to duplicate the observed statistics, based on the minimax entropy theory studied in two previous papers. The resulting Gibbs distributions have potentials of the form $U\left( {{\schmi{\bf I}};\,\Lambda ,\,S} \right)=\sum\nolimits_{\alpha =1}^K {\sum\nolimits_{ \left( {x,y} \right)} {\lambda ^{\left( \alpha \right)}}}\left( {\left( {F^{\left( \alpha \right)}*{\schmi{\bf I}}} \right)\left( {x,y} \right)} \right)$ with being a set of filters and the potential functions. The learned Gibbs distributions confirm and improve the form of existing prior models such as line-process, but, in contrast to all previous models, inverted potentials (i.e., (x) decreasing as a function of |x|) were found to be necessary. We find that the partial differential equations given by gradient descent on U(I; , S) are essentially reaction-diffusion equations, where the usual energy terms produce anisotropic diffusion, while the inverted energy terms produce reaction associated with pattern formation, enhancing preferred image features. We illustrate how these models can be used for texture pattern rendering, denoising, image enhancement, and clutter removal by careful choice of both prior and data models of this type, incorporating the appropriate features.
texture pattern rendering, denoising, image enhancement and clutter removal by careful
choice of both prior and data models of this type, incorporating the appropriate features.
Song Chun Zhu is now with the Computer Science Department, Stanford University,
Stanford, CA 94305, and David Mumford is with the Division of Applied Mathematics,
Brown University, Providence, RI 02912. This work started when the authors were at
Harvard University.
1 Introduction and motivation
In computer vision, many generic prior models have been proposed to capture the
universal low level statistics of natural images. These models presume that surfaces
of objects be smooth, and adjacent pixels in images have similar intensity values
unless separated by edges, and they are applied in vision algorithms ranging from
image restoration, motion analysis, to 3D surface reconstruction.
For example, in image restoration general smoothness models are expressed as
probability distributions [9, 4, 20, 11]:
where I is the image, Z is a normalization factor, and r x I(x;
r y I(x; are differential operators. Three typical forms of
the potential function /() are displayed in figure 1. The functions in figure 1b, 1c
have flat tails to preserve edges and object boundaries, and thus they are said to
have advantages over the quadratic function in figure 1a.

Figure

Three existing forms for /(). a, Quadratic:
2 . b, Line process:
These prior models have been motivated by regularization theory [26, 18], 1 phys-
1 Where the smoothness term is explained as a stabilizer for solving "ill-posed" problems [32].
ical modeling [31, 4], 2 Bayesian theory [9, 20] and robust statistics [19, 13, 3]. Some
connections between these interpretations are also observed in [12, 13] based on
effective energy in statistics mechanics. Prior models of this kind are either generalized
from traditional physical models [37] or chosen for mathematical convenience.
There is, however, little rigorous theoretical or empirical justification for applying
these prior models to generic images, and there is little theory to guide the construction
and selection of prior models. One may ask the following questions.
1. Why are the differential operators good choices in capturing image features?
2. What is the best form for p(I) and /()?
3. A relevant fact is that real world scenes are observed at more or less arbitrary
scales, thus a good prior model should remain the same for image features
at multiple scales. However none of the existing prior models has the scale-invariance
property on the 2D image lattice, i.e., is renormalizable in terms
of renormalization group theory [36].
In previous work on modeling textures, we proposed a new class of Gibbs distributions
of the following form [40, 41],
e \GammaU (I; ;S) ; (2)
In the above equation, is a set of linear filters, and
is a set of potential functions on the features extracted
by S. The central property of this class of models is that they can reproduce the
marginal distributions of F (ff)   I estimated over a set of the training images I - while
having the maximum entropy - and the best set of features fF (1) ; F (2) ; :::; F (K) g is
If /() is quadratic, then variational solutions minimizing the potential are splines, such as
flexible membrane or thin plate models.
selected by minimizing the entropy of p(I) [41]. The conclusion of our earlier papers
is that for an appropriate choice of a small set of filters S, random samples from
these models can duplicate very general classes of textures - as far as normal human
perception is concerned. Recently we found that similar ideas of model inference
using maximum entropy have also been used in natural language modeling[1].
In this paper, we want to study to what extent probability distributions of this
type can be used to model generic natural images, and we try to answer the three
questions raised above.
We start by studying the statistics of a database of 44 real world images, and then
we describe experiments in which Gibbs distributions in the form of equation (2)
were constructed to duplicate the observed statistics. The learned potential functions
can be classified into two categories: diffusion terms
which are similar to figure 1c, and reaction terms which, in contrast to all previous
models, have inverted potentials (i.e. -(x) decreasing as a function of jxj).
We find that the partial differential equations given by gradient descent on
are essentially reaction-diffusion equations, which we call the Gibbs Reaction
And Diffusion Equations (GRADE). In GRADE, the diffusion components
produce denoising effects which is similar to the anisotropic diffusion [25], while
reaction components form patterns and enhance preferred image features.
The learned prior models are applied to the following applications.
First, we run the GRADE starting with white noise images, and demonstrate
how GRADE can easily generate canonical texture patterns such as leopard blobs
and zebra stripe, as the Turing reaction-diffusion equations do[34, 38]. Thus our
theory provides a new method for designing PDEs for pattern synthesis.
Second, we illustrate how the learned models can be used for denoising, image
enhancement and clutter removal by careful choice of both prior and noise models
of this type, incorporating the appropriate features extracted at various scales
and orientations. The computation simulates a stochastic process - the Langevin
equations - for sampling the posterior distribution.
This paper is arranged as follows. Section (2) presents a general theory for prior
learning. Section (3) demonstrates some experiments on the statistics of natural
images and prior learning. Section (4) studies the reaction-diffusion equations.
Section (5) demonstrates experiments on denoising, image enhancement and clutter
removal. Finally section (6) concludes with a discussion.
2 Theory of prior learning
2.1 Goal of prior learning and two extreme cases
We define an image I on an N \Theta N lattice L to be a function such that for any
pixel (x; y), I(x; y) 2 L, and L is either an interval of R or L ae Z. We assume
that there is an underlying probability distribution f(I) on the image space L N 2
for general natural images - arbitrary views of the world. Let NI
2::; Mg be a set of observed images which are independent samples from f(I).
The objective of learning a generic prior model is to look for common features and
their statistics from the observed natural images. Such features and their statistics
are then incorporated into a probability distribution p(I) as an estimation of f(I),
so that p(I), as a prior model, will bias vision algorithms against image features
which are not typical in natural images, such as noise distortion and blurring. For
this objective, it is reasonable to assume that any image features have equal chance
to occur at any location, so f(I) is translation invariant with respect to (x; y). We
will discuss the limits of this assumption in section (6).
To study the properties of images fI obs
we start from exploring
a set of linear filters which are characteristic
of the observed images. The statistics extracted by S are the empirical marginal
distributions (or histograms) of the filter responses.
Given a probability distribution f(I), the marginal distribution of f(I)
with respect to F (ff) is,
Z
Z
F (ff) \LambdaI(x;y)=z
where 8z 2 R and ffi() is a Dirac function with point mass concentrated at 0.
Given a linear filter F (ff) and an image I, the empirical marginal
distribution (or histogram) of the filtered image F (ff)   I(x; y) is,
We compute the histogram averaged over all images in NI obs as the observed statistics

obs (z) =M
H (ff) (z; I obs
If we make a good choice of our database, then we may assume that - (ff)
obs (z) is
an unbiased estimate for f (ff) (z), and as M ! 1, - (ff)
obs (z) converges to f (ff)
Now, to learn a prior model from the observed images fI obs
immediately we have two simple solutions. The first one is,
Y
obs is the observed average histogram of the image intensities, i.e., the filter
used. Taking / 1
obs (z), we rewrite equation (4) as,
The second solution is:
Let kI obs
to the Potts model [37].
These two solutions stand for two typical mechanisms for constructing probability
models in the literature. The first is often used for image coding [35], and the
second is a special case of the learning scheme using radial basis functions (RBF)
[27]. 3 Although the philosophies for learning these two prior models are very differ-
ent, we observe that they share two common properties.
1. The potentials / 1 (), / 2 () are built on the responses of linear filters. In equation
(7), I obs
are used as linear filters of size N \Theta N pixels,
which we denote by F
n .
2. For each filter F (ff) chosen, p(I) in both cases duplicates the observed marginal
distributions. It is trivial to prove that E p [H (ff) (z;
obs (z), thus as M
increases,
This second property is in general not satisfied by existing prior models in equation
(1). p(I) in both cases meets our objective for prior learning, but intuitively
they are not good choices. In equation (5), the ffi() filter does not capture spatial
structures of larger than one pixels, and in equation (7), filters F (obsn) are too
specific to predict features in unobserved images.
In fact, the filters used above lie in the two extremes of the spectrum of all linear
filters. As discussed by Gabor [7], the ffi filter is localized in space but is extended
uniformly in frequency. In contrast, some other filters, like the sine waves, are well
3 In RBF, the basis functions are presumed to be smooth, such as a Gaussian function. Here,
using ffi () is more loyal to the observed data.
localized in frequency but are extended in space. Filter F (obsn) includes a specific
combination of all the components in both space and frequency. A quantitative
analysis of the goodness of these filters is given in table 1 at section (3.2).
2.2 Learning prior models by minimax entropy
To generalize the two extreme examples, it is desirable to compute a probability
distribution which duplicates the observed marginal distributions for an arbitrary
set of filters, linear or nonlinear. This goal is achieved by a minimax entropy theory
studied for modeling textures in our previous papers [40, 41].
Given a set of filters fF (ff) observed statistics f- (ff)
obs
Kg, a maximum entropy distribution is derived which has the following
Gibbs form:
In the above equation, we consider linear filters only, and
is a set of potential functions on the features extracted by S. In practice, image
intensities are discretized into a finite gray levels, and the filter responses are divided
into a finite number of bins, therefore - (ff) () is approximated by a piecewisely
constant functions, i.e., a vector, which we denote by - (ff)
The - (ff) 's are computed in a non-parametric way so that the learned p(I;  ; S)
reproduces the observed statistics:
Therefore as far as the selected features and their statistics are concerned, we cannot
distinguish between p(I;  ; S) and the "true" distribution f(I).
Unfortunately, there is no simple way to express the - (ff) 's in terms of the - (ff)
obs 's
as in the two extreme examples. To compute - (ff) 's, we adopted the Gibbs sampler
[9], which simulates an inhomogeneous Markov chain in image space L jN 2 j . This
Monte Carlo method iteratively samples from the distribution p(I;  ; S), followed
by computing the histogram of the filter responses for this sample and updating the
- (ff) to bring the histograms closer to the observed ones. For a detailed account of
the computation of - (ff) 's, the readers are referred to [40, 41].
In our previous papers, the following two propositions are observed.
Proposition 1 Given a filter set S, and observed statistics f- (ff)
there is an unique solution for Kg.
Proposition 2 f(I) is determined by its marginal distributions, thus
if it reproduces all the marginal distributions of linear filters.
But for computational reasons, it is often desirable to choose a small set of filters
which most efficiently capture the image structures. Given a set of filters S, and
an ME distribution p(I;  ; S), the goodness of p(I;  ; S) is often measured by the
Kullback-Leibler information distance between p(I;  ; S) and the ideal distribution
Z
Z
f(I) log f(I)
Then for a fixed model complexity K, the best feature set S   is selected by the
following criterion,
where S is chosen from a general filter bank B such as Gabor filters at multiple
scales and orientations.
Enumerating all possible sets of features S in the filter bank and comparing
their entropies is computational too expensive. Instead, in [41] we propose a step-wise
greedy procedure for minimizing the KL-distance. We start from
a uniform distribution, and introduce one filter at a time. Each added
filter is chosen to maximally decrease the KL-distance, and we keep doing this until
the decrease is smaller than a certain value.
In the experiments of this paper, we have used a simpler measure of the "infor-
mation gain" achieved by adding one new filter to our feature set S. This is roughly
an L 1 -distance (vs. the L 2 -measure implicit in the Kullback-Leibler distance), the
readers are referred to [42] for a detailed account).
and p(I;  ; S) defined above, the information
criterion (IC) for each filter F (fi) 2 B=S at step K
kH (fi) (z; I obs
kH (fi) (z; I obs
obs (z)k
we call the first term average information gain (AIG) by choosing F (fi) , and the
second term average information fluctuation (AIF ).
Intuitively, AIG measures the average error between the filter responses in the
database and the marginal distribution of the current model p(I;  ; S). In practice,
we need to sample p(I;  ; S), thus synthesize images fI syn
estimate E p(I; ;S) [H (fi) (z; I)] by - (fi)
n ). For a filter F (fi) , the
bigger AIG is, the more information F (fi) captures as it reports the error between
the current model and the observations. AIF is a measure of disagreement between
the observed images. The bigger AIF is, the less their responses to F (ff) have in
common.
3 Experiments on natural images
This section presents experiments on learning prior models, and we start from exploring
the statistical properties of natural images.

Figure

6 out of the 44 collected natural images.
3.1 Statistic of natural images
It is well known that natural images have statistical properties distinguishing them
from random noise images [28, 6, 24]. In our experiments, we collected a set of 44
natural images, six of which are shown in figure 2. These images are from various
sources, some digitized from books and postcards, and some from a Corel image
database. Our database includes both indoor and outdoor pictures, country and
urban scenes, and all images are normalized to have intensity between 0 and 31.
As stated in proposition (2), marginal distributions of linear filters alone are
capable of characterizing f(I). In the rest of this paper, we shall only study the
following bank B of linear filters.
1. An intensity filter ffi().
2. Isotropic center-surround filters, i.e., the Laplacian of Gaussian filters.
const
2oe stands for the scale of the filter. We denote these filters by LG(s).
A special filter is LG(
), which has a 3 \Theta 3 window [0;
we denote it by \Delta.
3. Gabor filters with both sine and cosine components, which are models for the
frequency and orientation sensitive simple cells.
const
It is a sine wave at frequency 2-
s modulated by an elongated Gaussian function, and
rotated at angle '. We denote the real and image parts of G(x;
and Gsin(s; '). Two special Gsin(s; ') filters are the gradients r x ; r y .
4. We will approximate large scale filters by filters of small window sizes on the
high level of the image pyramid, where the image in one level is a "blown-down"
version (i.e., averaged in 2 \Theta 2 blocks) of the image below.
We observed three important aspects of the statistics of natural images.
First, for some features, the statistics of natural images vary widely from image
to image. We look at the ffi() filter as in section (2.1). The average intensity
histogram of the 44 images - (o)
obs is plotted in figure 3a, and figure 3b is the intensity
histogram of an individual image (the temple image in figure 2). It appears that
obs (z) is close to a uniform distribution (figure 3c), while the difference between
figure 3a and figure 3b is very big. Thus IC for filter ffi() should be small (see table
1).
Second, for many other filters, the histograms of their responses are amazingly
consistent across all 44 natural images, and they are very different from the histograms
of noise images. For example, we look at filter r x . Figure 4a is the average
histogram of 44 filtered natural images, figure 4b is the histogram of an individual
filtered image (the same image as in figure 3b), and figure 4c is the histogram of a
filtered uniform noise image.
The average histogram in figure 4a is very different from a Gaussian distribution.

Figure

3 The intensity histograms in domain [0; 31], a, averaged over 44 natural images, b, an individual
natural image, c, a uniform noise image.

Figure

4 The histograms of r x I plotted in domain [-30, 30], a. averaged over 44 natural images, b,
an individual natural image, c, a uniform noise image.

Figure

5 a. The histogram of r x I plotted against Gaussian curve (dashed) of same mean and variance
in domain [\Gamma15; 15]. b, The logarithm of the two curves in a.
To see this, figure 5a plots it against a Gaussian curve (dashed) of the same mean
and same variance. The histogram of natural images has higher kurtosis and heavier
tails. Similar results are reported in [6]. To see the difference of the tails, figure 5b
plots the logarithm of the two curves.
Third, the statistics of natural images are essentially scale invariant with respect
to some features. As an example, we look at filters r x and r y . For each image
I obs
build a pyramid with I [s]
n being the image at the s-th layer. We set
I [0]
n , and let
I [s+1]
I [s]
The size of I [s]
n is N=2 s \Theta N=2 s .
-22
a b c

Figure

2. b. log - x;s
c. histograms of a filtered uniform noise image at scales:
curve), and (dashed curve).
For the filter r x , let - x;s (z) be the average histogram of r x I [s]
Figure 6a plots - x;s (z), for and they are almost identi-
cal. To see the tails more clearly, we display log - x;s (z); in figure 6c.
The difference between them are still small. Similar results are observed for - y;s (z)
the average histograms of r y I obs
n . In contrast, figure 6b plots the histograms
of r x I [s] with I [s] being a uniform noise image at scales 2.
Combining the second and the third aspects above, we conclude that the histograms
of r x I [s]
are very consistent across all observed natural images and
across scales 2. The scale invariant property of natural images is largely
caused by the following facts: 1). natural images contains objects of all sizes, 2).
natural scenes are viewed and made into images at arbitrary distances.
3.2 Empirical prior models
In this section, we learn the prior models according to the theory proposed in
section (2), and analyze the efficiency of the filters quantitatively.
Experiment I.
a b c

Figure

7 The three learned potential functions for filters a. \Delta, b. r x , and c. r y . Dashed curves are the
fitting functions: a. / 1
and c. / 3
We start from We compute the AIF ,
AIG and IC for all filters in our filter bank. We list the results for a small number
of filters in table 1. The filter \Delta has the biggest IC (= 0:642), thus is chosen as
F (1) . An ME distribution p 1 (I;  ; S) is learned, and the information criterion for
each filter is shown in the column headed p 1 (I) in table 1. We notice that the IC for
the filter \Delta drops to near 0, and IC also drops for other filters because these filters
are in general not independent of \Delta. Some small filters like LG(1) have smaller ICs
than others, due to higher correlations between them and \Delta.

Figure

8 A typical sample of p 3 (I) (256 \Theta 256 pixels).
The big filters with larger IC are investigated in Experiment II. In this experi-
ment, we choose both r x and r y to be F (2) ; F (3) as in other prior models. Therefore
a prior model p 3 (I) is learned with potential:
are plotted in figure 7. Since - (1)
we only plot - (1) (z) for z 2 [\Gamma9:5; 9:5] and
- (2) (z); - (3) (z) for z 2 [\Gamma22; 22]. These three curves are fitted with the functions
synthesized image sampled
from p 3 (I) is displayed in figure 8.
So far, we have used three filters to characterize the statistics of natural images,
and the synthesized image in figure 8 is still far from natural ones. Especially,
even though the learned potential functions - (ff) (z); tails to
4 In fact, - (1)
obs
obs
, with N \Theta N being the size of synthesized image.
Filter Filter Size AIF AIG IC AIG IC AIG IC AIG IC
I

Table

1 The information criterion for filter selection.
preserve intensity breaks, they only generate small speckles instead of big regions
and long edges as one may expect. Based on this synthesized image, we compute
the AIG and IC for all filters, and the results are list in table 1 in column p 3 (I).
Experiment II.
It is clear that we need large-scale filters to do better. Rather than using the
large scale Gabor filters, we chose to use r x and r y on 4 different scales and
impose explicitly the scale invariant property that we find in natural images. Given
an image I defined on an N \Theta N lattice L, we build a pyramid in the same way as
before. Let I 3 be four layers of the pyramid. Let H x;s (z; x; y) denote
the histogram of r x I [s] (x; y) and H y;s (z; x; y) the histogram of r y I [s] (x; y).
We ask for a probability model p(I) which satisfies,
3:
3:
where L s is the image lattice at level s, and -
-(z) is the average of the observed
histograms of r x I [s] and r y I [s] on all 44 natural images at all scales. This results
in a maximum entropy distribution p s (I) with energy of the following form,
U s (I) =X

Figure

3. At the beginning of the learning process,
all - x;s (); are of the form displayed in figure 7 with low values around
zero to encourage smoothness. As the learning proceeds, gradually - x;3 () turns "up-
side down" with smaller values at the two tails. Then - x;2 () and - x;1 () turn upside
down one by one. Similar results are observed for - y;s (); 3. Figure 11 is
a typical sample image from p s (I). To demonstrate it has scale invariant properties,
in figure 10 we show the histograms H x;s and log H x;s of this synthesized image for
3.
The learning process iterates for more than 10; 000 sweeps. To verify the learned
-()'s, we restarted a homogeneous Markov chain from a noise image using the
learned model, and found that the Markov chain goes to a perceptually similar
image after 6000 sweeps.
Remark 1. In figure 9, we notice that - x;s () are inverted, i.e. decreasing
functions of j z j for distinguishing this model from other prior models
in computer vision. First of all, as the image intensity has finite range [0; 31], r x I [s]
is defined in [\Gamma31; 31]. Therefore we may define - x;s
still well-defined. Second, such inverted potentials have significant meaning in visual
computation. In image restoration, when a high intensity difference r x I [s] (x; y)
is present, it is very likely to be noise if However this is not true for
3. Additive noise can hardly pass to the high layers of the pyramid because
at each layer the 2 \Theta 2 averaging operator reduces the variance of the noise by 4
times. When r x I [s] (x; y) is large for it is more likely to be a true
a b
-5
c d

Figure

9 Learned potential functions - x;s (); 3. The dashed curves are fitting functions:
-22
a b

Figure

a. The histograms of the synthesized image at 4 scales-almost indistinguishable. b. The
logarithm of histograms in a.

Figure

11 A typical sample of p s (I) (384 \Theta 384 pixels).
edge and object boundary. So in p s (I), - x;0 () suppresses noise at the first layer,
while - x;s (); encourages sharp edges to form, and thus enhances blurred
boundaries. We notice that regions of various scales emerge in figure 11, and the
intensity contrasts are also higher at the boundary. These appearances are missing
in figure 8.
Remark 2. Based on the image in figure 11, we computed IC and AIG for
all filters and list them under column p s (I) in table 1. We also compare the two
extreme cases discussed in section (2.1). For the ffi() filter, AIF is very big, and AIG
is only slightly bigger than AIF . Since all the prior models that we learned have
no preference about the image intensity domain, the image intensity has uniform
distribution, but we limit it inside [0; 31], thus the first row of table 1 has the same
value for IC and AIG. For filter I (obsi) ,
M i.e. the biggest among all
filters, and AIG ! 1. In both cases, ICs are the two smallest.
4 Gibbs reaction-diffusion equations
4.1 From Gibbs distribution to reaction-diffusion equations
The empirical results in the previous section suggest that the forms of the potentials
learned from images of real world scenes can be divided into two classes:
upright curves -(z) for which -() is an even function increasing as jzj increases and
inverted curves for which the opposite happens. Similar phenomenon was observed
in our learned texture models [40].
In figure 9, - x;s (z) are fit to the family of functions (see the dashed curves),
are respectively the translation and scaling constants, and kak weights the
contribution of the filter.
In general the Gibbs distribution learned from images in a given application has
potential function of the following form,
OE (ff)
ff=n d +1
Note that the filter set is now divided into two parts
Kg. In most cases S d consists
of filters such as r x ; r y ; \Delta which capture the general smoothness of images, and S r
contains filters which characterize the prominent features of a class of images, e.g.
Gabor filters at various orientations and scales which respond to the larger edges
and bars.
Instead of defining a whole distribution with U , one can use U to set up a
variational problem. In particular, one can attempt to minimize U by gradient
descent. This leads to a non-linear parabolic partial differential equation:
I
F (ff)
ff=n d +1
F (ff)
with F (ff)
\Gammay). Thus starting from an input image I(x;
I in , the first term diffuses the image by reducing the gradients while the second term
forms patterns as the reaction term. We call equation (14) the Gibbs Diffusion And
Reaction Equation (GRADE).
Since the computation of equation (14) involves convolving twice for each of the
selected filters, a conventional way for efficient computation is to build an image
pyramid so that filters with large scales and low frequencies can be scaled down into
small ones in the higher level of the image pyramid. This is appropriate especially
when the filters are selected from a bank of multiple scales, such as the Gabor filters
and the wavelet transforms. We adopt this representation in our experiments.
For an image I, let I [s] be an image at level of a pyramid, and
I I, the potential function becomes,
s
ff
s
s
ff
s
s   I [s] (x; y))
We can derive the diffusion equations similarly for this pyramidal representation.
4.2 Anisotropic diffusion and Gibbs reaction-diffusion
This section compares GRADEs with previous diffusion equations in vision.
In [25, 23] anisotropic diffusion equations for generating image scale spaces are
introduced in the following form,
I
where div is the divergence operator, i.e., div( ~
and Malik defined the heat conductivity c(x; as functions of local gradients,
for example:
I
I x
I y ); (16)
Equation (16) minimizes the energy function in a continuous form,
Z Z
are plotted in figure 12. Similar
forms of the energy functions are widely used as prior distributions [9, 4, 20, 11],
and they can also be equivalently interpreted in the sense of robust statistics [13, 3]
x
-0.4
a

Figure

In the following, we address three important properties of the Gibbs reaction-diffusion
equations.
First, we note that equation (14) is an extension to equation (15) on a discrete
lattice by defining a vector field,
~
and a divergence operator,
Thus equation (14) can be written as,
I
Compared to equation (15) which transfers the "heat" among adjacent pixels, equation
transfers the "heat" in many directions in a graph, and the conductivities
are defined as functions of the local patterns not just the local gradients.
Second, in figure 13, OE(-) has round tip for fl - 1, and a cusp occurs at
(0) can be any value in (\Gamma1; 1) as
shown by the dotted curves in figure 13d. An interesting fact is that the potential
function learned from real world images does have a cusp as shown in figure 9a
where the best fit is 0:7. This cusp forms because large part of objects in real
world images have flat intensity appearances, and OE(-) with produce
piecewise constant regions with much stronger forces than fl - 1.
By continuity, OE 0
(-) can be assigned any value in the range [\Gamma!; !] for - 2 [\Gammaffl; ffl]
In numerical simulations, for - 2 [\Gamma!; !] we take
\Gammaoe if oe 2 [\Gamma!; !]
a
-0.4
-0.4
x

Figure

13 The potential function
(-). a,c,
where oe is the summation of the other terms in the differential equation whose
values are well defined. Intuitively when
(0) forms an attractive basin in its neighborhood N (ff) (x; y) specified by the
filter window of F (ff) . For a pixel (u; v) 2 N (ff) (x; y), the depth of the attractive
basin is k!F (ff)
a pixel is involved in multiple zero filter responses,
it will accumulate the depth of the attractive basin generated by each filter. Thus
unless the absolute value of the driving force from other well-defined terms is larger
than the total depth of the attractive basin at (u; v), I(u; v) will stay unchanged. In
the image restoration experiments in later sections, performance
in forming piecewise constant regions.
Third, the learned potential functions confirmed and improved the existing prior
models and diffusion equations, but more interestingly reaction terms are first dis-
covered, and they can produce patterns and enhance preferred features. We will
demonstrate this property in the experiments below.
4.3 Gibbs reaction-diffusion for pattern formation
In the literature, there are many nonlinear PDEs for pattern formation, of which the
following two examples are interesting. (I) The Turing reaction-diffusion equation
which models the chemical mechanism of animal coats [33, 21]. Two canonical patterns
that the Turing equations can synthesize are leopard blobs and zebra stripes
[34, 38]. These equations are also applied to image processing such as image halftoning
[29] and a theoretical analysis can be found in [15]. (II) The Swindale equation
which simulates the development of the ocular dominance stripes in the visual cortex
of cats and monkey [30]. The simulated patterns are very similar to the zebra
stripes.
In this section, we show that these patterns can be easily generated with only
2 or 3 filters using the GRADE. We run equation (14) starting with I(x; as a
uniform noise image, and GRADE converges to a local minimum. Some synthesized
texture patterns are displayed in figure 14.
For all the six patterns in figure 14, we choose F (1)
the Laplacian of
Gaussian filter at level 0 of the image pyramid as the only diffusion filter, and we
fix
(-). For the three patterns in figure 14
a,b,c we choose isotropic center-surround filter LG(
of widow size 7 \Theta 7 pixels
as the reaction filter F (2)
1 at level 1 of the image pyramid, and we set (a = \Gamma6:0;
(-). The differences between these three patterns are caused
by - forms the patterns with symmetric appearances for both
black and white part as shown in figure 14a. As - negative, black blobs
begin to form as shown in figure 14b where - positive -
blobs in the black background as shown in figure 14c where - 6. The general
smoothness appearance of the images is attributed to the diffusion filter. Figure 14d
is generated with two reaction filters - LG(
2) at level 1 and level 2 respectively,
a b c

Figure

14 Leopard blobs and zebra stripes synthesized by GRADEs.
therefore the GRADE creates blobs of mixed sizes. Similarly we selected one cosine
Gabor filter Gcos(4; pixels oriented at 1 as the reaction
filter F (2)

Figure

14f is generated with two reaction filters Gcos(4;
It seems that the leopard blobs and zebra stripes are among the most canonical
patterns which can be generated with easy choices of filters and parameters. As
shown in [40], the Gibbs distribution are capable of modeling a large variety of
texture patterns, but filters and different forms for /(-) have to be learned for a
given texture pattern.
5 Image enhancement and clutter removal
So far we have studied the use of a single energy function U(I) either as the log
likelihood of a probability distribution at I or as a function of I to be minimized by
gradient descent. In image processing, we often need to model both the underlying
images and some distortions, and to maximize a posterior distribution. Suppose the
distortions are additive, i.e., an input image is,
I in = I +C:
In many applications, the distortion images C are often not i.i.d. Gaussian noise,
but clutter with structures such as trees in front of a building or a military target.
Such clutter will be very hard to handle by edge detection and image segmentation
algorithms.
We propose to model clutter by an extra Gibbs distribution, which can be learned
from some training images by the minimax entropy theory as we did for the underlying
image I. Thus an extra pyramidal representation for I in \Gamma I is needed in a
Gibbs distribution form as shown in figure 15. The resulting posterior distributions
are still of the Gibbs form with potential function,
U
where UC () is the potential of the clutter distribution.
Thus the MAP estimate of I is the minimum of U   . In the experiments which
we use the Langevin equation for minimization, a variant of simulated annealing

where w(x; is the standard Brownian motion process, i.e.,
w(x;
T (t) is the "temperature" which controls the magnitude of the random fluctuation.
Under mild conditions on U   , equation (19) approaches a global minimum of U   at
target features
image pyramid
for clutters
image pyramid
for targets
clutter features
observed image
target image clutter image

Figure

15 The computational scheme for removing noise and clutter.
a low temperature. The analyses of convergence of the equations can be found in
[14, 10, 8]. The computational load for the annealing process is notorious, but for
applications like denoising, a fast decrease of temperature may not affect the final
result very much.
Experiment I
In the first experiment, we take UC to be quadratic, i.e. C to be an i.i.d.
Gaussian noise image. We first compare the performance of the three prior models
potential functions are respectively:
U l
the 4-scale energy in equation (12) (22)
l () and / t () are the line-process and T-function displayed in figure 1b and 1c
respectively.

Figure

demonstrates the results: the original image is the lobster boat displayed
in figure 2. It is normalized to have intensity in [0; 31] and Gaussian noise
from N(0; 25) are added. The distorted image is displayed in figure 16a, where
we keep the image boundary noise-free for the convenience of boundary condition.
The restored images using p l (I), p t (I) and p s (I) are shown in figure 16b, 16c, 16d
respectively. p s (I), which is the only model with a reaction term, appears to have
the best effect in recovering the boat, especially the top of the boat, but it also
enhances the water.
Experiment II
In many applications, i.i.d. Gaussian models for distortions are not sufficient.
For example, in figure 17a, the tree branches in the foreground will make image
segmentation and object recognition extremely difficult because they cause strong
edges across the image. Modeling such clutter is a challenging problem in many
applications. In this paper, we only consider clutter as two dimensional pattern
despite its geometry and 3D structure.
We collected a set of images of buildings and a set of images of trees all against
clean background - the sky. For the tree images, we translate the image intensities to
sky. In this case, since the trees are always darker than the build,
thus the negative intensity will approximately take care of the occlusion effects.
We learn the Gibbs distributions for each set respectively in the pyramid, then
such models are respectively adopted as the prior distribution and the likelihood as
in equation (18). We recovered the underlying images by maximizing a posteriori
distribution using the stochastic process.
For example, figure 17b is computed using 6 filters with 2 filters for I: fr
and 4 filters for I C i.e., the potential for I C is,
In the above equation, OE   (-) and /   (-) are fit to the potential functions learned
from the set of tree images,
a
a
c d

Figure

a. The noise distorted image, b. c. d. are respectively restored images by prior models p l
(I).

Figure

17 a. the observed image, b, the restored image using 6 filters.
So the energy term OE   (I(x; y)) forces zero intensity for the clutter image while
allowing for large negative intensities for the dark tree branches.
Figure

18b is computed using 8 filters with 4 filters for I and 4 filters for I C . 13
filters are used for figure 19 where the clutter is much heavier.
As a comparison, we run the anisotropic diffusion process [25] on figure 19a, and
images at iterations are displayed in figure 20. As we can see that
as becomes a flat image. A robust anisotropic diffusion equation is
recently reported in [2].
6 Conclusion
In this paper, we studied the statistics of natural images, based on which a novel
theory is proposed for learning the generic prior model - the universal statistics of
real world scenes. We argue that the same strategy developed in this paper can
be used in other applications. For example, learning probability models for MRI

Figure

a. an observed image, b. the restored image using 8 filters.
a b

Figure

19 a. the observed image, b. the restored image using 13 filters.
a

Figure

20 Images by anisotropic diffusion at iteration
images and 3D depth maps.
The learned prior models demonstrate some important properties such as the
"inverted" potentials terms for patterns formation and image enhancement. The
expressive power of the learned Gibbs distributions allow us to model structured
noise-clutter in natural scenes. Furthermore our prior learning method provides a
novel framework for designing reaction-diffusion equations based on the observed
images in a given application, without modeling the physical or chemical processes
as people did before [33].
Although the synthesized images bear important features of natural images, they
are still far from realistic ones. In other words, these generic prior models can do very
little beyond image restoration. This is mainly due to the fact that all generic prior
models are assumed to be translation invariant, and this homogeneity assumption
is unrealistic. We call the generic prior models studied in this paper the first level
prior. A more sophisticated prior model should incorporate concepts like object
geometry, and we call such prior models second level priors. Diffusion equations
derived from this second level priors are studied in image segmentation [39], and
in scale space of shapes [16]. A discussion of some typical diffusion equations is
given in [22]. It is our hope that this article will stimulate further investigations
on building more realistic prior models as well as sophisticated PDEs for visual
computation.



--R

"A maximum entropy approach to natural language processing"
"Robust anisotropic diffusion"
"On the unification of line processes, outlier re- jection, and robust statistics with applications in early vision"
Visual Reconstruction.

"Relations between the statistics of natural images and the response properties of cortical cells"
"Theory of communication."
"On sampling methods and annealing algorithms"
"Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images"
"Diffusion for global optimization"
"Constrained restoration and the recover of discontinu- ities"
"Parallel and deterministic algorithms for MRFs: surface reconstruction"
"A common framework for image segmentation"
"A renormalization group approach to image processing problems"
The theory and applications of reaction-diffusion equations
"Shapes, shocks, and deformations I: the components of two-dimensional shape and the reaction-diffusion space"
"On information and sufficiency"
"Probabilistic solution of ill-posed problems in computational vision"
"Robust regression methods for computer vision: a review"
"Optimal approximations by piecewise smooth functions and associated variational problems."
"A pre-pattern formation mechanism for mammalian coat markings"
"A general framework for geometry-driven evolution equations"
"Nonlinear image filtering with edge and corner enhance- ment"
"Natural image statistics and efficient coding"
"Scale-space and edge detection using anisotropic diffusion"
"Computational vision and regularization theory"
"Networks for approximation and learning"
"Statistics of natural images: scaling in the woods"
"M-lattice: from morphogenesis to image processing"
"A model for the formation of ocular dominance stripes"
"Multilevel computational processes for visual surface reconstruc- tion"
Solutions of Ill-posed Problems
"The chemical basis of morphogenesis"
"Generating textures on arbitrary surfaces using reaction-diffusion"
"Efficiency of model human image code"
"The renormalization group: critical phenonmena and the Knodo prob- lem,"
Image Analysis
"Reaction-diffusion textures"
"Region Competition: unifying snakes, region grow- ing, Bayes/MDL for multi-band image segmentation"
"Filters, Random Fields, and Minimax Entropy Towards a unified theory for texture modeling"
"Minimax entropy principle and its application to texture modeling"
"Learning generic prior models for visual computa- tion"
--TR

--CTR
Katy Streso , Francesco Lagona, Hidden Markov random field and frame modelling for TCA image analysis, Proceedings of the 24th IASTED international conference on Signal processing, pattern recognition, and applications, p.310-315, February 15-17, 2006, Innsbruck, Austria
Ulf Grenander , Anuj Srivastava, Probability Models for Clutter in Natural Images, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.23 n.4, p.424-429, April 2001
Yufang Bao , Hamid Krim, Smart Nonlinear Diffusion: A Probabilistic Approach, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.26 n.1, p.63-72, January 2004
Dmitry Datsenko , Michael Elad, Example-based single document image super-resolution: a global MAP approach with outlier rejection, Multidimensional Systems and Signal Processing, v.18 n.2-3, p.103-121, September 2007
Giuseppe Boccignone , Mario Ferraro , Terry Caelli, Generalized Spatio-Chromatic Diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.10, p.1298-1309, October 2002
Thang V. Pham , Arnold W. M. Smeulders, Object recognition with uncertain geometry and uncertain part detection, Computer Vision and Image Understanding, v.99 n.2, p.241-258, August 2005
Song-Chun Zhu, Stochastic Jump-Diffusion Process for Computing Medial Axes in Markov Random Fields, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.21 n.11, p.1158-1169, November 1999
Alan L. Yuille , James M. Coughlan, Fundamental Limits of Bayesian Inference: Order Parameters and Phase Transitions for Road Tracking, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.2, p.160-173, February 2000
Ying Nian Wu , Song Chun Zhu , Xiuwen Liu, Equivalence of Julesz Ensembles and FRAME Models, International Journal of Computer Vision, v.38 n.3, p.247-265, July-August 2000
Marc Sigelle, A Cumulant Expansion Technique for Simultaneous Markov Random Field Image Restoration and Hyperparameter Estimation, International Journal of Computer Vision, v.37 n.3, p.275-293, June 2000
Song-Chun Zhu, Embedding Gestalt Laws in Markov Random Fields, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.21 n.11, p.1170-1187, November 1999
Ann B. Lee , David Mumford , Jinggang Huang, Occlusion Models for Natural Images: A Statistical Study of a Scale-Invariant Dead Leaves Model, International Journal of Computer Vision, v.41 n.1-2, p.35-59, January-February 2001
Daniel Cremers , Florian Tischhuser , Joachim Weickert , Christoph Schnrr, Diffusion Snakes: Introducing Statistical Shape Knowledge into the Mumford-Shah Functional, International Journal of Computer Vision, v.50 n.3, p.295-313, December 2002
Song Chun Zhu , Xiu Wen Liu , Ying Nian Wu, Exploring Texture Ensembles by Efficient Markov Chain Monte Carlo-Toward a 'Trichromacy' Theory of Texture, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.6, p.554-569, June 2000
Scott Konishi , Alan L. Yuille , James M. Coughlan , Song Chun Zhu, Statistical Edge Detection: Learning and Evaluating Edge Cues, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.25 n.1, p.57-74, January
J. Sullivan , A. Blake , M. Isard , J. MacCormick, Bayesian Object Localisation in Images, International Journal of Computer Vision, v.44 n.2, p.111-135, September 2001
Norberto M. Grzywacz , Rosario M. Balboa, A Bayesian framework for sensory adaptation, Neural Computation, v.14 n.3, p.543-559, March 2002
Hedvig Sidenbladh , Michael J. Black, Learning the Statistics of People in Images and Video, International Journal of Computer Vision, v.54 n.1-3, p.181-207, August-September
Matthias Heiler , Christoph Schnrr, Natural Image Statistics for Natural Image Segmentation, International Journal of Computer Vision, v.63 n.1, p.5-19, June      2005
Kwang In Kim , Matthias O. Franz , Bernhard Scholkopf, Iterative Kernel Principal Component Analysis for Image Modeling, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.9, p.1351-1366, September 2005
T. Freeman , Egon C. Pasztor , Owen T. Carmichael, Learning Low-Level Vision, International Journal of Computer Vision, v.40 n.1, p.25-47, Oct. 2000
Stefan Roth , Michael J. Black, On the Spatial Statistics of Optical Flow, International Journal of Computer Vision, v.74 n.1, p.33-50, August    2007
Charles Kervrann , Mark Hoebeke , Alain Trubuil, Isophotes Selection and Reaction-Diffusion Model for Object Boundaries Estimation, International Journal of Computer Vision, v.50 n.1, p.63-94, October 2002
Jens Keuchel , Christoph Schnrr , Christian Schellewald , Daniel Cremers, Binary Partitioning, Perceptual Grouping, and Restoration with Semidefinite Programming, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.25 n.11, p.1364-1379, November
Rosario M. Balboa , Norberto M. Grzywacz, The Minimal Local-Asperity Hypothesis of Early Retinal Lateral Inhibition, Neural Computation, v.12 n.7, p.1485-1517, July 2000
A. Srivastava , A. B. Lee , E. P. Simoncelli , S.-C. Zhu, On Advances in Statistical Modeling of Natural Images, Journal of Mathematical Imaging and Vision, v.18 n.1, p.17-33, January
clutter modeling;reaction-diffusion;anisotropic diffusion;gibbs distribution;image restoration;texture synthesis;visual learning