On Discontinuity-Adaptive Smoothness Priors in Computer Vision.
AbstractA variety of analytic and probabilistic models in connection to Markov random fields (MRFs) have been proposed in the last decade for solving low level vision problems involving discontinuities. This paper presents a systematic study of these models and defines a general discontinuity adaptive (DA) MRF model. By analyzing the Euler equation associated with the energy minimization, it shows that the fundamental difference between different models lies in the behavior of interaction between neighboring points, which is determined by the a priori smoothness constraint encoded into the energy function. An important necessary condition is derived for the interaction to be adaptive to discontinuities to avoid oversmoothing. This forms the basis on which a class of adaptive interaction functions (AIFs) is defined. The DA model is defined in terms of the Euler equation constrained by this class of AIFs. Its solution is C1 continuous and allows arbitrarily large but bounded slopes in dealing with discontinuities. Because of the continuous nature, it is stable to changes in parameters and data, a good property for regularizing ill-posed problems. Experimental results are shown.
Introduction
MOOTHNESS is a generic assumption underlying a
wide range of physical phenomena. It characterizes the
coherence and homogeneity of matter within a scope of
space (or an interval of time). It is one of the most common
assumptions in computer vision models, in particular,
those formulated in terms of Markov random fields (MRFs)
[1], [2], [3] and regularization [4]. Its applications are seen
widely in image restoration, surface reconstruction, optical
flow and motion, shape from X, texture, edge detection,
region segmentation, visual integration and so on.
The assumption of the uniform smoothness implies the
smoothness everywhere. However, improper imposition of
it can lead to undesirable, oversmoothed, solutions. This
occurs when the uniform smoothness is violated, for ex-
ample, at discontinuities where abrupt changes occur. It
is necessary to take care of discontinuities when using
smoothness priors. Therefore, how to apply the smoothness
constraint wile preserving discontinuities has been one
of the most active research areas in level vision (see, e.g.
[5], [6], [1], [7], [3], [8], [9], [10], [11], [12], [13], [14], [15],
[16], [17], [18]).
This paper presents a systematic study on smoothness
priors involving discontinuities. The results are based on
an analysis of the Euler equation associated with the energy
minimization in MRF and regularization models. Through
S. Z. Li is currently with the School of Electrical and Electronic
Engineering, Nanyang Technological University, Singapore 2263. E-
mail: szli@ntuix.ntu.ac.sg .
the analysis, it is identified that the fundamental difference
among different models for dealing with discontinuities lies
in their ways of controlling the interaction between neighboring
points. Thereby, an important necessary condition
is derived for any regularizers or MRF prior potential functions
to be able to deal with discontinuities.
Based on these findings, a so-called discontinuity adaptive
(DA) smoothness model is defined in terms of the Euler
equation constrained by a class of adaptive interaction
functions (AIFs). The DA solution is C 1 continuous, allowing
arbitrarily large but bounded slopes. Because of
the continuous nature, it is stable to changes in parameters
and data. This is a good property for regularizing
ill-posed problems. The results provide principles for the
selection of a priori clique potential functions in stochastic
MRF models and regularizers in deterministic regularization
models. It is also shown that the DA model includes
as special instances most of the existing models, such as the
line process (LP) model [1], [3], weak string and membrane
[10], approximations of the LP model [19], [20], minimal description
length [13], biased anisotropic diffusion [16], and
mean field theory approximation [17].
The study of discontinuities is most sensibly carried out
in terms of analytical properties, such as derivatives. For
this reason, analytical regularization, a special class of
MRF models, is used as the platform for it. If we consider
that regularization contains three parts [21]: the data, the
class of solution functions and the regularizer, the present
work addresses mainly the regularizer part. In Section II,
regularization models are reviewed in connection to dis-
continuities. In Section III, a necessary condition for the
discontinuity adaptivity is made explicit; based on this,
the DA model is defined and compared with other models.
In Section IV, an algorithm for finding the DA solution is
presented, some related issues are discussed. Experimental
results are shown in Section V. Finally, conclusions are
drawn.
II. Smoothness, Regularization and
Discontinuities
In MRF vision modeling, the smoothness assumption
can be encoded into an energy via one of the two routes:
analytic and probabilistic. In the analytic route, the encoding
is done in the regularization framework [4], [22].
From the regularization viewpoint, a problem is said to be
"ill-posed" if it fails to satisfy one or more of the following
criteria: the solution exists, is unique and depends continuously
on the data. Additional, a priori, assumptions have
to be imposed on the solution to convert an ill-posed problem
into a well-posed one. An important assumption of
IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-17, NO.6, PP.576-586, JUNE 1995
such assumptions is the smoothness [23]. It is incorporated
into the energy function whereby the cost of the solution
is defined.
From the probabilistic viewpoint, a regularized solution
corresponds to the maximum a posteriori (MAP) estimate
of an MRF [1], [24]. Here, the prior constraints are encoded
into the a priori MRF probability distribution. The MAP
solution is obtained by maximizing the posterior probability
or equivalently minimizing the corresponding energy.
The MRF model is more general than the regularization
model in (1) that it can encode prior constraints other than
the smoothness and (2) that it allows arbitrary neighborhood
systems other than the nearest ones. However, the
analytic regularization model provides a convenient platform
for the study of smoothness priors because of close
relationships between the smoothness and the analytical
continuity.
A. Regularization and Discontinuities
Consider the problem of restoring a signal f from the
data denotes the noise. The regularization
formulation defines the solution f   to be the global
minimum of an energy function E(f
The energy is the sum of two terms
The closeness term, U(d j f ), measures the cost caused by
the discrepancy between the solution f and the data d
a
where -(x) is a weighting function and a and b are the
bounds of the integral. The smoothness term, U(f ), measures
the cost caused by the irregularities of the solution f ,
the irregularities being measured by the derivative magnitudes
jf (n) (x)j. With identical independent additive Gaussian
noise, U(f j d), U(d j f) and U(f) correspond to the
energies in the posterior, the likelihood the prior Gibbs
distributions of an MRF, respectively [24].
The smoothness term U(f ), also called a regularizer, is
the object of study in this work. It penalizes the irregularities
according to the a priori smoothness constraint
encoded in it. It is generally defined as
-n
a
g(f (n) (x))dx (3)
where Un (f) is the n th order regularizer, N is the highest
order to be considered and -n - 0 is a weighting factor.
A potential function g(f (n) (x)) is the penalty against the
irregularity in f (n\Gamma1) (x) and corresponds to prior clique
potentials in MRF models. Regularizers differ in the definition
of Un (f ), more specifically in the selection of g.
A.1 Standard Regularization
In the standard regularization [23], [4], the potential
function takes the pure quadratic form
With g q , the more irregular f (n\Gamma1) (x) is at x, the larger
jf (n) j, and consequently the larger potential g(f (n) ) contributed
to Un (f ). The standard quadratic regularizer can
have a more general form
Un (f; wn
a
wn (x)[f (n) (x)] 2 dx (5)
where wn (x) are the pre-specified non-negative continuous
functions [23]. It may also be generalized to multi-dimensional
cases and to include cross derivative terms.
The quadratic regularizer imposes the smoothness constraint
everywhere. It determines the constant interaction
between neighboring points and leads to smoothing strength
proportional to jf (n) j, as will be shown in the next section.
The homogeneous or isotropic application of the smoothness
constraint inevitably leads to oversmoothing at discontinuities
at which the derivative is infinite.
If the function wn (x) can be pre-specified in such a way
that wn at x where f (n) (x) is infinite, then the
oversmoothing can be avoided. In this way, wn (x) act as
continuity-controllers [6]. It is further suggested that wn (x)
may be discontinuous and not pre-specified [9]. For exam-
ple, by regarding wn (x) as unknown functions, one could
solve these unknowns using variational methods. But how
well wn (x) can thus be derived remains unclear. The introduction
of line processes [1], [3] or weak continuity constraints
[10] provides a solution to this problem.
A.2 Line Process Model and Its Approximations
The LP model assumes piecewise smoothness whereby
the smoothness constraint is switched off at points where
the magnitude of the signal derivative exceeds certain
threshold. It is defined on a lattice rather than on a continuous
domain. Quantize the continuous interval [a; b] into
uniformly spaced points x 1 ; :::; xm so that f
Introduce a set of binary line
process variables l i 2 f0; 1g into the smoothness term. If
also takes on a value in f0; 1g, then l i
is related it by l . The on-state, l of the line
process variable indicates that a discontinuity is detected
between neighboring points the off-state, l
indicates that the signal between the two points is contin-
uous. Each turn-on of a line process variable is penalized
by a quantity -ff. These give the LP regularizer
l i (6)
The energy for the LP model is
l i
This is the weak string model and its extension to 2D the
weak membrane model [10]. Eq.(7) corresponds to the energy
in the posterior distribution of the surface field f and
LI: ON DISCONTINUITY-ADAPTIVE SMOOTHNESS PRIORS IN COMPUTER VISION 3
the line process field l, the distribution being of the Gibbs
e \GammaU (f;l j d) (8)
where Z is a normalizing constant called the partition function

The line process variables are determined as follows: If
then it is cheaper to pay the price [f
it is more economical to
turn on the variable l to insert a discontinuity with the
cost of ff. This is an interpretation of the LP model based
on the concept of the weak continuity constraint introduced
by Blake [5] for edge labeling. An earlier idea of weak
constraints can be found in Hinton's thesis work [25]. In
the LP model, the interaction is piecewise constant (1 or
and the smoothing strength at i is either proportional to
see the next section. The concept of
discontinuities can be extended to model-based recognition
of overlapping objects [26], [27]. There, the relational bond
between any two features in the scene should be broken if
the features are ascribed to two different objects.
Finding f   2 IR m and l   2 f0; 1g m such that U(f; l j d)
is minimized is a mixture of real and combinatorial op-
timization. Algorithms for this can be classified as cate-
gories: stochastic [1], [3] and deterministic [19], [20], [10],
[17]. Some annealing techniques are often combined into
them to obtain global solutions.
In stochastic approaches, f and l are updated according
to some probability distribution parameterized by a
temperature parameter. For example, Geman and Geman
propose to use simulated annealing with the Gibbs sampler
[1] to find the global MAP solution. Marroquin [3] minimizes
the energy by a stochastic update in l together with
a deterministic update in f using gradient descent.
Deterministic approaches often use some classical gradient
based methods. Before these can be applied, the
combinatorial minimization problem has to be converted
into one of real minimization. By eliminating the line pro-
cess, Blake and Zisserman [10] convert the previous minimization
problem into one which minimizes the following
function containing only real variables
where the truncated quadratic potential function
shall be referred to as the line process potential function.
Blake and Zisserman introduce a parameter p into g ff (j) to
control the convexity of E, obtaining g (p)
ff (j). The parameter
p varies from 1 to 0, which corresponds to the variation
from a convex approximation of the function to its original
form.
Koch et al. [19] and Yuille [20] perform the conversion
using the Hopfield approach [28]. Continuous variables - l i
in the range [0; 1] are introduced to replace the binary line
process variables l i in f0; 1g. Each - l i is related to an internal
variable by a sigmoid function
with - ? as the parameter whereby lim -!0
. The
energy with this treatment is
It is shown that at stationary points where
there are v hence the approximated
line process variables [20]
This gives the effective potential function as
As the temperature decreases toward zero, - l i approaches
Geiger and Girosi [17] approximate the line process using
mean field theory. They introduce a parameter fi into (8),
giving an approximated posterior probability
Using the saddle point approximation method, they derive
mean field equations which yield the approximated line process
variables which are identical to (13). The solution is
found in the limit when fi !1.
proposes a continuous adaptive regularizer model.
There, the smoothness constraint is applied without the
switch-off as the LP model. Its effect is decreased as the
derivative magnitude becomes larger and is completely off
only at the true discontinuities where the derivative is infi-
nite. This is an earlier form of the DA model in this work.
B. Other Regularization Models
Grimson and Pavlidis [7] propose an approach in which
the degree of interaction between pixels across edges is adjusted
in order to detect discontinuities. Lee and Pavlidis
[11] investigate a class of smoothing splines which are piece-wise
polynomials. Errors of fit are measured after each
successive regularization and used to determine whether
discontinuities should be inserted. This process iterates
until convergence is reached. Besl et al. [29] propose a
smoothing window operator to prevent smoothing across
discontinuities based on robust statistics. Liu and Harris
4 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-17, NO.6, PP.576-586, JUNE 1995
[30] develop, based on a previous work [31], a computational
network in which surface reconstruction, discontinuity
detection and estimation of first and second derivatives
are performed cooperatively.
Mumford and Shah [8] define an energy on a continuous
domain
a
Z a i+1
a i
where - and ff are constants, k is the number of discontinuities
and the sequence a = a
indicates the locations of discontinuities. The minimal
solution (f   ; fa i g   ; k   ) is found by minimizing over each
value of the integer k, every sequence fa k g, and every
function f(x) continuously differentiable on each interval
a . The minimization over k is a hard problem.
Using the minimal description length principle, Leclerc
[13] presents the following function for restoration of piece-wise
constant image f from noisy data d
where ffi(\Delta) 2 f0; 1g is the Kronecker delta function. To
minimize the function, he approximates the delta function
with the exponential function parameterized by -
and approaches the solution by continuation in - toward 0.
III. The Discontinuity Adaptive MRF Model
By analyzing the smoothing mechanism in terms of the
Euler equation, it will be clear that major difference between
different models lies in their way of controlling the
interaction between neighboring points and adjusting the
smoothing strength. The DA model is defined based on
the principle that wherever a discontinuity occurs, the interaction
should diminish.
A. Defining the DA Model
We focus on the models which involve only the first order
derivative and consider the general string model
a
where
Solutions f   minimizing U(f j d) must satisfy the following
associated Euler-Lagrange differential equation or simply
the Euler equation [32]
dx
with the boundary conditions
where f a and f b are prescribed constants. In the following
discussion of solutions to the differential equation, the
following assumptions are made: Both -(x) and d(x) are
continuous and f(x) is continuously differentiable 1 . Writing
the Euler equation out yields
d
dx
A potential function g is usually chosen to be (a) even
such that (b) the derivative of g can be
expressed as the following form
where h is called an interaction function, . Obviously, h
thus defined is also even. With these assumptions, the
Euler equation can be expressed as
d
dx
The magnitude jg 0 (f 0 relates to the strength
with which a regularizer performs smoothing; and h(f 0 (x))
determines the interaction between neighboring pixels.
necessary condition for any regularization
model to be adaptive to discontinuities is
lim
where C 2 [0; 1) is a constant. The above condition with
prohibits smoothing at discontinuities where
limited (bounded)
smoothing. In any case, however, the interaction h(j) must
be small for large jjj and approaches 0 as jjj goes to 1.
This is an important guideline for selecting g and h for the
purpose of the adaptation.
Definition 1. An adaptive interaction function (AIF) h fl
parameterized by fl (? 0) is a function that satisfies:
The class of AIFs, denoted by IHI fl , is defined as the collection
of all such h fl . 2
The continuity requirement (i) guarantees the twice differentiability
of the integrand u(f j d) in (20) with respect
1 In this work, the continuity of -(x) and d(x) and differentiability of
are assumed for the variational problems defined on continuous
domains [32], [33]. However, they are not necessary for discrete problems
where [a; b] is quantized into discrete points. For example, in
the discrete case, -(x i ) is allowed to take a value in f1,0g, indicating
whether datum d(x i ) is available or not.
LI: ON DISCONTINUITY-ADAPTIVE SMOOTHNESS PRIORS IN COMPUTER VISION 5


I
Four choices of AIFs, the corresponding APFs and bands.
AIF APF Band
to f 0 , a condition for the solution f to exist [32]. How-
ever, this can be relaxed to h fl 2 C 0 for discrete problems.
The evenness of (ii) is usually assumed for spatially unbiased
smoothing. The positive definiteness of (iii) keeps
the interaction positive such that the sign of jh fl (j) will
not be altered by h fl (j). The monotony of (iv) leads to
decreasing interaction as the magnitude of the derivative
increases. The bounded asymptote property of (v) provides
the adaptive discontinuity control as stated earlier. Other
properties
the last meaning zero interaction at dis-
continuities. The above definition characterizes the properties
AIFs should possess rather than instantiates some
particular functions. Therefore, the following definition of
the DA model is rather broad.
Definition 2. The DA solution f is defined by the Euler
equation (25) constrained by
The DA solution is C 1 continuous 2 . Therefore the DA
solution and its derivative never have discontinuities in
them. The DA overcomes oversmoothing by allowing its
solution to be steep, but C 1 continuous, at point x where
data d(x) is steep. With every possible data configuration
d, every f 2 C 1 is possible.
There are two reasons for defining the DA in terms of the
constrained Euler equation: First, it captures the essence of
the DA problem; the DA model should be defined based on
the h fl therein. Second, some satisfying h fl may not have
their corresponding g fl (and hence the energy) in closed-
form, where g fl is defined below.
Definition 3. The adaptive potential function (APF) corresponding
to an h fl 2 IHI fl is defined by
is called an adaptive
string. The following are some properties of
sically, g fl is one order higher than h fl in continuity; it
2 See [33] for a comprehensive discussion about the continuity of
solutions of the class of problems to which the DA belongs.
s
s
(1) (2) (3) (4)
Fig. 1. The qualitative shapes of the four DA functions.
is even, g fl its derivative function is odd,
however, it is not necessary for g fl (1)
to be bounded. Furthermore, g fl is strictly monotonically
increasing as jjj increases because g fl
means larger jjj
leads to larger penalty g fl (j). It conforms to the original
spirit of the standard quadratic regularizers determined by
q . The line process potential function g ff does not have
such a property: Its penalty is fixed and does not increase
as jjj increases beyond
ff. The idea that large values of
are equally penalized is questionable [14].
In practice, it is not always necessary to know the explicit
definition of g fl . The most important factors are the
Euler equation and the constraining function h fl . Nonethe-
less, knowing g fl is helpful for analyzing the convexity of
E(f ).
For a given g fl (j), there exists a region of j within which
the smoothing strength jg 0
increases monotonically
as jjj increases and the function g fl is convex:
(b
The region B fl is referred to as the band. The lower and
upper bounds b L ; b H correspond to the two extrema of
(j), which can be obtained by solving g 00
we have b l = \Gammab H when g is even. When b L
thus g fl (j) is strictly convex. For h fl defined
with C ? 0, the bounds are b

Table

I instantiates four possible choices of AIFs, the
corresponding APFs and the bands. Fig.1 shows their qualitative
shapes (a trivial constant may be added to g fl (j)).
Fig.2 gives a graphical comparison of the g(j)'s for the
quadratic, the LP model and the first three DA models
listed in Table I and their derivatives, g 0
The fourth AIF, h allows bounded but
non-zero smoothing at discontinuities: lim j!1 jh 4fl
fl. It is interesting because g 00
all j (except at and leads to strictly convex mini-
mization. In fact, a positive number C in (27) leads to a
convex AIF and hence energy function 3 . The convex subset
of models for discontinuity adaptivity and M estimation is
3 This is due to the following theorem: If g(\Delta) is convex on IR, a real-valued
energy function
convex w.r.t. f for all f 2 C 1 and fixed -; d 2 C 1 .
6 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-17, NO.6, PP.576-586, JUNE 1995
Quadratic
Line-Process
Fig. 2. Comparison of different potential functions g(j) (top) and
their first derivatives g 0
appealing because they have some inherent advantages over
nonconvex models, both in stability and computational efficiency

Fig.2 also helps us visualize how the DA performs
smoothing. Like the quadratic g q the DA allows
the smoothing strength jg 0
(j)j to increase monotonically
as j increases within the band B fl . Outside the band, the
smoothing decreases as j increases and becomes zero as
This differs from the quadratic regularizer which
allows boundless smoothing when j !1. Also unlike the
LP model which shuts down smoothing abruptly just beyond
its band B
ff), the DA decreases smoothing
continuously towards zero.
B. Relations with Previous Models
The first three instantiated models behave in a similar
way to the quadratic prior model when
noticed by [15]. This can be understood by
looking at the power series expansion of g fl (j), g fl
are constants with c 1 ? 0 (the
expansion can also involve a trivial additive constant c 0 ).
Thus in this situation of sufficiently small j 2 =fl, the adaptive
model inherits the convexity of the quadratic model.
The interaction function h also well explains the differences
between various regularizers. For the quadratic
regularizer, the interaction is constant everywhere
and the smoothing strength is proportional to jjj. This
is why the quadratic regularizer leads to oversmoothing at
discontinuities where j is infinite. In the LP model, the
interaction is piecewise constant
Obviously, it inhibits oversmoothing by switching off
smoothing when jjj exceeds
ff in a binary manner.
In the LP approximations using the Hopfield approach
[19], [20] and mean field theory [17], the line process variables
are approximated by (13). This approximation effectively
results in the following interaction function
As the temperature - decreases toward zero, the above
approaches h ff =2, that is,
Obviously, h ff;- (j) with nonzero
- is a member of the AIF family, i.e. h ff;- (j) 2 IHI fl and
therefore the approximated LP models are instances of the
DA model.
It is interesting to note an observation made by Geiger
and Girosi: "sometimes a finite fi solution may be more desirable
or robust" ([17], pages 406-407) where
T . They
further suggest that there is "an optimal (finite) temperature
(fi)" for the solution. An algorithm is presented in
[34] for estimating an optimal fi. The LP approximation
with finite fi ! +1 or nonzero T ? 0 is more an instance
of the DA than the LP model which they aimed to approx-
imate. It will be shown in Section III-D that the DA model
is indeed more stable than the LP model.
Anisotropic diffusion [35] is a scale-space method for
edge-preserving smoothing. Unlike fixed coefficients in the
traditional isotropic scale-space filtering [36], anisotropic
diffusion coefficients are spatially varying according to the
gradient information. A so-called biased anisotropic diffusion
[16] model is obtained if anisotropic diffusion is combined
with a closeness term. Two choices of APFs are used
in those anisotropic diffusion models: g 1fl and g 2fl .
Shulman and Herve [14] propose to use the following
Huber's robust error penalty function [37] as the adaptive
potential
similar role to ff in g ff . The above is
a convex function and has the first derivative as: g 0
for other j. Comparing g 0
fi (j) with (24), we find that
the corresponding AIF is h fi
This function allows bounded
but nonzero smoothing at discontinuities. The same function
has also been applied by [38] to curve fitting. A comparative
study on the DA model and robust statistics can
be found in [39], [27].
The approximation (18) of Leclerc's minimal length
model [13] is in effect the same as the DA with APF 1.
may be one of the best cost functions for the piece-wise
constant restoration; for more general piecewise continuous
restoration, one needs to use (18) with a nonzero - ,
which is a DA instance. Regarding the continuity property
of domains, Mumford and Shah's model [8] and Terzopou-
los' continuity-controlled regularization model [9] can also
be defined on continuous domains as the DA model.
LI: ON DISCONTINUITY-ADAPTIVE SMOOTHNESS PRIORS IN COMPUTER VISION 7
C. Discrete Data and 2D Cases
When the data d is available at discrete points: d
m, the Euler equation for the adaptive
string is
d
dx
where ffi(\Delta) is the Dirac delta function. Integrating the
above equation once yields
where c is a constant. From the above, the solution f is
determined by the following
with appropriate boundary conditions at
xm . Obviously, f in this case is piecewise C 1 continuous;
more exactly, it is composed of consecutively joined line
segments.
The adaptive string model can be extended to 2D and
higher order equivalents. The 2D adaptive membrane has
the following
@
@x
@
@y
[f y h(f y
corresponding to the Euler equation (25), where data
images. Extending
the DA to second order, one obtains the following for
the adaptive rod on 1D
d
and for the adaptive plate on 2D
@
@xy
[f xy h(f xy
@
[f yy h(f yy
D. Solution Stability
The DA solution depends continuously on its parameters
and the data whereas the LP solution does not. An informal
analysis follows. Consider an LP solution f ff obtained
with ff. The solution is a local equilibrium satisfying
In the above, h ff (f 0
0, depending on ff and the
configuration of f ff . When jf 0
ff (x)j 2 is close to ff for some
a small change \Deltaff may flip h ff over from one state to the
other. This is due to the binary non-linearity of h ff . The
flip-over leads to a significantly different solution. This can
be expressed as
lim
\Deltaff!0
where 0(x) is a function which is constantly zero in the
domain [a; b]. The variation ffif ff with respect to \Deltaff may
not be zero for some f ff and ff, which causes instability.
However, the DA solution, denoted f fl , is stable
lim
\Deltafl !0
where f fl denotes the DA solution. Conclusion on the stability
due to changes in parameter - can be drawn similarly.
The same is also true with respect to the data. Given ff
and - fixed, the solution depends on the data, i.e.
Assume a small variation ffid in the data d. The solution
f [d] must change accordingly to reach a new equilibrium
to satisfy the Euler equation. However, there
always exist possibilities that h ff (f 0
ff (x)) may flip over for
some x when jf 0
ff (x)j 2 is near the ff, resulting in an abrupt
change in the LP solution f ff . This can be represented by
lim
That is, the variation ffif ff with respect to ffid may not be
zero for some f ff and d. However, the DA model is stable
to such changes, i.e.
lim
[f
because of its continuous nature. From the analysis, it
can be concluded that the DA better regularizes ill-posed
problems than the LP.
IV. Computation of DA Solutions
A. Solving the Euler Equation
The Euler equation (25) can be treated as a boundary
value problem. It can also be solved by minimizing
the corresponding energy (19) because a minimum of the
energy is (sufficiently) a solution of the equation (With
the form of the energy needs not to be
known in order to minimize it - see below). The energy
minimization approach is chosen here.
Because both are
bounded below, so is the energy. This means that a minimal
solution, and hence a solution to the Euler equation,
exists. If h fl is chosen with the corresponding energy
E(f) is non-convex with respect to f , and the minimization
is subject local minima. See [10] for an analysis of
convexity in string and membrane models. The following
presents a discrete method for finding a local minimum.
Sample the integral interval [a; uniformly
spaced points: x b. For clarity and without
loss of generality, let us assume that the bounds a and b are
re-scaled in such a way that the point spacing is one unit 4
4 More advanced numerical methods using varying spacing, such as
[40] may be advantageous in obtaining more accurate solutions.
8 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. PAMI-17, NO.6, PP.576-586, JUNE 1995
i.e. 1. Approximate the first derivative f 0
i by
the first order backward difference f 0
Then the energy (19) is approximated by
consists of the neighbors of i. Using the gradient-descent
method, we can obtain the following updating
(f (t)
The solution f (t) takes prescribed values at the boundary
points at to meet the boundary condition
(22) and the values may be estimated from data d i near the
boundaries. With initial f (0) , the solution is in the limit
Equation (49) helps us see more of how the DA works.
The smoothing at i is due to -
1g. The contributions to smoothing
are from the two neighboring points. That from site i 0
is proportional to the product of the two factors, f
and h fl (f This is why we relate jg 0
to the strength of smoothing performed by regularizers.
On the other hand, h(j) acts as an adaptive weighting
function to control the smoothing due to the difference
Therefore we regard it as the interaction.
Two more remarks are made: First, the contributions from
the two sides, are treated separately or
non-symmetrically. Second, the sum of the contributions
to smoothing is zero if the three points are aligned when
and h(j) is even since in this situation
The updating rule for the adaptive membrane on 2D can
be easily derived. In the 2D case (40) , there is an additional
term in the y dimension. This leads to the following
energy
Letting @f i;j
@f i;j
leads to the following updating rule
(f (t)
is the
set of the four neighboring points 5 of (i; j). The updating
on 2D grid can be performed on the white and black sites
of a checkerboard alternatively to accelerate convergence.
5 With the 4-neighborhood system, the model considers derivatives
in the horizontal and vertical directions. With the 8-neighborhood
system, the regularizer also includes the diagonal derivatives weighted
by 1=
2.
Choose a convex
using Eq.(49);
set
Until (f
Fig. 3. A GNC algorithm for finding the DA solution.
There are three parameters in (49) which shall be determined
for the DA model: -, fl and -. Parameter -
is related to the convergence of the relaxation algorithm.
There is an upper bound for - for the system to be sta-
ble. There also exists an optimal value for - [10]. Optimal
choices of - for quadratic regularization may be made using
cross-validation [41], [42]. In [43], a least squares method
[44] is presented for estimating MRF clique potentials for
a LP model. Automated selection of the - and fl parameters
for the DA model is an unsolved problem. They are
currently chosen in an ad hoc way.
B. Non-Convex Minimization
The DA model with leads to non-convex dynamic
systems and direct minimization using gradient descent
only guarantees to find a local minimum. A GNC-like
algorithm can be constructed for approximating the global
solution.
Because a convex g guarantees a convex energy function
E(f ), it is useful to study the convexity of E by analyzing
the convexity of g fl . The expansion (30) illustrates that
when
the DA model behaves in a similar way to the quadratic
regularizer and hence is convex. An analysis in [10] shows
that it is sufficient to guarantee the convexity if fl is chosen
large enough to satisfy
where c   is some real number (see an analysis in [10] for
determining c   ). The exact value of c   for the convexity
depends on g and -. It can be safely assumed that c
in any case.
Therefore, a convex fl (0) can be chosen for which
This is equivalent to choosing
a fl (0) such that f (0)
is the
band introduced in the section. For APFs 1, 2 and 3, the
must be larger than 2v, 3v and v, respectively, where
The graduation from an initially convex approximation
of g fl to its target form can be implemented by continuation
in fl from a big value to the target value. A GNC-like
algorithm using this heuristic is outlined in Fig.3. Given
d, - and a target value fl target for fl, the algorithm aims to
construct a sequence ffl (t) g with fl (1) ! fl target , and thus
LI: ON DISCONTINUITY-ADAPTIVE SMOOTHNESS PRIORS IN COMPUTER VISION 9
ff (t)
g to approach the global minimum f
for which E(f   In the algorithm, ffl is a constant
for judging the convergence, and - is a factor for decreasing
1. The choice of - controls the
balance between the quality of the solution and the computational
time. In principle, fl (t) should vary continuously
to keep a good track of the global minimum. In discrete
computation, we choose 0:9 - 0:99. More rapid decrease
of fl (t) (with smaller -) is likely to lead the system to
an unfavorable local minimum. Witkin et al. [45] present
a more sophisticated scheme for decreasing fl by relating
the step to the energy change:
are constants. This seems reasonable.
C. Analog Network
The computation can be performed using an analog net-
work. Let f i be the potential of neural cell i. Let C
1=2- be the membrane capacitance and R
the membrane resistance. Let
be the conductance or synaptic efficacy between neurons i
. If the exponential
g 1fl is used, then
Let d i be the external current input to i, with d
Now (49) can be written as
@t
The above is the dynamic equation at neuron i of the net-
work. The diagram of the network circuit is shown in Fig.4.
The synaptic current from i to i 0 is
I
If the exponential g 1fl is used, then
I
A plot of current I i;i 0 versus potential difference f
was shown at the bottom of Fig.2. The voltage-controlled
nonlinear synaptic conductance T i;i 0 , characterized by the
h function defined in (27), realizes the adaptive continuity
control; the corresponding nonlinear current I i;i 0 realizes
the adaptive smoothing. The current I i;i 0 diminishes
asymptotically to zero as the potential difference between
neurons i and i 0 reaches far beyond the band B fl .
Fig.7 shows the behavior of the analog network under
component defects such as manufacturing inadequacy,
quality changes, etc. The defects are simulated by adding
\Sigma25% evenly distributed random noise into R, C, and T
in (55). The data d is shown in triangles with 50% missing
rate; the locations of the missing data, for which - are
indicated by triangles at the bottom. The noise in the data
is white Gaussian with standard deviation
r
R
f
l, s
r r
Fig. 4. Schematic diagram of the analog network circuit for the DA
model.
and (right). The interaction function is chosen to
be h 2fl
Solutions obtained with simulated
component defects are shown in dashed lines in comparison
with those obtained without such noise shown in thicker
solid lines. The ideal signal is shown in thinner solid lines.
As can be seen, there is only a little difference between
the solutions obtained with and without such noise. This
demonstrates not only the stability of the network circuit
but also the error-tolerance property of the DA model.
V. Experiments
Two experimental results are presented in the following 6 .
The first is the reconstruction of a real image of size
256 \Theta 256 (Fig.5). Here, APF 1 (g 1fl ) is used and the parameters
are empirically chosen as 2. The
result shows that the reconstructed image is much cleaner
with discontinuities well preserved.
The second experiment is the detection of step and roof
edges from a simulated noisy pyramid image of size 128 \Theta
128 (Fig.6). The detection process runs in three stages:
regularizing the input image and computing images of
first derivatives in the two directions from the regularized
image using finite difference, 2) regularizing the derivative
images and detecting steps and roofs by thresholding
the regularized derivative images. APF 2 (g 2fl ) is used and
the parameters are empirically chosen as
for the first stage of the regularization and
for the second stage.
Edges in the horizontal and vertical directions are best
detected while those in the diagonal directions are not so
well done. This is because only derivatives in the two
axes directions are considered in the DA discussed so far;
changes in the diagonal directions are largely ignored. Regularizers
using the 8-neighborhood system (see the footnote
for Eq.(51) should help improve the detection of diagonal
changes.
Results in Fig.7 show the behavior of the analog DA
network under component defects such as manufacturing
inadequacy, quality changes, etc. The defects are simulated
by adding \Sigma25% evenly distributed random noise into R,
C, and T in (55). The data d is shown in triangles with 50%
missing rate; the locations of the missing data, for which
are indicated by triangles at the bottom. The noise
6 More results can be found in Chapter 3 of [18].
Fig. 5. 3D plots of an image (left) and its reconstruction using DA
(right). The plots are drawn after sampling the images at every
4 pixels in both directions into the size of 64 \Theta 64.
in the data is white Gaussian with standard deviation
Solutions obtained with simulated
component defects are shown in dashed lines in comparison
with those obtained without such noise shown in thicker
solid lines. The ideal signal is shown in thinner solid lines.
As can be seen, there is only a little difference between
the solutions obtained with and without such noise. This
demonstrates not only the stability of the network circuit
but also the error-tolerance property of the DA model.
Fig. 6. Step and roof edges (right) detected from a pyramid image
(left). Step edges are shown in dots and roof edges in crosses.
VI. Conclusion
Through an analysis of the associated Euler equation,
a necessary condition is made explicit for MRF or regularization
models to be adaptive to discontinuities. On
this basis, the DA model is defined by the Euler equation
constrained by the class of adaptive interaction functions
(AIFs). The definition provides principles for choosing
deterministic regularizers and MRF clique potential func-
tions. It also includes many existing models as special instances

LI: ON DISCONTINUITY-ADAPTIVE SMOOTHNESS PRIORS IN COMPUTER VISION 11
Fig. 7. Stability of the DA solution under disturbances in parameters.
The DA model has its solution in C 1 and adaptively
overrides the smoothness assumption where the assumption
is not valid, without the switching-on/off of discontinuities
in the LP model. The DA solution never contains true dis-
continuities. The DA model "preserves" discontinuities by
allowing the solution to have arbitrarily large but bounded
slopes. The LP model "preserves true discontinuities" by
switching between small and unbounded (or large) slopes.
Owing to its continuous properties, the DA model possesses
some theoretical advantages over the LP model. Unlike
the LP model, it is stable to changes in parameters and
in the data. Therefore it is better than the LP model in
solving ill-posed problems. In addition, it is able to deal
with problems on a continuous domain. Furthermore, it is
better suited for analog VLSI implementation.

Acknowledgments

The author is grateful to Yihong Huang, Eric Sung, Wei
Yun Yau and Han Wang for their helpful comments.



--R

""
""
Probabilistic Solution of Inverse Problems
""
""
""
""
""
""

""
""
""
""
""
""
""
Towards 3D Vision from Range Images: An Optimisation framework and Parallel Distributed Networks
""
""
""
""
Solutions of Ill-posed Prob- lems
""
Relaxation and Its Role in Vision
""
Markov Random Field Modeling in Computer Vision
""
""
""
""
Methods of Mathematical Physics

""
""
""
Robust Statistics
""
""
""
""
""
""
""
""
--TR

--CTR
Jian-Feng Cai , Raymond H. Chan , Carmine Fiore, Minimization of a Detail-Preserving Regularization Functional for Impulse Noise Removal, Journal of Mathematical Imaging and Vision, v.29 n.1, p.79-91, September 2007
A. Tonazzini , L. Bedini, Monte Carlo Markov chain techniques for unsupervised MRF-based image denoising, Pattern Recognition Letters, v.24 n.1-3, p.55-64, January
Michele Ceccarelli, A Finite Markov Random Field approach to fast edge-preserving image recovery, Image and Vision Computing, v.25 n.6, p.792-804, June, 2007
Stan Z. Li , Han Wang , William Y. C. Soh, Robust Estimation of Rotation Angles from Image Sequences Usingthe Annealing M-Estimator, Journal of Mathematical Imaging and Vision, v.8 n.2, p.181-192, March, 1998
David W. Jacobs , Daphna Weinshall , Yoram Gdalyahu, Classification with Nonmetric Distances: Image Retrieval and Class Representation, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.6, p.583-600, June 2000
energy functions;regularization;markov random fields;computer vision;minimization;euler equation;discontinuities