Image iterative �ltering and P.D.E.’s

Frédéric Guichard, J.M.M., Bob Ryan

0-0 Figure I.1: A noisy image magni�ed to show the pixels.

0-1 Digital image : �rst the underlying image u_{S}, then a smoothed by optics version u_{f}, �nally the digital one u_{d}, discrete on a grid. It is widely assumed that the underlying “real image” u_{S} is either a measure or, for more optimistic authors, a function that has strong discontinuities. Rudin in 1987 and De Giorgi and Ambrosio in 1988 proposed independently the space BV (ℝ^{2}) of functions with bounded variation as the correct function space for modeling the images u_{S}.

A function f is in BV (ℝ^{2}) if its partial derivatives ∂f/∂x and ∂f/∂y, taken as distributions, are Radon measures with �nite total mass. BV (ℝ^{2}) looked at �rst well adapted to modelling digital images because it contains functions having step discontinuities.

I.1 Image processing

Three disciplines with different scopes and different techniques : compression, restoration, analysis.

Figure I.2: Compression. Left to right: the original image and its increasingly compressed versions. The compression factors are roughly 7, 10, and 25.

1-1 Figure I.3: Denoising. Left: an image with up to 75% of its pixels contaminated by simulated noise. Right: a denoised version.

1-2 Figure I.4: Analysis of a shape. The original scanned shape is on the left. Simpli�ed versions are to the right.

1-3

I.2 PDEs and image processing

Figure I.5: Heat equation and smoothing. The original image is on the left; the heat equation has been applied at some scale, and the resulting blurred image is on the right. k is any positive kernel such that k(x) = k(|x|) and such that k is localized in the sense that k(x) → 0 suf�ciently rapidly as |x| → ∞. If k is normalized properly and if k_{h}(x) = (1/h)k(x/h^{1/2}), then

\frac{k_{h} ∗ u_{0}(x) − u_{0}(x)}{h} → Δu_{0}(x)

as h → 0 whenever the image u_{0} is suf�ciently smooth. We write this as

k_{h} ∗ u_{0}(x) − u_{0}(x) = hΔu_{0}(x) + o(h).  (I.1)

Now let u(t, x) denote the solution of the heat equation

\frac{∂u}{∂t} = Δu, u(0, x) = u_{0}(x).

If u_{0} is suf�ciently smooth, then we can write

u(t, x) − u(0, x) = tΔu_{0}(x) + o(t).  (I.2)

Equations (I.1) and (I.2) suggest that blurring u_{0} with a kernel k_{h} for small h is equivalent to applying the heat equation to u_{0} at some small scale t. Lindenbaum, Fischer, and Bruckstein that Kovasznay and Joseph introduced in 1955 the notion that a slightly blurred image could be deblurred by subtracting a small amount of its Laplacian. Numerically, this amounts to subtracting a fraction λ of the Laplacian of the observed image from itself (later, Dennis Gabor found the optimal parameter for the method):

u_{restored} = u_{observed} − λΔu_{observed}.

Figure I.6: Kovasznay–Joseph–Gabor deblurring. Left to right: original image; three iterations of the algorithm; ten iterations of the algorithm.

Figure I.7: Kovasznay–Joseph–Gabor deblurring. This is the same deblurring experiment as in Figure I.6, but it is applied to a much more blurred image.

Rudin in 1987 and Osher and Rudin in 1990 proposed a pseudoinverse for the heat equation, where the propagation term |Du| = |(u_{x}, u_{y})| is controlled by the sign of the Laplacian:

\frac{∂u}{∂t} = −sign(Δu)|Du|.  (I.3)

A shock �lter, propagates the level lines of an image with a constant speed and in the same direction as the reverse heat equation would propagate these lines : pseudoinverse for the heat equation.

An improvement (Guichard et al.), Kramer’s �lter (seventies) not formulated as a PDE, but equivalent to a shock �lter :

\frac{∂u}{∂t} = −sign(D^{2}u(Du, Du))|Du|,  (I.4)

where the Laplacian has been replaced by

D^{2}u(Du, Du) = u_{xx}(u_{x})^{2} + 2u_{xy}u_{x}u_{y} + u_{yy}(u_{y})^{2}.

D^{2}u(Du, Du)/|Du|^{2} is the second derivative of u in the direction of its gradient Du, Haralick’s edge detector. Kramer’s equation yields a slightly better version of a shock �lter.

A third deblurring method, the Rudin–Osher–Fatemi algorithm, is de�nitely better. It poses the deblurring problem as an inverse problem. Given the observed image u_{0}, one tries to �nd a restored version u such that k ∗ u is as close as possible to u_{0} and such that the oscillation of u is nonetheless bounded. This is done by �nding u that minimizes

∫ (|Du(x)| + λ(k ∗ u(x) − u_{0}(x))^{2}) dx.  (I.5)

Figure I.8: Deblurring with shock �lters and a variational method. Left to right: blurred image; Rudin–Osher shock �lter; Kramer’s improved shock �lter; Rudin–Osher–Fatemi restoration method.

9-1 In 1983, Burt and Adelson : Laplacian pyramid. A convolution followed by subsampling. After each convolution, one keeps only the difference k_{n} ∗ u_{n} − u_{n} ≃ Δu_{n} The image is then coded by the (�nite) sequence of these differences which tend to be sparse for natural images.

Figure I.9: The Laplacian pyramid of Burt and Adelson. Left to right: the original image; the image blurred by Gaussian convolution; the difference between the original image and the blurred version, which approximates the Laplacian of the original image; the points where this Laplacian image is large.

10-1

Edge detection

Marr and Hildreth proposed computing the points where Δu crosses zero, the now-famous zero-crossings. A signi�cant improvement was made by Haralick who de�ned the boundaries, or edges, of an image as those points where |Du| attains a local maximum along the gradient lines Canny implemented Haralick’s detector Gaussian smoothing followed by computing the (edge) points where D^{2}u(Du, Du) = 0 and |Du| is above some threshold Figure I.10: Heat equation and Haralick’s edge detector. Left to right: original image; edge points found in the original image using Haralick’s detector; blurred image; edges found in the blurred image using the Haralick–Canny detector.

11-1 Perona and Malik : to smooth what needs to be smoothed, namely, the irrelevant homogeneous regions, and to enhance the boundaries. With this in mind, the diffusion should look like the heat equation when |Du| is small, but it should act like the inverse heat equation when |Du| is large.

\frac{∂u}{∂t} = div(g(|Du|)Du),  (I.6)

where g(s) = 1/(1 + λ^{2}s^{2}). It is easily checked that we have a diffusion equation when λ|Du| ≤ 1 and an inverse diffusion equation when λ|Du| > 1.

To see this, consider the second derivative of u in the direction of Du,

u_{ξξ} = D^{2}u(\frac{Du}{|Du|}, \frac{Du}{|Du|}),

and the second derivative of u in the orthogonal direction,

u_{ηη} = D^{2}u(\frac{Du^{⊥}}{|Du|} , \frac{Du^{⊥}}{|Du|} ),

where Du = (u_{x}, u_{y}) and Du^{⊥} = (−u_{y}, u_{x}). The Laplacian can be rewritten in the intrinsic coordinates (ξ, η) as Δu = u_{ξξ} + u_{ηη}. The Perona–Malik equation then becomes

\frac{∂u}{∂t} = \frac{1}{1 + λ^{2}|Du|^{2}} u_{ηη} + \frac{1 − λ^{2}|Du|^{2}}{(1 + λ^{2}|Du|^{2})^{2}} u_{ξξ}.

This model mixes the heat equation and the reverse heat equation.

Figure I.11: A Perona–Malik equation and edge detection. This is the same experiment as in Figure I.11, but here the Perona–Malik equation is used in place of the heat equation. Notice that the edge map looks slightly better in this case.

13-1

A proliferation of anisotropic diffusions...

Rudin–Osher–Fatemi model ∫ |Du(x)| dx, which is the BV norm of u. The gradient descent for ∫ |Du(x)| dx translates into the equation

\frac{∂u}{∂t} = div(\frac{Du}{|Du|}) = \frac{1}{|Du|}u_{ηη}.

Written this way, the method appears as diffusion in the direction orthogonal to the gradient, tuned by the size of the gradient. Andreu, Caselles and Mazon proved that this equation is well posed. Alvarez, Lions, et al. :

\frac{∂u}{∂t} = \frac{1}{|k ∗ Du|} |Du| div(\frac{Du}{|Du|}) = \frac{1}{|k ∗ Du|}u_{ηη},

and again the diffusion is in the direction Du^{⊥} orthogonal to the gradient. Note that the rate of diffusion depends on the average value k ∗ Du of the gradient in a neighborhood of x, whereas the direction of diffusion, Du^{⊥}(x)/|Du(x)|, depends on the value of Du(x) at x.

Kimia, Tannenbaum, and Zucker, in a more general shape-analysis framework taken from Osher-Sethian:

\frac{∂u}{∂t} = |Du| div(\frac{Du}{|Du|}) = D^{2}u(\frac{Du^{⊥}}{|Du|} , \frac{Du^{⊥}}{|Du|} ) = u_{ηη}.  (I.7)

This equation is a “pure" diffusion in the direction orthogonal to the gradient. We call this equation the curvature equation.

The Weickert equation

\frac{∂u}{∂t} = D^{2}u(v_{1}^{⊥}, v_{1}^{⊥}) = u_{v_{1}v_{1}}

uses a nonlocal estimate of the direction orthogonal to the gradient for the diffusion direction.The eigenvector corresponding to the smallest eigenvalue of k ∗ (Du ⊗ Du), where (y ⊗ y)(x) = (x · y)y. If the convolution kernel is removed, then this eigenvector is simply Du^{⊥}.

Carmona and Zhong proposed a diffusion in the direction of the eigenvector corresponding to the smallest eigenvalue of D^{2}u Sochen, Kimmel, and Malladi proposed a nondegenerate diffusion associated with a minimal surface variational formulation [?]. Their idea was to make a gradient descent for the area, ∫ \sqrt{1 + |Du(x)|^{2}} dx, of the graph of u. This leads to the diffusion equation

\frac{∂u}{∂t} = div(\frac{Du}{\sqrt{1 + |Du|^{2}}} ).

Caselles, Morel, and Sbert proposed a diffusion to interpolate between level lines

\frac{∂u}{∂t} = D^{2}u(Du, Du) = |Du|^{2}u_{ξξ}.

Among the models mentioned, only the curvature motion proposed by Kimia, Tannenbaum, and Zucker was speci�cally introduced as a shape analysis tool.

Figure I.12: Diffusion models I. Left to right: Osher, Sethian 1988: the curvature equation; Rudin, Osher, Fatemi 1992: minimization of the image’s total variation; Alvarez, Lions, Morel 1992: nonlocal variant of the preceding; Weickert 1994: nonlocal variant of the curvature equation. All of these models diffuse only in the direction orthogonal to the gradient, using a more or less local estimate of this direction.

17-1 Figure I.13: Diffusion models II. Left to right: original image; Perona–Malik equation 1987; Carmona, Zhong 1998: diffusion along the least eigenvector of D^{2}u; Sochen, Kimmel, Malladi 1998: minimization of the image graph area.

17-2 Figure I.14: Diffusion models III. Left to right: original image; quantized image (only 10 levels are kept - 3.32 bits/pixel); the quantized image reinterpolated using the Caselles–Sbert algorithm 1998. They apply a diffusion on the quantized image with values on the remaining level lines as boundary conditions.

17-3

Contrast-invariant classes. Two images u and v are said to be

(perceptually) equivalent if there is a continuous increasing function g such that v = g(u). In this case, u and v are said to belong to the same contrast-invariant class. (“Increasing” always means “strictly increasing.”)

Contrast invariance requirement. An image analysis operator T must act directly on the equivalence class. As a consequence, we ask that T (g(u)) = g(T u), which means that the image analysis operator commutes with contrast changes. (Note that contrast changes g are not limited to increasing functions.)

Matheron in 1975 formulates image analysis as set analysis, namely, the analysis of the level sets of an image. The upper level set of an image u at level λ is the set

X_{λ}u = {x | u(x) ≥ λ}.

We de�ne in exactly the same way the lower level sets by changing “≥” into “≤.” The main point to retain here is the global invariance of level sets under contrast changes. Thus, if g is a continuous increasing contrast change, then

X_{g(λ)}g(u) = X_{λ}u.

According to mathematical morphology, all of an image’s shape information is contained in its level sets (forgetting about their levels).

Figure I.15: An image and one of its level sets. On the right is level set 140 of the left image. This experiment illustrates Matheron’s thesis that the main shape information is contained in the level sets of an image. Level sets are contrast invariant.

Extrema killer (Guichard), grain �lter (Guichard-Monasse), area �lter (Vincent) de�ned by

Vincent in 1993 This image operator simply removes all connected components of upper and lower level sets with areas smaller than some �xed value in the �rst and third de�nition. More cleverly, it removes all level lines who surround an area less than some threshold in the middle de�nition.

Figure I.16: The extrema killer �lter. Left to right: original image; extrema killer applied with area threshold equal 80 pixels; 75% salt and pepper noise added to the original image; the same �lter applied.

23-1 In 1996, Caselles, Coll et al. proposed as the main objects of analysis the level lines of an image, that is, the boundaries of its level sets. Level sets and isolevel sets {x | u(x) = λ}, computed for an interpolation of the digital image.

For an interpolation method to be useful, the level lines should have certain minimal properties: They should be composed of a �nite number of recti�able Jordan curves, and they should be nested. This means that they do not cross, and thus that they form a tree by inclusion.

Kronrod 1950 shows that if the function u is continuous, then the isolevels sets {x | u(x) = λ} are nested and thus form a tree when ordered by inclusion. These isolevel sets are not necessarily curves; they are curves, however, if u has continuous �rst derivatives. Monasse proved Kronrod’s result for lower semicontinuous and upper semicontinuous functions in 2000.

Figure I.17: Level lines of an image. Level lines, de�ned as the boundaries of level sets, can be de�ned to be a nested set of Jordan curves. They provide a contrast-invariant representation of the image. On the right are the level lines at level 183 of the left image.

Lecture notes on Image iterative �ltering and P.D.E.’s by F. Guichard, J.M.M. and Bob Ryan to be downloaded chapter by chapter at :

http://www.cmla.ens-cachan.fr/Utilisateurs/morel/GMRChap01.pdf http://www.cmla.ens-cachan.fr/Utilisateurs/morel/GMRChap02.pdf ........ ....... http://www.cmla.ens-cachan.fr/Utilisateurs/morel/GMRChap19.pdf http://www.cmla.ens-cachan.fr/Utilisateurs/morel/GMRIntro.pdf http://www.cmla.ens-cachan.fr/Utilisateurs/morel/GMRBib.pdf http://www.cmla.ens-cachan.fr/Utilisateurs/morel/GMRApndxA.pdf Many (and the most important) mentioned image processing methods are accessible in the software http://www.cmla.ens-cachan.fr.

Chen, Giga, and Goto and Alvarez et al. proved that if one adds contrast invariance to the usual isotropic invariance then all multiscale image analyses should have a curvature-dependent motion of the form

\frac{∂u}{∂t} = F (curv(u), t)|Du|,  (I.8)

where F is increasing with respect to its �rst argument. This equation can be interpreted as follows: Consider a point x on a given level curve C of u at time t. Let n(x) denote the unit vector normal to C at x and let curv(x) denote its curvature. Then the preceding equation is associated with the curve motion equation

\frac{∂x}{∂t} = F (curv(x), t)n(x)

that describes how the point x moves in the direction of the normal.

Two main examples : F (curv(u), t) = curv(u), the curvature equation. F (curv(u), t) = (curv(u))^{1/3}, the af�ne morphological scale space (AMSS)

This particular one-third power form for the curvature dependence provides an important additional invariance, namely, af�ne invariance. We would like to have complete projective invariance, but a theorem proved by Alvarez et al. shows that this is impossible. The best we can have is invariance with respect to the so-called Chinese perspective, which preserves parallelism. Most of these equations, particularly when F is a power of the curvature, have a viscosity solution in the sense of Crandall and Lions. This was shown in 1995 by Ishii and Souganidis.

As we have mentioned, contrast-invariant processing can be reduced to level set processing and, �nally, to level curve processing. The equations mentioned above are indeed equivalent to curve evolution models if existence and regularity have been established. These results exist for the most important cases, namely, for F (curv(u), t) = curv(u), called curve shortening, and for F (curv(u), t) = (curv(u))^{1/3}, known as af�ne curve shortening (Sapiro and Tannenbaum). Grayson proved existence, uniqueness, and analyticity for the curve shortening equation,

\frac{∂x}{∂t} = curv(x)n(x),  (I.9)

and Angenent, Sapiro, and Tannenbaum proved the same results for the af�ne shortening equation [?],

\frac{∂x}{∂t} = (curv(x))^{\frac{1}{3}} n(x).  (I.10)

These results are very important for image analysis : why ?

Figure I.18: Experimental veri�cation of the af�ne invariance of the af�ne shortening (AMSS). The �rst panel contains three shapes, X, A(X), and A^{−1}A(X). The second panel contains S(X), SA(X), and A^{−1}SA(X). The congruence of the �rst and third shapes implies that S and A commute. In the third panel, the same procedure has been applied using the curvature equation. Here the �rst and third shapes are not congruent, which shows that the curve shortening is not af�ne invariant, as expected.

30-1 Evans and Spruck and Chen, Giga, and Goto proved in 1991 that a continuous function moves by the curvature motion with F (curv(u), t) = curv(u)) if and only if almost all of its level curves move by curve shortening. The same result is true for the af�ne invariant curve evolution (equation (I.8) with F (curv(u), t) = (curv(u))^{1/3}) and af�ne shortening.

In the case of the curvature motion, this result provides a mathematical justi�cation for the now-classic Osher–Sethian numerical method for moving fronts. They associate with some curve or surface C its signed distance function u(x) = ±d(x, C), and the curve or surface is handled indirectly as the zero isolevel set of u. Then u is evolved by, say, the curvature motion with a classic numerical difference scheme. Thus, the evolution of the curve C is dealt with ef�ciently and accurately as a by-product of the evolution of u. The point of view that we adopt is slightly different from that of Osher and Sethian. We view the image as a generalized distance function to each one of its level sets, since we are interested in all of them.

We show in Figure I.19 how the level lines are simpli�ed by evolving the image numerically using af�ne invariant curvature motion. For clarity, we display only sixteen levels of level curves. Notice that the aim here is not subsampling; we keep the same resolution. Nor is the aim restoration; the processed image is clearly worse than the original. The aim is invariant simpli�cation leading to shape recognition.

Figure I.19: The af�ne and morphological scale space (AMSS model). Left to right: original image; level lines of this image (16 levels only); original image smoothed using the AMSS equation; level lines of the third image.

33-1 Figure I.20: Curvature scale space I. Top, left to right: original sea bird image smoothed by af�ne curvature motion at calibrated scale 1; the absolute value of the curvature. Bottom, left to right: the positive part of the curvature; the negative part. Compare with Figure I.22, where the calibrated smoothing scale is 4.

33-2 Figure I.21: Curvature scale space II. Top, left to right: original sea bird image smoothed by af�ne curvature motion at calibrated scale 4; the absolute value of the curvature. Bottom, left to right: the positive part of the curvature; the negative part. Compare with Figure I.21, where the calibrated smoothing scale is 1.

33-3 A variant of the curvature equation can be used for shape detection. A well-known method of contour detection, initially proposed by Kass, Witkin, and Terzopoulos. A better method is a variant of curvature motion proposed by Caselles, Catté, Coll, and Dibos and improved simultaneously by Caselles, Kimmel, and Sapiro and Malladi, Sethian, and Vemuri. The user draws roughly the desired contour in the image, and the algorithm then �nds the best possible contour in terms of some variational criterion E. Given an original image u_{0} containing some closed contour that we wish to approximate, we start with an edge map

g(x) = \frac{1}{1 + |Du_{0}(x)|^{2}} ,

that is, a function that vanishes on the edges of the image. The user then designates the contour of interest by drawing a polygon γ_{0} roughly following the desired contour. The geodesic snake algorithm then builds a distance function v_{0} to this initial contour, so that γ_{0} is the zero level set of v_{0}. The energy to be minimized is

E(γ) = ∫_{γ} g(x(s)) ds,

where g is the edge map associated with the original image u_{0} and s denotes the parameter measuring the length along γ. The motion of the “analyzing image" v is governed by

\frac{∂v}{∂t} (x, t) = g(x)|Dv(x)|curv(v)(x) − Dv(x) · Dg(x).

Figure I.22: Active contour, or “snake.” Left to right: original image; initial contour; evolved distance function; �nal contour.

35-1 Figure I.23: A shape parser based on level lines. The two left images are of a desk and the backs of chairs viewed from different angles. In the far left panel, one level line has been put in white. In the second panel we show, also in white, the matching level lines. The match is ambiguous, as must be expected when the same object is repeated in the scene. In the two panels on the right, we display all the matching pairs of level lines in white.

35-2

Chapter 1 : The Heat Equation

I.3 Linear smoothing and the Laplacian

Consider a function u_{0}, which we interpret to be an image, de�ned on a domain Ω of ℝ^{2}, and assume that u_{0} ∈ C^{2}(Ω). If we wish to smooth u_{0}, then the simplest way to do so without favoring a particular direction is to replace u_{0}(x) with the unweighted average of the values of u_{0} in a disk D(x, h) of radius h centered at x. This means that we replace u_{0}(x) with

M_{h}u_{0}(x) = \frac{1}{πh^{2}} ∫_{D(x,h)} u_{0}(y) dy = \frac{1}{πh^{2}} ∫_{D(0,h)} u_{0}(x + y) dy. (I.11)

35-3 Although the operator M_{h} is quite simple, it exhibits important characteristics of a general linear isotropic smoothing operator. For example, it is localizable: As h becomes small, M_{h} becomes more localized, that is, M_{h}u_{0}(x) depends only on the values of u_{0}(x) in a small neighborhood of x.

Figure I.24: Local averaging algorithm. Left to right: original image; result of replacing the grey level at each pixel by the average of the grey levels over the neighboring pixels. The shape of the neighborhood is shown by the black spot displayed in the upper right-hand corner.

36-1

Our objective is to show the relation between the action of M_{h} and the action of the Laplacian, or the heat equation. We begin by expanding u_{0} around the point x using Taylor’s formula:

u_{0}(x + y) = u_{0}(x) + Du_{0}(x) · y + \frac{1}{2}D^{2}u_{0}(x)(y, y) + o(|y|^{2}). (I.12)

Now apply M_{h} to both sides of this expansion one has

\frac{1}{πh^{2}} ∫_{D(0,h)} u_{0}(x + y) dy = M_{h}u_{0}(x).

The �rst term on the right-hand side is

\frac{1}{πh^{2}} ∫_{D(0,h)} u_{0}(x) dy = u_{0}(x).

The second (linear) term on the right-hand side vanishes because

∫_{D(0,h)} x dy = ∫_{D(0,h)} y dy = 0.  (I.13)

36-1 Similarly, the cross term in the third (quadratic) term on the right-hand side vanishes:

∫_{D(0,h)} xy dy = 0.  (I.14)

This leaves the problem of computing the other quadratic terms, but since

∫_{D(0,h)} x^{2} dy = ∫_{D(0,h)} y^{2} dy,  (I.15)

we have

∫_{D(0,h)} x^{2} dy = \frac{1}{2} ∫_{D(0,h)}(x^{2} + y^{2}) dy = \frac{1}{2} ∫_{0}^{h} ∫_{0}^{2πr^{3}} dr dθ = \frac{πh^{4}}{4} .

Putting this all together, we have

M_{h}u_{0}(x) = u_{0}(x) + \frac{h^{2}}{8} Δu_{0}(x) + h^{2}ε(x, h).  (I.16)

36-2 Let u(t, x) be the solution of the heat equation with initial condition u_{0}. We shall prove :

lim_{{n→+∞ \atop nh^{2}→t}} M_{h}^{n}u_{0}(x) = u(t, x).  (I.17)

This says that the asymptotic limits of M_{h}^{n} operate on u_{0} exactly as the heat equation does; it thus establishes a relation between smoothing with M_{h} and smoothing with the Gaussian.

Proposition 1 (The Gaussian and the heat equation). For all

t > 0, the function x ↦ G_{t}(x) = (1/(4πt)^{N/2})e^{−|x|^{2}/4t} belongs to S and satis�es the heat equation

\frac{∂G_{t}}{∂t} − ΔG_{t} = 0.

Figure I.25: The Gaussian in two dimensions. 37-1

Standard extension of an image Assume that the real-valued

image u is originally de�ned on the hypercube [0, 1]^{N}. We �rst extend u to the hypercube C = [−1, 1]^{N} by symmetry across the coordinate hyperplanes: For each point x = (x_{1}, x_{2}, …, x_{N}) in C, de�ne u(x) by

u(x) = u(|x_{1}|, |x_{2}|, …, |x_{N}|).

This extension satis�es the following relation: For each x = (x_{1}, x_{2}, …, x_{N}) in C,

u(x_{1}, x_{2}, …, x_{N}) = u(ε_{1}x_{1}, ε_{1}x_{2}, …, ε_{N}x_{N}),  (I.18)

where (ε_{1}, ε_{2}, …, ε_{N}) is any vector in {−1, 1}^{N}. Next, we extend u periodically to all of ℝ^{N}. Since the period must be two in each of the variables x_{i},

u(x_{1}, x_{2}, …, x_{N}) = u(x_{1} + 2n_{1}, x_{2} + 2n_{2}, …, x_{N} + 2n_{N}), (I.19)

for all (x_{1}, x_{2}, …, x_{N}) ∈ ℝ^{N} and all (n_{1}, n_{2}, …, n_{N}) ∈ ℤ^{N}.

The purpose of these successive extensions of u to all of ℝ^{N} is to facilitate the de�nition of certain operations on u, such as convolution with smoothing kernels, and, at the same time, to preserve the continuity of u if u is continuous on [0, 1]^{N}. This method of extending a function is widely used in image processing; in particular, it is used in most compression and transmission standards. It is easy to check that the extension of u is continuous on ℝ^{N} if the original u is continuous on [0, 1]^{N}.

38-1 Figure I.26: Image extension by symmetry, followed by periodization.

38-2 Theorem[Existence and uniqueness of solutions of the heat equation] Assume that u_{0} ∈ L_{C}^{1}. Then we have the following two results:

1. For all t > 0 and x ∈ ℝ^{N}, u(t, ·) = G_{t} ∗ u_{0} satis�es the heat equation with initial value u_{0} in the sense that

\frac{∂u}{∂t} = Δu and ∫_{C} |u(t, x) − u_{0}(x)| dx → 0 as t → 0. (I.20)

Furthermore, for each t > 0, u(t, ·) ∈ C^{∞}(ℝ^{N}) ∩ L_{C}^{∞}, and for any t_{0} > 0 it is uniformly bounded for t ∈ [t_{0}, +∞), speci�cally,

sup_{x∈ℝ^{N},t≥t_{0}} |u(t, x)| ≤ C(t_{0})∥u_{0}∥_{L^{1}(C)}.  (I.21)

2. Conversely, given u_{0} ∈ L_{C}^{1}, G_{t} ∗ u_{0} is the unique solution u of (I.20) that belongs to L_{C}^{∞} for each t > 0, that is uniformly bounded on [t_{0}, +∞) for each t_{0} > 0, and that is C^{2} on (0, +∞) × ℝ^{N}.

Figure I.27: Level lines and the heat equation. Top, left to right: original 410×270 grey level image; level lines of original image for levels at multiples of 12. Bottom, left to right: original image smoothed by the heat equation (convolution with the Gaussian). The standard deviation of the Gaussian is 4, which means that its spatial range is comparable to a disk of radius 4. The image gets blurred by the convolution, which averages grey level values and removes all sharp edges. This can be appreciated on the right, where we have displayed all level lines for levels at multiples of 12. Note how some level lines on the boundaries of the image have split into parallel level lines that have drifted away from each other. The image has become smooth, but it is losing its structure.

39-1

I.4 Applications of the Heat Equation to Image Analysis

I.4.1 Convergence of iterated �lters to the heat equation

We wish to prove that iterated linear smoothing is asymptotic to applying the heat equation, or equivalently, to convolution with the Gaussian at different scales.

De�nition 1. We say that a real-valued kernel g ∈ L^{1}(ℝ^{N}) is Laplacian consistent if it satis�es the following moment conditions:

1. ∫_{ℝ^{N}} g(x) dx = 1.

2. For i = 1, 2, …, N, ∫_{ℝ^{N}} x_{i}g(x) dx = 0.

3. For each pair i, j = 1, 2, …, N, i ≠ j, ∫_{ℝ^{N}} x_{i}x_{j}g(x) dx = 0.

4. For i = 1, 2, …, N, ∫_{ℝ^{N}} x_{i}^{2}g(x) dx = σ, where σ > 0.

5. ∫_{ℝ^{N}} |x|^{3}|g(x)| dx < +∞.

39-2 We consider rescalings of a kernel g de�ned by

g_{h}(x) = \frac{1}{h^{N/2}} g(\frac{x}{h^{1/2}} )  (I.22)

for h > 0. Figure I.28: The rescalings g_{t}(x) = (1/t^{2})g(x/t) of a kernel for t=4, 3, and 2.

We denote the convolution of g with itself n times by g^{n∗}. The main result of this section concerns the behavior of g_{h}^{n∗} as n → +∞ and h → 0.

Theorem 1. If g is Laplacian consistent, then for every u ∈ L_{C}^{∞} ∩ C^{3},

g_{h} ∗ u(x) − u(x) = h\frac{σ}{2} Δu(x) + ε(h, x)  (I.23)

where |ε(h, x)| ≤ Ch^{3/2}.

Proof.

g_{h} ∗ u(x) − u(x) = ∫_{ℝ^{N}} \frac{1}{h^{N/2}} g(\frac{x}{h^{1/2}} )(u(x − y) − u(x)) dy
= ∫_{ℝ^{N}} g(z)(u(x − h^{1/2}z) − u(x)) dz.

Using Taylor’s formula with the Lagrange remainder, we have u(x−h^{1/2}z)−u(x) = −h^{1/2}Du(x)·z+\frac{h}{2}D^{2}u(x)(z, z)−\frac{1}{6}h^{3/2}D^{3}u(x−h^{1/2}θz)(z, z, z), where θ = θ(x, z, h) ∈ [0, 1]. By condition (1), ∫_{ℝ^{N}} g(z)Du(x) · z dz = 0; by conditions (3) and (4), ∫_{ℝ^{N}} g(z)D^{2}u(x)(z, z) dz = σΔu(x). Thus, g_{h}∗u(x)−u(x) = h\frac{σ}{2} Δu(x)−\frac{1}{6}h^{3/2} ∫_{ℝ^{N}} g(z)D^{3}u(x−h^{1/2}θz)(z, z, z) dz.

41-1 Theorem 2. Let g be a nonnegative Laplacian-consistent kernel with σ = 2 and de�ne g_{h} by (I.22). Write T_{h}u_{0} = g_{h} ∗ u_{0} for u_{0} ∈ L_{C}^{1}, and let u(t, ·) = G_{t} ∗ u_{0} be the solution of the heat equation (I.20). Then, for each t > 0,

T_{h}^{n}u_{0} → u(t, ·) in L^{1}(C) as n → +∞ and nh → t.  (I.24)

Figure I.29: Iterated linear smoothing converges to the heat equation. In this experiment with one-dimensional functions, it can be appreciated how fast an iterated convolution of a positive kernel converges to a Gaussian. On the left are displayed nine iterations of the convolution of the characteristic function of an interval with itself, with appropriate rescalings. On the right, the same experiment is repeated with a much more irregular kernel. The convergence is almost as fast as the �rst case.

42-1

Proof

T_{h}u(t, x) − u(t, x) = hΔu(t, x) + O(h^{3/2}),  (I.25)

Since u is a solution of the heat equation, we also have

u(t + h, x) − u(t, x) = hΔu(t, x) + O(h^{2}).  (I.26)

By subtracting (I.26) from (I.25) we see that

T_{h}u(t, x) = u(t + h, x) + O(h^{3/2}).  (I.27)

This shows that applying T_{h} to a solution of the heat equation at time t advances the solution to time t + h, plus an error term.

Proof We can apply T_{h} to both sides of (I.27) and obtain

T_{h}^{2}u(t, x) = T_{h}u(t + h, x) + O(h^{3/2}). (I.28)

If we write equation (I.27) with t + h in place of t and substitute the expression for T_{h}u(t + h, x) in equation (I.28), we have

T_{h}^{2}u(t, x) = u(t + 2h, x) + 2O(h^{3/2}). (I.29)

We can iterate this process and get the relation

T_{h}^{n}u(t, x) = u(t + nh, x) + nO(h^{3/2}). (I.30)

Take t = t_{1} and h = τ/n. Then

T_{h}^{n}u(t_{1}, x) = u(t_{1} + τ, x) + O((\frac{τ}{n})^{1/2}) (I.31)

and

lim_{{n→+∞ \atop nh=τ}} T_{h}^{n}u(t_{1}, x) = u(t_{1} + τ, x).  (I.32)

I.5 Directional averages and directional heat equations

The purpose of this section is to show how to approximate directional second derivatives using local averages. Consider a C^{2} function from ℝ^{N} into ℝ and a vector z ∈ ℝ^{N} with |z| = 1. We wish to compute the mean value of u along a segment of the line through x parallel to the vector z. To do this, we de�ne the operator T_{h}^{z}, h ∈ [−1, 1], by

T_{h}^{z}u(x) = \frac{1}{2h} ∫_{−h}^{h} u(x + sz) ds.

We use Taylor’s formula to expand u at the point x along the line through x parallel to the vector z:

u(x + sz) = u(x) + sDu(x) · z + \frac{s^{2}}{2} D^{2}u(x)(z, z) + o(s^{2}). (I.33)

By averaging both sides of (I.33) over s ∈ [−h, h], we obtain the next result.

45-1

Proposition 2.

T_{h}^{z}u(x) = \frac{1}{2h} ∫_{−h}^{h} u(x + sz) = u(x) + \frac{h^{2}}{6} D^{2}u(x)(z, z) + o(h^{2}).

This suggests that iterations of the operator T_{h}^{z} are associated with the directional heat equation

\frac{∂u}{∂t}(t, x) = \frac{1}{6}D^{2}u(t, x)(z, z)  (I.34)

This is correct under the assumption that z is �xed. We are able to approximate the directional second derivative by taking directional averages when z varies from point to point. The main choice considered in the book is z = Du^{⊥}/|Du|, where Du = (u_{x}, u_{y}) and Du^{⊥} = (−u_{y}, u_{x}).

45-2 1. Average in the direction of the gradient. By choosing z = Du/|Du|,

u_{ξξ} = \frac{1}{|Du|^{2}} D^{2}u(Du, Du) = 6 lim_{h→0} \frac{T_{h}^{Du/|Du|}u − u}{h^{2}}.

We will interpret this differential operator as Haralick’s edge detector.

2. Average in the direction orthogonal to the gradient. By choosing z = Du/|Du^{⊥}|,

\frac{1}{|Du|^{2}} D^{2}u(Du^{⊥}, Du^{⊥}) = 6 lim_{h→0} \frac{T_{h}^{Du^{⊥}/|Du|}u − u}{h^{2}}.

This differential operator appears as the second term of the curvature equation.

Exercise 1. Suppose that u ∈ C^{2}(ℝ). If Du(x) ≠ 0, prove that

u_{ξξ}(x) = lim_{h→0} \frac{1}{h^{2}} ( max_{y∈D(0,h)} u(x+y)+ min_{y∈D(0,h)} u(x+y)−2u(x)), (I.35)

where D(0, h) is the disk of radius h centered at the origin. 46-1

I.6 Edge detection and linear scale space

The general analysis framework in which an image is associated with smoothed versions of itself at several scales is called scale space. We are �rst going to discuss some of the ways linear smoothing has been used in image analysis. We will then indicate the problems associated with the linear theory.

I.6.1 The edge detection doctrine

Algorithm 1 (Edge detection: Marr–Hildreth zero-crossings).

1. Create the multiscale images u(t, ·) = G_{t} ∗ u_{0} for increasing values of t.

2. At each scale t, compute all the points where Δu changes sign and Du ≠ 0. These points are called zero-crossings of the Laplacian, or simply zero-crossings.

3. (Optional) Eliminate the zero-crossings where the gradient is small.

Figure I.30: A three-dimensional representation of the Laplacian of the Gaussian. This convolution kernel, which is a wavelet, is used to estimate the Laplacian of an image at different scales of linear smoothing.

48-1 Figure I.31: Zero-crossings of the Laplacian at different scales (David Marr). To extract more global structure, the image is convolved with Gaussians whose variances are powers of two. One computes the Laplacian of the smoothed image and displays the lines along which this Laplacian changes sign: the zerocrossings of the Laplacian. According to Marr, these zero-crossings represent the “raw primal sketch" of the image, or the essential information on which further vision algorithms should be based. Above, left to right: the results of smoothing and the associated Gaussian kernels at scales 1, 2, and 4. Below, left to right: the zero-crossings of the Laplacian and the corresponding kernels, which are the Laplacians of the Gaussians used above.

48-2 Haralick gives up linearity and de�ne edge points as those points where the gradient has a local maximum in the direction of the gradient. In other words, an edge point x satis�es g′(0) = 0, where g(t) = |Du(x + tDu(x)|/|Du(x)|. This implies that D^{2}u(x)(Du(x), Du(x)) = 0.

Exercise 2. Check that g′(0) = 0 implies that D^{2}u(x)(Du(x), Du(x)) = 0.

Algorithm 2 (Edge detection: The Haralick–Canny edge detector).

1. As before, create the multiscale images u(t, ·) = G_{t} ∗ u_{0} for increasing values of t.

2. At each scale t, �nd all points x where Du(x) ≠ 0 and D^{2}u(x)(z, z) crosses zero, z = Du/|Du|. At such points, the function s ↦ u(x + sz) changes from concave to convex, or conversely, as s passes through zero.

3. At each scale t, �x a threshold θ(t) and retain as edge points at scale t only those points found above that satisfy |Du(x)| > θ(t).

In practice, edges are computed for a �nite number of dyadic scales, t = 2^{n}, n ∈ ℤ.

Figure I.32: Canny’s edge detector. Left column: result of the Canny �lter without the threshold on the gradient. Middle column: result with a visually “optimal” scale and an image-dependent threshold (from top to bottom: 15, 0.5, 0.6). Right column: result with a �xed gradient threshold equal to 0.5. Note that this edge detection “theory” depends on no fewer than three parameters that must be �xed by the user.

50-1 Figure I.33: Zero-crossings of the Laplacian of a synthetic image. Left to right: the original image; the image linearly smoothed by convolution with a Gaussian; the sign of the Laplacian of the �ltered image (the gray color corresponds to values close to 0, black to clear-cut negative values, white to clear-cut positive values); the zero-crossings of the Laplacian. This experiment clearly shows a drawback of the Laplacian as edge detector.

50-2 Figure I.34: The heat equation creates structure. This experiment shows that linear scale space can create new structures and thus increase the complexity of an image. Left to right: The original synthetic image (a) contains three gray levels. The black disk is a regional and absolute minimum. The “white” ring around the black disk is a regional and absolute maximum. The outer gray ring has a gray value between the other two and is a regional minimum. The second image (b) shows what happens when (a) is smoothed with the heat equation: New extrema have appeared. Image (c) illustrates the action on (a) of a contrast-invariant local �lter, the iterated median �lter.

50-3 Figure I.35: Violation of the inclusion by the linear scale space. Top, left: an image that contains a black disk enclosed by a white disk. Top, right: At a certain scale, the black and white circles mix together. Bottom level lines of the above images: the inclusion is no longer preserved.

50-4

Contrast invariance

If image analysis is to be robust, it must be invariant under changes in lighting that produce contrast changes. It must also be invariant under the nonlinear response of the sensors used to capture an image. These, and perhaps other, contrast changes are modeled by g strictly increasing. A problem with linear theory is that linear smoothing, that is, convolution, is not generally contrast invariant:

g(k ∗ u) ≠ k ∗ (g(u)).

In the same way, the operator T_{t} that maps u_{0} into the solution of the heat equation, u(t, ·) is not generally contrast invariant. In fact, if g is C^{2}, then

\frac{∂(g(u))}{∂t} = g′(u)\frac{∂u}{∂t}

and

Δ(g(u)) = g′(u)Δu + g′′(u)|Du|^{2}.

Of course, if g is linear, then g(u) satis�es the heat equation if u does.

Exercise 3. Construct simple functions u, g, and k such that g(k ∗ u) ≠ k ∗ (g(u)).

The heat equation : fours ways to smooth a shape

53

I.7 Dynamic shape

In 1986, Koenderink and van Doorn de�ned a shape in ℝ^{N} to be a closed subset X of ℝ^{N} and proposed to smooth the shape by applying the heat equation ∂u/∂t − Δu = 0 directly to 1_{X}, the characteristic function of X. Of course, the solution G_{t} ∗ 1_{X} is not a characteristic function. The authors de�ned the evolved shape at scale t to be

X_{t} = {x | u(t, x) ≥ 1/2}.

The value 1/2 is chosen so the following simple requirement is satis�ed: Suppose that X is the half-plane X = {(x, y) | (x, y) ∈ ℝ^{2}, x ≥ 0}. The requirement is that

X = X_{t} = {(x, y) | G_{t} ∗ 1_{X}(x, y) ≥ λ},

and this is true only if λ = 1/2.

Figure I.36: Nonlocal interactions in the dynamic shape method and the creation of singularity (Alan Yuille : “The creation of structure in dynamic shape"). Left to right: Two close disks interact as the scale increases. This creates a new, qualitatively different, shape. The change of topology, at the scale where the two disks merge into one shape, also entails the appearance of a singularity (a cusp) on the shape(s) boundaries.

54-1 Figure I.37: Nonlocal behavior with the dynamic shape method. Smoothing of two irregular shapes by the dynamic shape method (Koenderink–van Doorn). Top left: initial image, made of two irregular shapes. From left to right, top to bottom: dynamic shape smoothing with increasing Gaussian variance. Notice how the shapes merge more and more. We do not have a separate analysis of each shape but rather a “joint analysis” of the two shapes. The way the shapes merge is of course sensitive to the initial distance between the shapes.

54-2

I.8 Curve evolution using the heat equation

We consider shapes in ℝ^{2} whose boundaries can be represented by a �nite number of simple closed recti�able Jordan curves. We will focus on smoothing one of these Jordan curves, which we call C_{0}. We assume that C_{0} is parameterized by s ∈ [0, L], where L is the length of the curve. Thus, C_{0} is represented as x_{0}(s) = (x(s), y(s)), where s is the length of the curve between x_{0}(0) and x_{0}(s).

At �rst glance, it might seem reasonable to smooth C_{0} by smoothing the coordinate functions x and y separately.

\frac{∂x}{∂t}(t, s) = \frac{∂^{2}x}{∂s^{2}}(t, s), x(0, s) = x_{0}(s).  (I.36)

By the maximum principle, the evolved curves C_{t} remain in the rectangle that held C_{0} and x(t, ·) and y(t, ·) are C^{∞}. There are, however, at least two reasons that argue against smoothing curves this way:

1. When t > 0, s is no longer a length parameter for the evolved curve C_{t}.

2. Although x(t, ·) and y(t, ·) are C^{∞} for t > 0, this does not imply that the curves C_{t} have similar smoothness properties. It is possible for an evolved curve to cross itself and to develop singularities.

Figure I.38: Curve evolution by the heat equation. The coordinates of the curves are parameterized by the arc length and then smoothed as real functions of the length using the heat equation. From A to D: the coordinates are smoothed with an increasing scale. Each coordinate function therefore is C^{∞}; the evolving curve can, however, develop self-crossings (as in C) or singularities (as in D).

56-1 Exercise 4. Construct a C^{∞} mapping f : [0, 1] → ℝ^{2} such that the image of [0, 1] is a square. This shows that a curve can have a C^{∞} parameterization without being smooth.

The fact that curves with corners can be described by a C^{2} parameterization is used in the seminal paper by Kass, Witkin and Terzopoulos on “snakes, active contour models" (1988).

I.9 Restoring locality and causality

I.9.1 Localizing the dynamic shape method Algorithm 3 (The Merriman–Bence–Osher algorithm).

1. Convolve the characteristic function of the initial shape X_{0} with G_{h}, where h is small.

2. De�ne X_{1} = {x | G_{h} ∗ 1_{X_{0}} ≥ 1/2}.

3. Set X_{0} = X_{1} and go back to (1).

The dynamic shape method itself is an example of a median �lter.

Figure I.39: The Bence–Merriman–Osher shape smoothing method is a localized and iterated version of the dynamic shape method. A convolution of the binary image with small-sized Gaussians is alternated with mid-level thresholding. It uses the same initial data (top, left) as in Figure I.37. From left to right, top to bottom: smoothing with increasing scales. Notice that the shapes remain separate. In fact, their is no interaction between the evolving shapes. Each one evolves as if the other did not exist.

58-1

I.9.2 Renormalized heat equation for curves

In 1992, Mackworth and Mokhtarian :

Algorithm 4 (Renormalized heat equation for curves).

1. Convolve the initial curve x_{0}, parameterized by its length parameter s_{0} ∈ [0, L_{0}], with the Gaussian G_{h}, where h is small.

2. Let L_{n} denote the length of the curve x_{n} obtained after n iterations and let s_{n} denote its length parameter. For n ≥ 1, write x̃_{n+1}(s_{n}) = G_{h} ∗ x_{n}(s_{n}). Then reparameterize x̃_{n+1} by its length parameter s_{n+1} ∈ [0, L_{n+1}], and denote it by x_{n+1}.

3. Iterate.

Theorem 3. Let x be a C^{2} curve parameterized by its length parameter s ∈ [0, L]. Then for small h,

G_{h} ∗ x(s) − x(s) = h\frac{∂^{2}x}{∂s^{2}} + o(h)(= h Curv(x) + o(h). (I.37)

It is reasonable to conjecture that, in the asymptotic limit, the renormalized shape scale space will yield the solution of following evolution equation:

\frac{∂x}{∂t} = \frac{∂^{2}x}{∂s^{2}},  (I.38)

where x_{0} = x(0, ·). It is important to note that (I.38) is not the heat equation. Indeed, from the Algorithm we see that s must denote the length parameter of the evolved curve x(t, ·) at time t. Since ∂^{2}x/∂s^{2} is the curvature vector Curv(x), this is easier to remember if we write (I.38) as

\frac{∂x}{∂t} = Curv(x).

60-1 Figure I.40: Curve evolution by the renormalized heat equation (Mackworth– Mokhtarian). After each smoothing step, the coordinates of the curve are reparameterized by the arc length of the smoothed curve. From A to D: the curve is smoothed with an increasing scale. Note that, in contrast with the linear heat equation (Figure I.38), the evolving curve shows no singularities and does not cross itself. 60-2

I.10 From an image to its level sets and back

\overline{ℝ} = ℝ ∪ \{−∞, +∞\}.

De�nition 2. Let u be a function de�ned on ℝ^{N} with values in \overline{ℝ}. The level set of u at level λ is denoted by X_{λ}u and de�ned by

X_{λ}u = {x | u(x) ≥ λ}.

We can reconstruct u from its level sets X_{λ}u by the formula

u(x) = sup\{λ | x ∈ X_{λ}u, λ ∈ \overline{ℝ}\}.

If a contrast change g : \overline{ℝ} → \overline{ℝ} is continuous and increasing. Then v = g(u) and u have the same level sets globally, that is, every level set of u is a level set of v, and conversely.

Figure I.41: Level sets of a digital image. Left to right, top to bottom: We �rst show an image with range of gray levels from 0 to 255. Then we show eight level sets in decreasing order from λ = 225 to λ = 50, where the grayscale step is 25. Notice how essential features of the shapes are contained in the boundaries of level sets, the level lines. Each level set (which appears as white) is contained in the next one, as guaranteed by Proposition 3.2.

61-1 Proposition 3. Let X_{λ}, be the level sets X_{λ}u of a function u. Then 1. If λ > μ, then X_{λ} ⊂ X_{μ}, and X_{−∞} = ℝ^{N}.

2. X_{λ} = ⋂_{μ<λ} X_{μ} for every λ > −∞. Conversely, if (X_{λ})_{λ∈\overline{ℝ}} is a family of subsets of ℝ^{N} that satis�es 1), and 2), then the level sets of the function u de�ned by

u(x) = sup\{λ | x ∈ X_{μ}, λ ∈ \overline{ℝ}\}  (I.39)

satisfy

X_{λ} = X_{λ}u  (I.40)

Proof X_{λ} = X_{λ}u.

(1) If λ > μ, then X_{λ} ⊂ X_{μ}, and X_{−∞} = ℝ^{N}.

(2) X_{λ} = ⋂_{μ<λ} X_{μ} for every λ > −∞.

Relation (1) follows directly from the de�nition of X_{λ}u. Relation (2) follows from u(x) ≥ λ iff u(x) ≥ μ for every μ < λ. Conversely, assume that a family of subsets (X_{λ})_{λ∈\overline{ℝ}} satis�es (1) and (2). Let u(x) = sup\{λ | x ∈ X_{μ}, λ ∈ \overline{ℝ}\} Suppose that x ∈ X_{λ}. By de�nition of u, u(x) ≥ λ, and hence x ∈ X_{λ}u. Thus, X_{λ} ⊂ X_{λ}u.

Assume now that x ∈ X_{λ}u. Then u(x) = sup{ν | x ∈ X_{ν}} ≥ λ. Consider any μ < λ. Then there exists a μ′ such that μ < μ′ ≤ sup{ν | x ∈ X_{ν}} and x ∈ X_{μ′}. It follows from (1) that x ∈ X_{μ}. Since μ < λ was arbitrary, we get x ∈ ⋂_{μ<λ} X_{μ} = X_{λ}.

De�nition 3. 1. If λ > μ, then X_{λ} ⊂ X_{μ}, and X_{−∞} = ℝ^{N}.

2. X_{λ} = ⋂_{μ<λ} X_{μ} for every λ > −∞.

A family of subsets (X_{λ})_{λ∈\overline{ℝ}} of ℝ^{N} that satis�es condition 1. is called nonincreasing; it is called closed nonincreasing if it also satis�es condition 2..

Figure I.42: A simple synthetic image and all of its level sets (in white) with decreasing levels, from left to right and from top to bottom.

64-1 Figure I.43: Reconstruction of an image from its level sets: an illustration of Proposition 3.2. We use four different subsets of the image’s level sets to give four reconstructions. Top, left: all level sets; top, right: all level sets whose gray level is a multiple of 8; bottom, left: multiples of 16; bottom, right: multiples of 32. Notice the relative stability of the image shape content under these drastic quantizations of the gray levels.

64-2

I.10.1 Functions and level sets de�ned almost everywhere

Lemma 1. Assume that (X_{λ})_{λ∈ℝ} is a nonincreasing family of measurable sets. Then,

X_{λ} = ⋂_{μ<λ} X_{μ}  (I.41)

almost everywhere, with the possible exception of at most a countable number of λ′s,

Proof Let h be the Gauss function G_{1}. Write m(X) = ∫_{X} h(x) dx. The function λ ↦ m(X_{λ}) is nonincreasing. Thus, it is continuous except possibly for a countable number of points where it “jumps.” Call this exceptional set Λ. Then lim_{μ→λ} m(X_{μ}) = m(X_{λ}) for λ ∈ ℝ \ Λ. Thus m(X_{μ}) ↓ m(⋂_{μ<λ} X_{μ}) as μ ↑ λ. Thus, m(⋂_{μ<λ} X_{μ}) = m(X_{λ}). Since X_{λ} ⊂ ⋂_{μ<λ} X_{μ}, this means that m(⋂_{μ<λ} X_{μ} \ X_{λ}) = 0, which implies (I.41).

There is a one-to-one correspondence between the points λ ∈ \overline{ℝ} where the function λ ↦ m(X_{λ}) is discontinuous and the sets {x | u(x) = λ} that have positive measure.

Corollary 1. Let (X_{λ})_{λ∈\overline{ℝ}} be a nondecreasing family in ℝ^{N} and de�ne u : ℝ^{N} → \overline{ℝ} by

u(x) = sup{λ | x ∈ X_{λ}}.

Then X_{λ} = X_{λ}u almost everywhere, except possibly for at most a countable number of λ′s.

65-1

I.11 Contrast changes and level sets

De�nition 4. Any nondecreasing function g : ℝ → ℝ will be called a contrast change.

De�nition 5. The pseudoinverse of any contrast change g : \overline{ℝ} → \overline{ℝ} is de�ned by

g^{(−1)}(λ) = inf{r | g(r) ≥ λ}.

Exercise 5. Compute g^{(−1)} for the following functions.

1. g(s) = max(0, s).

2. g(s) = 1 if s ≥ 0, and g(s) = 0 if s < 0.

Figure I.44: The histogram of the image Bird. For each i ∈ {0, 1, …, 255}, we display (above, right) the function h(i) = Card{x | u(x) = i}. The function below is the is given by g(i) = Card{x | u(x) ≤ i}, an integral of h. It provides an indication about the overall contrast of the image and about the contrast change imposed by the sensors. The inverse function g^{−1} can be used as a contrast change to create an image g^{−1}(u) with a �at histogram.

66-1 Lemma 2. Suppose that g : \overline{ℝ} → \overline{ℝ} is a contrast change. Then for λ ∈ \overline{ℝ}, we have the following two properties:

If λ ≤ g(g^{(−1)}(λ)), then g(s) ≥ λ if and only if s ≥ g^{(−1)}(λ). (I.42)
If λ > g(g^{(−1)}(λ)), then g(s) ≥ λ if and only if s > g^{(−1)}(λ). (I.43)

Theorem 4. Let u be a real-valued function de�ned on ℝ^{N} and let g be a contrast change. Then for λ ∈ \overline{ℝ}, the level set X_{λ}g(u) satis�es one of the following relations: Either

X_{λ}g(u) = {x | u(x) ≥ g^{(−1)}(λ)} = X_{g^{(−1)}(λ)}u,  (I.44)

or

X_{λ}g(u) = {x | u(x) > g^{(−1)}(λ)} = ⋃_{μ>g^{(−1)}(λ)} X_{μ}u.  (I.45)

Figure I.45: Contrast changes and an equivalence class of images. The three images have exactly the same level sets and level lines, but their level sets are mapped onto three different gray-level scales. The graphs on the right are the graphs of the contrast changes u ↦ g(u) that have been applied to the initial gray levels. The �rst one is concave; it enhances the darker parts of the image. The second one is the identity; it leaves the image unaltered. The third one is convex; it enhances the brighter parts of the image. Software allows one to manipulate the contrast of an image to obtain the best visualization. From the image analysis viewpoint, image data should be considered as an equivalence class under all possible contrast changes.

67-1 Remark Let u : ℝ^{N} → ℝ be measurable and let g be a contrast change. Then for λ ∈ ℝ and ε > 0, there exists a μ ∈ ℝ such that X_{μ}u ⊂ X_{λ}g(u) and meas(X_{λ}g(u) \ X_{μ}u) < ε.

Figure I.46: The two images (left) have the same set of level sets. The contrast change that maps the upper image onto the lower image is displayed on the right. It corresponds to one of the possible g functions whose existence is stated in Corollary 3.14. The function g may be locally constant on intervals where the histogram of the upper image is zero (see top, middle graph). Indeed, on such intervals, the level sets are invariant.

68-1 Theorem 5. Suppose that u and v are two real-valued functions de�ned on ℝ^{N} such that for every λ ∈ \overline{ℝ} the following alternative holds: Either there exists a μ ∈ \overline{ℝ} such that

X_{λ}v = {x | u(x) ≥ μ} = X_{μ}u,  (I.46)

or there exists a μ ∈ \overline{ℝ} such that

X_{λ}v = {x | u(x) > μ} = ⋃_{ν>μ} X_{ν}u.  (I.47)

Then the contrast change g : \overline{ℝ} → \overline{ℝ} de�ned by

g(μ) = sup{λ | X_{μ}u ⊂ X_{λ}v}  (I.48)

maps u onto v, that is, v = g(u).

Corollary 2. If u and v are two real-valued functions such that all the level sets of v are level sets of u, then the contrast change g : \overline{ℝ} → \overline{ℝ} de�ned by

g(μ) = sup{λ | X_{μ}u ⊂ X_{λ}v}  (I.49)

maps u onto v.

Figure I.47: The original image (top, left) has a strictly positive histogram (all gray levels between 0 and 255 are represented). Therefore, if any contrast change g that is not strictly increasing is applied, then some data will be lost. Every level set of the transformed image g(u) is a level set of the original image; however, the original image has more level sets than the transformed image.

70-1 Linear smoothing can introduce new gray levels in an image. One of the main properties of contrast-invariant operators is that they do not create new grey levels in an image. This means in practice that the contrast in images �ltered by contrast-invariant operators is not degraded.

I.12 Properties of Contrast-invariant function operators

De�nition 6. We say that a function operator T is monotone if

u ≥ v ⟹ T u ≥ T v.  (I.50)

We say that a set operator T is monotone if

X ⊂ Y ⟹ T X ⊂ T Y.  (I.51)

De�nition 7. A function operator T de�ned on F is said to be contrast invariant with respect to a set G of contrast changes if g ∈ G and u ∈ F imply that g(u) ∈ F and that

g(T u) = T g(u).  (I.52)

Proposition 4. Let T be a monotone, Lipschitz contrast invariant operator. Then for every u ∈ F, R(T u) ⊂ \overline{Ru}. In particular, if Ru is �nite, then the range of T u is a �nite subset of the range Ru of u.

Proof: Given u ∈ F, we construct a continuous contrast change g such that g(s) = s for s ∈ \overline{Ru} and g(s) > s otherwise. One such function is given by

g(s) = s + \frac{1}{2}d(s, \overline{Ru}),

where d(s, X) denotes the distance from s to X, that is, d(s, X) = inf_{x∈X} |s − x|. Function g is nondecreasing and \frac{3}{2}-Lipschitz. Furthermore, g(s) = s if and only if s ∈ \overline{Ru}. In particular, g(u) = u. This and the contrast invariance of T imply that

T u = T g(u) = g(T u).

Hence, for every x ∈ ℝ^{N}, T u(x) is a �xed point of g. This means that T u(x) ∈ \overline{Ru}.

In practice, image operators are also translation invariant. For x ∈ ℝ^{N} we are going to use the notation τ_{x} to denote the translation operator for both sets and functions: For X ∈ L, τ_{x}X = {x + y | y ∈ X}, and for u ∈ F, τ_{x}u is de�ned by τ_{x}u(y) = u(y − x).

De�nition 8. A set operator T is said to be translation invariant if X ∈ L and x ∈ ℝ^{N} imply that τ_{x}X ∈ L and that

τ_{x}T X = T τ_{x}X.

A function operator T is said to be translation invariant if u ∈ F and x ∈ ℝ^{N} imply that τ_{x}u ∈ F and that

τ_{x}T u = T τ_{x}u.

We say that a function operator T commutes with the addition of constants if u ∈ F and c ∈ ℝ imply that u + c ∈ F and that T (u + c) = T u + c. Contrast-invariant operators clearly commute with the addition of constants: Consider the contrast change de�ned by g(s) = s + c.

Lemma 3. Let T be a translation-invariant monotone function operator de�ned on F that commutes with the addition of constants. If u ∈ F satis�es the Lipschitz condition |u(x) − u(y)| ≤ K|x − y|, then |T u(x) − Tu(y)| ≤ K|x − y|.

Proof. For any x, y, and z, we can write

u(y + z) − K|x − y| ≤ u(x + z) ≤ u(y + z) + K|x − y|. (I.53)

Consider x and y �xed and the terms in (I.53) to be functions of z. Since T is monotone, we can apply T to the functions in (I.53) and preserve the inequalities. If we do this and evaluate the result at z = 0, we have T (u(y+·)−K|x−y|)(0) ≤ Tu(x+·)(0) ≤ T(u(y+·)+K|x−y|)(0). Now use the fact that T commutes with the addition of constants to obtain T (u(y + ·)(0) − K|x − y| ≤ Tu(x + ·)(0) ≤ T(u(y + ·)(0) + K|x − y|. Next, use the translation invariance of T to write

T u(y) − K|x − y| ≤ Tu(x) ≤ Tu(y) + K|x − y|,

which is the announced result.

Exercise 6. Assume that the function u : ℝ^{N} → ℝ is uniformly continuous and write

|u(x + y) − u(x)| ≤ ε(|y|),

where x, y ∈ ℝ^{N} are arbitrary and ε is the modulus of continuity of u. Assume that T is a translation-invariant monotone operator that commutes with the addition of constants. Prove that T u is uniformly continuous and that the modulus of continuity of T u is dominated by ε. (Hint: Adapt the preceding proof.)

Exercise 7. Let T be a contrast-invariant monotone operator. Prove the following statements:

1. T u = c for every constant function u : ℝ^{N} → c.

2. u ≥ c implies T u ≥ c, and u ≤ c implies T u ≤ c.

3. sup |T u(x) − Tv(x)| ≤ sup |u(x) − v(x)|, where the supremum is taken over x ∈ ℝ^{N}. (Hint: Write − sup |u(x) − v(x)| ≤ u(x) − v(x) ≤ sup |u(x) − v(x)|.)

4. Show that (i), (ii), and (iii) remain true under the weaker assumption that T is monotone and commutes with the addition of constants.

I.13 From contrast-invariant function operators to set operators: The threshold superposition principle

We begin with a function operator T de�ned on a class of functions F, and we assume that F contains the characteristic functions 1_{X_{λ}u} for all u ∈ F. Then we can associate with T a set operator T de�ned on L = {X | X = X_{λ}u, u ∈ F, λ ∈ ℝ} as follows:

T X = X_{1}T 1_{X}.  (I.54)

If T is a monotone function operator, then T is a monotone set operator. Indeed, X ⊂ Y if and only if 1_{X} ≤ 1_{Y}, in which case T 1_{X} ≤ T 1_{Y} and X_{1}T 1_{X} ⊂ X_{1}T 1_{Y}.

Lemma 4. Let T be a monotone function operator de�ned on F and assume that T is Lipschitz contrast invariant. Consider the threshold function γ_{λ} de�ned by γ_{λ}(s) = 1 if s ≥ λ and γ_{λ}(s) = 0 otherwise. Then, for u ∈ F,

T γ_{λ}(u) = γ_{λ}(T u)

almost everywhere, except possibly for at most a countable number of λ’s.

Figure I.48: Approximation of the threshold function γ_{λ} from above and from below

80-1 Proof. γ_{λ}^{+ε} and γ_{λ}^{−ε} are Lipschitz-continuous contrast changes, and γ_{λ}^{−ε} ≤ γ_{λ} ≤ γ_{λ}^{+ε}. Since T is monotone and since T is invariant under Lipschitz-continuous contrast changes, for any u ∈ F:

γ_{λ}^{−ε}(T u) = T γ_{λ}^{−ε}(u) ≤ T γ_{λ}(u) ≤ T γ_{λ}^{+ε}(u) = γ_{λ}^{+ε}(T u).

As ε → 0, γ_{λ}^{−ε}(s) ↑ γ_{λ}^{−}(s) and γ_{λ}^{+ε}(s) ↓ γ_{λ}(s) for all s ∈ ℝ. Thus, in the limit we have

γ_{λ}^{−}(T u) ≤ T γ_{λ}(u) ≤ γ_{λ}(T u),

where γ_{λ}^{−}(s) = 1 if s > λ and γ_{λ}^{−}(s) = 0 otherwise. The crux of the proof is to show that

γ_{λ}^{−}(T u(x)) = γ_{λ}(Tu(x))

almost everywhere, except possibly for a countable number of λ’s. The only x for which γ_{λ}^{−}(T u(x)) ≠ γ_{λ}(Tu(x)) are the points where T u(x) = λ. Let Λ = Λ(T u) be the (countable) set of λ for which meas({x | Tu(x) = λ}) > 0. If λ ∉ Λ, then γ_{λ}^{−}(T u(x)) = γ_{λ}(Tu(x)).

De�nition 9. We say that the function operator T is obtained from the set operator T by the threshold superposition principle almost everywhere if

T u(x) = sup{λ | x ∈ T X_{λ}u}  (I.55)

for almost every x ∈ ℝ^{N}.

The relation (I.57) has practical implications. It means that T u can be computed by applying T separately to each characteristic function of the level sets X_{λ}u. This leads to the following stack �lter algorithm:

u X_{λ}u → T X_{λ}u ; ↗  ↘ ; ⋮ ; ↘  ↗ ; X_{μ}u → T X_{μ}u T u(x) = sup{λ | x ∈ T X_{λ}u}.

We assume that F contains 1_{X_{λ}u} for u ∈ F, L is the set of all level sets X_{λ}u, u ∈ F.

Proposition 5. Let T be a monotone contrast invariant function operator de�ned on F. De�ne the associated set operator T on L by T X = X_{1}T 1_{X}. Then

1. If u ∈ F, then

T X_{λ}u = X_{λ}T u  (I.56)

almost everywhere, except possibly for at most a countable number of λ’s.

2. The threshold superposition principle holds almost everywhere:

T u(x) = sup{λ | x ∈ T X_{λ}u} for a.e. x ∈ ℝ^{N}.  (I.57)

3. T is monotone, T ∅ = ∅, and T ℝ^{N} = ℝ^{N}.

I.14 From set operators to contrast invariant function operators.

Can we extend a set monotone operator into a contrast invariant and monotone operator de�ned on an adequate set of functions ? We start from a set operator T , de�ned on a subset L of P(ℝ^{N}) and with values in L. We assume that T is a monotone operator, that is, X ⊂ Y ⇒ T (X) ⊂ T (Y). We de�ne F as the set of all real functions de�ned on ℝ^{N} whose level sets X_{λ}u belong to L. As a consequence, F contains the characteristic functions of elements of L.

A natural de�nition for the function operator T associated with T is given by the threshold superposition principle,

T u(x) = sup{λ, x ∈ T (X_{λ}u)}

We have now to check whether T is contrast invariant.

Proposition 6. Let T , T → T be a monotone operator satisfying

T (∅) = ∅, T (ℝ^{N}) = ℝ^{N}

Then the operator T de�ned on F by superposition principle,

T u(x) = sup{λ, x ∈ T (X_{λ}u)}

satis�es, for almost every λ ∈ ℝ,

X_{λ}(T u) = T (X_{λ}(u)) almost everywhere in x

and for any nondecreasing continuous contrast change g,

g(T u) = T (g(u)).

Proof The relation

T u(x) = sup{λ, x ∈ T (X_{λ}u)}

implies that for almost every λ

X_{λ}(T u) = T (X_{λ}u) almost everywhere

that is for all λ in a subset Λ of ℝ such that meas(ℝ \ Λ) = 0. We notice that u ≤ v if and only if X_{λ}u ⊂ X_{λ}v for all λ in a dense countable subset of ℝ. We deduce immediately that T is monotone, ie u ≥ v ⇒ T u ≥ T v.

Let us now show that T commutes with contrast changes. Assume �rst that g is strictly increasing. By the superposition principle:

T (g(u(x))) = sup{λ, g(−∞) ≤ λ ≤ g(+∞), x ∈ T (X_{λ}g(u))}

= sup{g(μ), x ∈ T (X_{g(μ)}g(u))}
= sup{g(μ), x ∈ T (X_{μ}u)} = g(Tu(x)).

Then, approximate uniformly nondecreasing g’s by increasing ones and use the preceding result and the commutation of T with constants.

We shall now examine stronger conditions on T , which ensure commutation of T with all thresholds X_{λ} (and not only with almost all of them).

De�nition 10. We say that a set operator T is upper semicontinuous if for every nonincreasing sequence of sets X_{n} ∈ T , we have

T (⋂_{n} X_{n}) = ⋂_{n} T (X_{n}).

Theorem 6. Let T : T → T be a monotone upper semicontinuous set operator satisfying T (∅) = ∅ and T (ℝ^{N}) = ℝ^{N}. Then the function operator de�ned by

T u(x) = sup{λ, x ∈ T (X_{λ}u)}.

is contrast invariant and satis�es, for every λ,

X_{λ}(T u) = T (X_{λ}(u)).

The upper semicontinuity is necessary to ensure that a monotone set operator de�nes a function operator such that the commutation with thresholds X_{λ}(T u) = T (X_{λ}(u)) holds for every λ. Let us choose for example the following set operator T ,

T (X) = X if meas(X) > a and T (X) = ∅ otherwise .

Let u be the function from ℝ into ℝ de�ned by u(x) = −|x|. Then

T (u)(x) = sup{λ, x ∈ T (X_{λ}u)} = min(−|x|, −a/2)

Therefore

X_{−a/2}T (u) = [−a/2, a/2].

Now, X_{−a/2}u = [−a/2, a/2] ; its measure is a. Thus

T (X_{−a/2}u) = ∅ ≠ X_{−a/2}T (u),

which means that T is not contrast invariant.

I.15 The extrema killer: An application

T_{a}X = ⋃_{meas(X_{i})≥a} X_{i}.  (I.58)

We are going to prove that the small component killer is upper semicontinuous on descending sequences of compact subsets.

Lemma 5. Let (Y_{n})_{n∈ℕ} be a descending sequence of nonempty compact subsets of ℝ^{N}, and consider the intersection Y = ⋂_{n∈ℕ} Y_{n}. If the Y_{n} are connected, then Y is connected.

Figure I.49: Extrema killer: maxima killer followed by minima killer. The extrema killer removes all connected components of upper and lower level sets with area less than some threshold, which here equals 20 pixels.

92-1 Lemma 6. The small component killer T_{a} is upper semicontinuous on L, the set of all compact subsets of ℝ^{N}.

De�nition 11. The function operator T_{a} de�ned by

T_{a}u(x) = sup{λ | x ∈ T_{a}X_{λ}u}

on the set F of functions in L_{C} whose level sets are compact (as subsets of 𝕋^{N}) is called a maxima killer.

Theorem 7. The maxima killer T_{a} is monotone, contrast invariant, and translation invariant, and

X_{λ}T_{a}u = T_{a}X_{λ}u  (I.59)

for all λ. As a consequence of (I.59), no connected component of a level set of T u has measure less than a. Furthermore, T_{a} maps the continuous C-periodic functions into continuous C-periodic functions.

A maxima killer T_{a} cuts off the peaks of continuous functions, but it does nothing for the “valleys.” We could de�ne a minima killer to �ll the valleys by developing a dual theory based on lower level sets X_{λ}^{−}u = {x | u(x) ≤ λ}. It is, however, simpler to de�ne a minima killer T_{a}^{−} in terms of T_{a} as

T_{a}^{−}u = −T_{a}(−u).

It is clear that both of these operators have the same continuity properties. A good denoising process is to alternate T_{a} and T_{a}^{−}, as illustrated in Figures 6.2 and 6.3 . We note, however, that T_{a} and T_{a}^{−} do not necessarily commute, as is shown in the next exercise.

Exercise 8. Let N = 1 and take u(x) = sin x. Compute T_{a}u and T_{a}^{−}u and show that they commute on u if a ≤ π and do not commute if a > π. Following the same idea, construct function u de�ned on [−π, π]^{2} such that T_{a}T_{a}^{−}u ≠ T_{a}^{−}T_{a}u.

Figure I.50: Extrema killer: maxima killer followed by minima killer. Above, left: original image. Above, right: image after extrema killer removes connected components of 20 pixels or less. Below: level lines (levels of multiples of 4) of the image before and after the application of the extrema killer.

94-1

I.16 Erosions and Dilations

T_{d}u(x) = sup_{y∈x+B} u(y) and T_{e}u(x) = inf_{y∈x+B} u(y).

The operator T_{d} is called a dilation, and T_{e} is called an erosion. These function operators have associated set operators, and in fact we are going to de�ne the set operators �rst.

De�nition 12. Let B be a subset of ℝ^{N} and let t ≥ 0 be a scale parameter. The set operators D_{t} and E_{t} are de�ned on subsets X ∈ P(ℝ^{N}) by

D_{t}X = X + tB = {z | z = x + y, x ∈ X, y ∈ tB} = ⋃_{y∈tB}(X + y),

(I.60)

and

E_{t}X = {x | x + tB ⊂ X} = ⋂_{y∈−tB}(X + y).  (I.61)

D_{t}X is called the dilation of X with structuring element B at scale t. E_{t}X is called the erosion of X with structuring element B at scale t.

Figure I.51: Dilation of a set. Left to right: A set; its dilation by a ball of radius 20; the difference set.

When B is symmetric with respect to zero, the operator D_{t}E_{t} is called an opening at scale t and E_{t}D_{t} is called a closing at scale t. These names have a topological origin. If B is the open ball centered at zero with radius one, then the opening at scale t of a set X is the union of all balls with radius t contained in X. The interior of X is the union of all open balls contained in X; it is also the largest open set contained in X. If we call the interior map T^{∘}X = X^{∘} the opening, then an opening at scale t appears as a quali�ed opening. The topological statement “the closure of the complement of X is the complement of the interior of X” has its counterpart for openings and closings at scale t.

Figure I.52: Opening of a set as curvature threshold from above. Left to right: A set X; its opening by a ball of radius 20; the difference set. This opening transforms X into the union of all balls of radius 20 contained in it. The resulting operation can be understood as a threshold from above of the curvature of the set boundary.

1. Let B = {x | |x| < 1}. Then D_{t}E_{t}X is the union of all open balls with radius t contained in X.

2. Let B be any structuring element that is symmetric with respect to zero. Writing X^{c} = ℝ^{N} \ X, one has D_{t}X^{c} = (E_{t}X)^{c}. This implies E_{t}D_{t}X^{c} = (D_{t}E_{t}X)^{c}.

Figure I.53: Erosion of a set. Left to right: A set; its erosion by a ball of radius 20; the difference set.

100-1 De�nition 13. For u ∈ F, de�ne D_{t}u by

D_{t}u(x) = sup_{y∈tB} u(x − y)  (I.62)

and de�ne E_{t}u by

E_{t}u(x) = inf_{y∈−tB} u(x − y).  (I.63)

Theorem 8. For u ∈ F,

D_{t}u(x) = sup{λ | x ∈ D_{t}X_{λ}u} for a.e. x ∈ ℝ^{N},  (I.64)
E_{t}u(x) = sup{λ | x ∈ E_{t}X_{λ}u} for a.e. x ∈ ℝ^{N},  (I.65)
D_{t}X_{λ}u = X_{λ}D_{t}u E_{t}X_{λ}u = X_{λ}E_{t}u  (I.66)

almost everywhere, except possibly for at most a countable number of λ’s.

Figure I.54: Closing of a set as a curvature threshold from below. Left to right: A set X; its closing by a ball of radius 20; the difference set. The closing of X is just the opening of X^{c}. It can be viewed as a threshold from below of the curvature of the set boundary.

102-1 We say that the family of dilations {D_{t} | t > 0} associated with a structuring element B is recursive if D_{t}D_{s} = D_{t+s} for all s, t > 0, and similarly for {E_{t} | t > 0}. (A recursive family is also called a semigroup.)

Lemma 7. B is convex if and only if (s + t)B = sB + tB for all s, t ≥ 0.

First row: a sea bird image and its level lines for all levels multiple of 12. Second row: an erosion with radius 4. On the right, the resulting level lines where the circular shape of the structuring element (a disk with radius 4) appears around each local minimum of the original image. Erosion removes local maxima (in particular, all small white spots) but expands minima. Thus, all dark spots, like the eye of the bird, are expanded. Third row: the effect of a dilation with radius 4 and the resulting level lines. We see how local minima are removed (for example, the eye of the bird) and how white spots on the tail expand. Here, in turn, circular level lines appear around all local maxima of the original image.

103-1 Proposition 7. The dilations D_{t} and D_{t} and the erosions E_{t} and E_{t} are recursive if and only the structuring element B is convex.

Figure I.55: Openings and closings of a natural image. First row: the original image and its level lines for all levels multiple of 12. Second row: an opening with radius 4 has been applied. Third row: a closing with radius has been applied. We can recognize the circular shape of the structuring element in the level lines displayed on the right.

104-1

I.17 The PDEs associated with erosions and dilations

Scaled dilations and erosions are associated with the equations ∂u/∂t = ±|Du|. To explain this connection, we begin with a bounded convex set B that contains the origin, and we de�ne the gauge ∥ · ∥_{B} on ℝ^{N} associated with B by ∥x∥_{B} = sup_{y∈B}(x · y). If B is a ball centered at the origin with radius one, then ∥ · ∥_{B} is the usual Euclidean norm, which we write simply as | · |.

Proposition 8 (Hopf–Lax formula). Assume that B is a bounded

convex set in ℝ^{N} that contains the origin. Given u_{0} : ℝ^{N} → ℝ, de�ne u : ℝ^{+} × ℝ^{N} → ℝ by u(t, x) = D_{t}u_{0}(x). Then u satis�es the equation

\frac{∂u}{∂t} = ∥Du∥_{−B}

at each point (t, x) where u has continuous derivatives in t and x. The same result hold when D_{t} is replaced with E_{t} and the equation is replaced with ∂u/∂t = −∥Du∥_{−B}.

Proof.

u_{0}(x − y) − u_{0}(x) = −Du_{0}(x) · y + o(|y|),

u(h, x) − u(0, x) = sup_{y∈hB}(−Du_{0}(x) · y + o(|y|)).
u(h, x) − u(0, x) = h sup_{z∈B}((−Du_{0}(x) · z) + o(|h|).

We can divide both sides by h and pass to the limit as |h| → 0 to obtain

\frac{∂u}{∂t}(0, x) = ∥Du_{0}(x)∥_{−B},

which is the result for t = 0. For an arbitrary t > 0, we have D_{t+h} = D_{t}D_{h} = D_{h}D_{t}, and we can write

u(t + h, x) − u(t, x) = D_{h}u(t, ·)(x) − u(t, x).

By repeating the argument made for t = 0 with u_{0} replaced with u(t, ·), we arrive at the general result.

Figure I.56: Denoising based on openings and closings. First row: scanned picture of the word “operator” with black dots and a black line added; a dilation with a 2 × 5 rectangle; an erosion with the same structuring element applied to the middle image. The resulting operator is a closing. Small black structures are removed by such a process. Second row: the word “operator” with a white line and white dots inside the letters; erosion with a rectangle 2 × 5; a dilation with the same structuring element applied to the middle image. The resulting operator is an opening. This time, small white structures are removed.

107-1 Figure I.57: Examples of denoising based on opening or closing, as in Figure I.56. Perturbations made with both black and white lines or dots have been added to the “operator” image. First column, top to bottom: original perturbed image; erosion with a 1 × 3 rectangle; then dilation with the same structuring element. (In other words, opening with this rectangle.) Then a dilation is applied with a rectangle 3×1, and �nally an erosion with the same rectangle. Second column: The same process is applied, but with erosions and dilations exchanging their roles. It does not work so well because closing expands white perturbations and opening expands black perturbations. These operators do not commute. See Figure ??, where an application of the median �lter is more successful.

107-2

I.18 From set operators to contrast invariant function operators

De�nition 14. We say that a set operator T is upper semicontinuous if for every strictly decreasing sequence of compact sets X_{n}, we have

T (⋂_{n} X_{n}) = ⋂_{n} T (X_{n}).

(Strictly decreasing : closure of X_{n+1} contained in X_{n}^{∘}).

Theorem 9. Let T : be a monotone upper semicontinuous set operator on compact periodic sets. Then the function operator de�ned on continuous periodic functions by superposition principle,

T u(x) = sup{λ, x ∈ T (X_{λ}u)}

is contrast invariant and satis�es, for every λ,

X_{λ}(T u) = T (X_{λ}(u)).

I.19 Conversely : the Evans–Spruck extension

Theorem 10 (Evans–Spruck). Let T be a contrast-invariant monotone operator de�ned on continuous C-periodic functions. Then there is a unique set operator T de�ned on C-periodic compact sets satisfying the commutation relation

T X_{λ}u = X_{λ}T u.  (I.67)

In addition, T is upper semicontinuous on compact sets.

Lemma 8. Let u and v be continuous periodic functions bounded above by zero and assume that X_{0}u = X_{0}v (≠ ∅). Then there is a continuous contrast change h such that h(0) = 0 and u ≥ h(v).

110-1 Proof. (Evans-Spruck) It is suf�cient to prove that X_{0}u = X_{0}v implies X_{0}T u = X_{0}T v when u and v are bounded from above by zero.

The proof follows directly from the preceding lemma : Let h be a continuous contrast change such that h(0) = 0 and u ≥ h(v). Since T is monotone and contrast invariant, T u ≤ 0 and T u ≥ T h(v) = h(T v). Using the fact that h(0) = 0, we see that T v(x) = 0 implies that T u(x) = 0. By interchanging the roles of u and v we see that T u(x) = 0 implies that T v(x) = 0 and conclude that X_{0}T u = X_{0}T v.

Let us prove that T is upper semicontinuous on compact sets.

Lemma 9. (Urysohn). Let X_{n} be a strictly decreasing sequence of periodic compact sets, namely X_{n+1} ⊂ X_{n}^{∘}. Then there exists a continuous periodic function u such that X_{1− \frac{1}{n}} u = X_{n} and X_{1}u = ∩_{n}X_{n}.

By assumption,

X_{1}(T u) = T (X_{1}u) = T (∩_{n}X_{1− \frac{1}{n}} u) = T (∩_{n}X_{n}).

On the other hand, X_{1}(T u) = ∩_{n}X_{1− \frac{1}{n}} T u = ∩_{n}T (X_{1− \frac{1}{n}} u) = ∩_{n}T (X_{n}). Combining both results, we get

T (∩_{n}X_{n}) = ∩_{n}T (X_{n}).

Figure I.58: Evans-Spruck extension and the chessboard dilemma. Left: a chessboard image. Next: result with the �nite difference scheme for the curvature motion, applied up to a �xed scale. The creation of a new gray level proves that the scheme is not fully contrast invariant. The next two images are obtained by the Evans–Spruck extension of the curvature motion, �rst applied to the white set considered as closed, second to the black set considered as closed. The median �lters Med and Med^{−} will yield exactly the same dilemma.

111-1

I.20 Contrast-Invariant Monotone Operators as Sup-Inf Operators

Two representation theorems: If T is a monotone contrast invariant translation invariant function operator, then it has a “sup-inf” representation of the form

T u(x) = sup_{B∈B} inf_{y∈B} u(x + y),

where B is a family of subsets of P(ℝ^{N}), the set of all subsets of ℝ^{N}. This theorem is a nonlinear version of the Riesz theorem that states that a continuous linear translation-invariant operator can be represented as a convolution

T u(x) = ∫_{ℝ^{N}} u(x − y)k(y) dy.

In this case, the kernel k is called the impulse response. In the same way, B is an impulse response for the nonlinear operator.

I.21 Translation-invariant monotone set operators

Theorem 11 (Matheron). Let T be a translation-invariant monotone set operator acting on the set of periodic compact sets L. Then

T X = ⋃_{B∈B} ⋂_{y∈B}(X − y) = {x | x + B ⊂ X for some B ∈ B}, (I.68)

where B = {X | 0 ∈ T X}. Conversely, if B is any nonempty subset of P(ℝ^{N}), then formula (I.68) de�nes a translation-invariant monotone set operator on P(ℝ^{N}).

B is called set of structuring elements in the language of mathematical morphology.

I.22 Sup-inf operators

Theorem 12. Let T be a monotone function contrast invariant operator on F (continuous periodic functions). Let T be the set operator associated with T by the Evans-Spruck extension T X = X_{1}T 1_{X}. For x ∈ ℝ^{N} and X ∈ L, let B_{x} = {X | x ∈ T X}. Then for u ∈ F,

T u(x) = sup_{B∈B_{x}} inf_{y∈B} u(y) for a.e. x ∈ ℝ^{N}.  (I.69)

If T is also translation invariant, then

T u(x) = sup_{B∈x+B_{0}} inf_{y∈B} u(y) for a.e. x ∈ ℝ^{N}.  (I.70)

Conversely, equation (I.69) de�nes a contrast-invariant monotone function operator, and equation (I.70) de�nes a contrast-invariant monotone function operator that is also translation invariant.

Proof. For u ∈ F, let T̃u(x) = sup_{B∈B_{x}} inf_{y∈B} u(y), where B_{x} = {X | x ∈ T X}. We know that X_{μ}T u = T X_{μ}u. We are going to prove that T̃u(x) = Tu(x) by showing that

T̃u(x) ≥ λ ⟺ Tu(x) ≥ λ  (I.71)

for λ ∈ Λ and x ∈ ℝ^{N} \ N.

T u(x) ≥ λ ⟺^{(1)} T u(x) ≥ μ for all μ < λ ⟺^{(2)} x ∈ X_{μ}Tu for all μ < λ

⟺^{(3)} x ∈ T X_{μ}u for all μ < λ ⟺^{(4)} X_{μ}u ∈ B_{x} for all μ < λ

⟺^{(5)} There is a B ∈ B_{x} such that inf_{y∈B} u(y) ≥ μ for all μ < λ

⟺^{(6)} sup_{B∈B_{x}} inf_{y∈B} u(y) ≥ λ ⟺^{(7)} T̃u(x) ≥ λ.

If T is also translation invariant, then B_{x} = x + B_{0}.

To prove the converse, de�ne T on {u | u : ℝ^{N} → ℝ} by

T u(x) = sup_{B∈B_{x}} inf_{y∈B} u(y).

The operator T is well de�ned, since we have assumed that B_{x} is not empty. Clearly, T is monotone: If u ≥ v, then inf_{y∈B} u(y) ≥ inf_{y∈B} v(y), and it follows that T u ≥ T v. If g is any continuous contrast change, then for any nonempty set I ⊂ ℝ, inf_{x∈I} g(x) = g(inf_{x∈I} x) and sup_{x∈I} g(x) = g(sup_{x∈I} x). This implies that T commutes with g. If we replace B_{x} with x + B_{0}, then it is easy to check that T is also translation invariant.

116-1 Exercise 9. Use these last results to show that the extrema killer T_{a} can be represented as a sup-inf function operator with the structuring elements

B = {B | B is compact, connected, meas (B) = a, and0 ∈ B}.

Sup-inf operators can be represented as inf-sup operators.

Proposition 9. Assume that T can be represented by T u(x) = sup_{B∈B} inf_{y∈x+B} u(y). Then T also has an inf-sup representation

T u(x) = inf_{B∈B′} sup_{y∈x+B} u(y),  (I.72)

where B′ is a new structuring set.

Proof. A map ϕ : B → ⋃_{B∈B} B with the property that ϕ(B) ∈ B is called a selection function. Let S_{B} denote the set of all selection functions. We then have

sup_{B∈B} inf_{y∈B} u(y) = inf_{ϕ∈S_{B}} sup_{B∈B} u(ϕ(B)) = inf_{ϕ∈S_{B}} sup_{y∈ϕ(B)} u(y).

For ϕ ∈ S_{B} write ϕ(B) = ⋃_{B∈B} ϕ(B) and let D = {ϕ(B) | ϕ ∈ S_{B}}. Then

sup_{B∈B} inf_{y∈B} u(y) = inf_{ϕ(B)∈D} sup_{y∈ϕ(B)} u(y).

I.23 Median Filters and Mathematical Morphology

All of the median operators (or �lters) will be de�ned in terms of a nonnegative measurable weight function k : ℝ^{N} → [0, +∞) that is normalized:

∫_{ℝ^{N}} k(y) dy = 1.

The k-measure of a subset B ⊂ ℝ^{N} is denoted by |B|_{k} and de�ned by

|B|_{k} = ∫_{B} k(y) dy = ∫_{ℝ^{N}} k(y)1_{B}(y) dy.

Clearly, 0 ≤ |B|_{k} ≤ 1. Example : k = c_{N}^{−1}(r)1_{B(0,r)}, where B(0, r) denotes the ball of radius r centered at the origin and c_{N}(r) is the Lebesgue measure of B(0, r).

I.24 The median operators Med_{k} and Med_{k}

As we did for the erosions and dilations, we �rst de�ne the set operators.

De�nition 15. Let X be a measurable subset of ℝ^{N} and let k be a weight function. The median set of X weighted by k, which we denote by Med_{k}X, is de�ned by

Med_{k}X = \{x | |X − x|_{k} ≥ \frac{1}{2}\}.  (I.73)

When k is a Gaussian : Koenderink and van Dorn. If k = (1/πr^{2})1_{B(0,r)}. Then x belongs to Med_{k}X if and only if the Lebesgue measure of X ∩ B(x, r) is greater than or equal to half the measure of B(0, r). Intuitively, x ∈ Med_{k}X if points of X are in the majority around x.

Lemma 10. The set operator Med_{k} is upper semicontinuous on the measurable sets L.

Proof. Med_{k}X = \{x | |X − x|_{k} ≥ \frac{1}{2}\}. Given a descending sequence (X_{n})_{n∈ℕ} in L, we must show that

Med_{k} ⋂_{n∈ℕ} X_{n} = ⋂_{n∈ℕ} Med_{k}X_{n}.

Since Med_{k} is clearly monotone, it is always true that Med_{k} ⋂_{n∈ℕ} X_{n} ⊂ ⋂_{n∈ℕ} Med_{k}X_{n}. To prove the other inclusion, assume that x ∈ ⋂_{n∈ℕ} Med_{k}X_{n}. Then by the de�nition of Med_{k}, |X_{n} − x| ≥ 1/2 for all n ∈ ℕ. Since X_{n} − x ↓ ⋂_{n∈ℕ}(X_{n} − x), we know from the countable additivity of Lebesgue measure that |X_{n} − x|_{k} ↓ | ⋂_{n∈ℕ}(X_{n} − x)|_{k}. This means that | ⋂_{n∈ℕ}(X_{n} − x)|_{k} ≥ 1/2, and hence that x ∈ Med_{k} ⋂_{n∈ℕ}(X_{n} − x).

De�nition 16. The function operator Med_{k} is de�ned on the set of all measurable functions F by the superposition principle:

Med_{k}u(x) = sup{λ | x ∈ Med_{k}X_{λ}u}.

Proposition 10. The median operator Med_{k} has the sup-inf representation

Med_{k}u(x) = sup_{B∈B} inf_{y∈x+B} u(y),  (I.74)

where B = {B | B ∈ L, |B|_{k} = 1/2}. Furthermore,

X_{λ}Med_{k}u = Med_{k}X_{λ}u  (I.75)

for all λ ∈ \overline{ℝ}.

Proof. We mention again that Med_{k} is monotone. It is also clearly translation invariant. Thus, by the sup-inf representation theorem, Med_{k} satis�es

Med_{k}u(x) = sup_{B∈x+B_{0}} inf_{y∈B} u(y),

where B_{0} = {B | 0 ∈ Med_{k}B}.

Figure I.59: Example of denoising with a median �lter. Left to right: scanned picture of the word “operator” with perturbations and noise made with black or white lines and dots; the image after one application of a median �lter with a circular neighborhood of radius 2; the image after a second application of the same �lter. Compare with the denoising using openings and closings (Figure I.57).

Exercise 10. Let k = (1/2)1_{[−1,1]}. This implies that |B|_{k} = 1/2 if and only if meas ([−1, 1] ∩ B) = 1.

Let u(x) = x if x ≤ 0, and u(x) = 0 if x > 0. Check that Med_{k}u(0) = 0, so that the median value coincides with the supremum and is different from the mean value on the neighborhood [−1, 1].

I.25 The median operators Med_{k}^{−} and Med_{k}^{−}

The median operator Med_{k}, as de�ned, is not invariant under “reverse contrast,” that is, it does not satisfy −Med_{k}u = Med_{k}(−u) for all measurable u. We de�ne Med_{k}^{−} by

Med_{k}^{−}u(x) = inf_{|B|_{k}≥ \frac{1}{2}} sup_{y∈x+B} u(y).  (I.76)

Then we clearly have the relation

Med_{k}^{−}u = −Med_{k}(−u),  (I.77)

If we now de�ne Med_{k}^{−} in terms of Med_{k}^{−} by Med_{k}^{−}X = X_{1}Med_{k}^{−}1_{X}, then using (I.78) and the relation 1_{X} + 1_{X^{c}} = 1,

Med_{k}^{−}X = (Med_{k}X^{c})^{c}.  (I.78)

It is easily checked that

Med_{k}^{−}X = \{x | |X − x|_{k} > \frac{1}{2}\}.  (I.79)

We say that k is not separable if |B|_{k} ≥ 1/2 and |B′|_{k} ≥ 1/2 imply that \overline{B} ∩ \overline{B′} ≠ ∅.

Proposition 11.

1. For every function u, Med_{k}u ≥ Med_{k}^{−}u.

2. Assume that k is not separable. Then Med_{k}u = Med_{k}^{−}u.

Proof. Let λ = Med_{k}u(0) = sup_{|B|_{k}≥1/2} inf_{y∈B} u(y). Consider the level set X_{λ+ε}u. Then inf_{y∈X_{λ+ε}} u(y) ≥ λ + ε. Thus |X_{λ+ε}u|_{k} < 1/2, since inf_{y∈B} < λ + ε for any set B such that |B| ≥ 1/2. Then |(X_{λ+ε}u)^{c}|_{k} ≥ 1/2, and by the de�nition of level sets, sup_{y∈(X_{λ+ε}u)^{c}} u(y) ≤ λ + ε. This implies that

inf_{|B|_{k}≥ \frac{1}{2}} sup_{y∈B} u(y) ≤ λ + ε.

k is not separable implies that for all B and B′ having k-measure greater than or equal to 1/2, we have inf_{y∈\overline{B}} u(y) ≤ sup_{y∈\overline{B′}} u(y). If u is continuous, then we have inf_{y∈B} u(y) ≤ sup_{y∈B′} u(y). Since B and B′ were arbitrary except for the conditions |B|_{k} ≥ 1/2 and |B′|_{k} ≥ 1/2, the last inequality implies that

sup_{|B|_{k}≥ \frac{1}{2}} inf_{y∈B} u(y) ≤ inf_{|B′|_{k}≥ \frac{1}{2}} sup_{y∈B′} u(y).

This yields Med_{k}u = Med_{k}^{−}u when u is continuous.

Exercise 11. The non separability condition is necessary. Let k = \frac{1}{2}1_{[−2,−1]∪[1,2]} and let u(x) = −1 if x ≤ −1, u(x) = 1 if x ≥ 1, u(x) = x elsewhere. Check that Med_{k}u(0) ≠ Med^{−}u(0).

Exercise 12. Saying that k is not separable is a fairly weak assumption. It corresponds roughly to saying that the support of k cannot be split into two disjoint connected components each having k-measure 1/2. Show that if k is continuous and if its support is connected, then it is not separable.

Figure I.60: Denoising based on a median �lter. Left: an image altered on 40% of its pixels with salt and pepper noise. Right: the same image after three iterations of a median �lter with a 3 × 3 square mask.

I.26 Discrete median �lters and the "usual" median value

Medu(x) = sup_{{P ∈P(N) \atop card(P)≥N/2}} inf_{i∈P} u(x − x_{i}),

Med^{−}u(x) = inf_{{P ∈P(N) \atop card(P)≥N/2}} sup_{i∈P} u(x − x_{i}).

This shows that Med = Med^{−} if and only if N is odd. In general, a median �lter with an odd number of pixels is preferred, since Med = Med^{−} in this case.

133-1 The de�nition of the discrete median �lter Med corresponds to the usual statistical de�nition of the median of a set of data: If the given data consists of the numbers y_{1} ≤ y_{2} ≤ ⋯ ≤ y_{N} and N = 2n + 1, the median is y_{n+1}. In case N = 2n, the median is (y_{n} + y_{n+1})/2. In both cases, half of the terms are greater than or equal to the median and half of the terms are less than or equal to the median. The usual median minimizes the functional ∑_{i=1}^{N} |y_{i} − y|.

\underbrace{y_{1} = ⋯ = y_{1}}_{k_{1} terms} ≤ ⋯ ≤ \underbrace{y_{i} = ⋯ = y_{i}}_{k_{i} terms} ≤ ⋯ ≤ \underbrace{y_{N} = ⋯ = y_{N}}_{k_{N} terms}. (I.80)

If K is odd, Med_{k}^{−} = Med_{k}, Med_{k}u(x) is equal to the ordinary median of the ordered set.

Exercise 13. 1. Show that

\frac{1}{N} ∑_{i=1}^{N} x_{i} = arginf_{m} ∑_{i=1}^{N}(x_{i} − m)^{2}.

2. Show that

Med^{−}(x_{i}) ≤ arginf_{m} ∑_{i=1}^{N} |x_{i} − m| ≤ Med(x_{i}).

Exercise 14. 1. Let k = 1_{B}, where B is set with Lebesgue measurable equal to one. Let Med_{B}u denote the median value of u in B, de�ned by Med_{B}u = Med_{k}u(0). Consider a bounded measurable function u de�ned on B. Show that

∫_{B} u(x) dx = arginf_{m} ∫_{B}(u(x) − m)^{2} dx

and that

Med_{B}^{−}u ≤ arginf_{m} ∫_{B} |u(x)−m| dx = \frac{Med_{B}^{−}u + Med_{B}u}{2} ≤ Med_{B}u.

2. Deduce from the above that the mean value is the best constant approximation in the L^{2} norm and that the median is the best constant approximation in the L^{1} norm.

136

Figure I.61: Smoothing effect of a median �lter on level lines. Above, left to right: original image; all of its level lines (boundaries of level sets) with levels multiple of 12; level lines at level 100. Below,left to right: result of two iterations of a median �lter with a disk with radius 2; corresponding level lines (levels multiple of 12); level lines at level 100.

I.27 Asymptotic Behavior of Contrast-Invariant Isotropic Operators, Dimension Two

We are going to study the monotone contrast invariant operators under the inf-sup form de�ned by

T u(x) = inf_{B∈B} sup_{y∈x+B} u(y),  (I.81)

where B is a given set of structuring elements.

We assume that the set of structuring elements B is bounded and isotropic. De�ne the scaled operators T_{h} by

T_{h}u(x) = inf_{B∈B} sup_{y∈x+hB} u(y).

We will prove that in the limit, as h tends to zero, the action of T_{h} on smooth functions is not as varied as one might expect given the possible sets of structuring elements. As an example, we will show that if T_{h} is a scaled median operator, then

T_{h}u(x) − u(x) = h^{2}C|Du(x)|curv(u)(x) + o(h^{2}),

where the constant C depends only on the function k used to de�ne the median operator. Thus, the operator |Du|curv(u) plays the same role for a class of contrast-invariant operators as the Laplacian Δu does for linear operators.

I.28 Asymptotic behavior theorem in ℝ^{2}

Let T be a monotone contrast invariant (and therefore an inf-sup operator)) generated by a set of structuring elements B that is bounded and isotropic, and de�ne the function H(s), s ∈ ℝ, by

H(s) = T [x + sy^{2}](0),  (I.82)

where T [x + sy^{2}] denotes “T u with u(x, y) = x + sy^{2}.” We also note the following relations for 0 < h ≤ 1:

T_{h}[x + sy^{2}](0) = hT [x + hsy^{2}](0) = hH(hs),  (I.83)
T_{h}[x](0) = hT [x](0) = hH(0)  (I.84)

Theorem 13. Let B be a set of structuring elements that is invariant under the group of linear isometries of ℝ^{2} and assume that B is bounded by M. De�ne the operator T by T u(x) = inf_{B∈B} sup_{y∈x+B} u(y), and de�ne the scaled operators T_{h}, 0 < h ≤ 1, by T_{h}u(x) = inf_{B∈B} sup_{y∈x+hB} u(y). Then, for any C^{2} function u : ℝ^{2} → ℝ,

T_{h}u(x) − u(x) = hT[x](0)|Du(x)| + O(h^{2}).

Proof. The invariance of B under the action of linear isometries implies that all of the T_{h} are invariant under rotations. Now suppose that p = |Du(x)| > 0. Then by a suitable rotation, we may assume that Du(x) = (|Du(x)|, 0), and the �rst-order expansion of u in a neighborhood of x can be written as

u(x + y) = u(x) + px + O(x, |y|^{2}),  (I.85)

where y = (x, y) ∈ D(0, M). Hence,

u(x+hy)−u(x) ≤ phx+Ch^{2}|y|^{2} and phx ≤ u(x+hy)−u(x)+Ch^{2}|y|^{2}

(I.86)

for all y ∈ D(0, M).

I.28.1 The asymptotic behavior of T_{h} when T [x](0) = 0

If T [x](0) = 0, then Theorem 10.1 is true but not very interesting. On the other hand, operators for which T [x](0) = 0 are interesting. If we consider T u(x) to be a kind of average of the values of u in a neighborhood of x, then assuming that T [x](0) = 0 makes sense and leads to interesting results. This means, however, that we must consider the next term in the expansion of T_{h}u; to do so we need to assume that u is C^{3}.

Theorem 14. Let B be a set of structuring elements that is isometric and bounded, and let M denote the bound. De�ne the operator T by T u(x) = inf_{B∈B} sup_{y∈x+B} u(y), and de�ne the scaled operator T_{h}, 0 < h ≤ 1, by T_{h}u(x) = inf_{B∈B} sup_{y∈x+hB} u(y). Also assume that T [x](0) = 0. Then for any C^{3} function u : ℝ^{2} → ℝ :

T_{h}u(x) − u(x) = h|Du(x)|H(\frac{1}{2}h curv(u)(x)) + O(h^{3})

Figure I.62: When s is small, the parabola P (s) with equation x+sy^{2} = m divides D(0, 1) into two components. The median value m(s) of x + sy^{2} on D(0, 1) simply is the value m for which these two components have equal area.

I.29 Median �lters and curvature motion in ℝ^{2}

Everything is now in place to investigate the asymptotic behavior of the scaled median �lter Med_{k_{h}}, which is represented by

Med_{k_{h}}u(x) = inf_{B∈hB} sup_{y∈x+B} u(y),

where hB = {B | |B|_{k_{h}} = 1/2, B ⊂ D(0, h)}.

145-1

Lemma 11.

Med_{k}[x + sy^{2}](0) = \frac{s}{3} + O(|s|^{3}).

Theorem 15. If u : ℝ^{2} → ℝ is C^{2}, then we have the following expansions:

1. On every compact set K ⊂ {x | Du(x) ≠ 0},

Med_{k_{h}}u(x) = u(x) + \frac{1}{6}|Du(x)|curv(u)(x)h^{2} + O(x, h^{3}),

where |O(x, h^{3})| ≤ C_{K}h^{3} for some constant C_{K} that depends only on u and K.

2. On every compact set K in ℝ^{2},

|Med_{k_{h}}u(x) − u(x)| ≤ \frac{1}{2}∥D^{2}u(x)∥h^{2} + C_{K}h^{3}

where the constant C_{K} depends only on u and K.

Our second example is called the Catté–Dibos scheme. Theorem 16. Let B be the set of all line segments of length 2 centered at the origin of ℝ^{2}. De�ne the operators SI_{h} and IS_{h} by SI_{h}u(x) = sup_{B∈hB} inf_{y∈x+B} u(y) and IS_{h}u(x) = inf_{B∈hB} sup_{y∈x+B} u(y). If u : ℝ^{2} → ℝ is C^{2} and |Du(x)| ≠ 0, then

\frac{1}{2}(IS_{h} + SI_{h})u(x) = u(x) + h^{2} \frac{1}{4}curv(u)(x)|Du(x)| + O(h^{3}).

Figure I.63: The result of smoothing with an erosion is independent of the curvature of the level lines. Left: image of a simple shape. Right: difference of this image and its eroded image. Note that the width of the difference is constant.

Figure I.64: Median �lter and the curvature of level lines. Smoothing with a median �lter is related to the curvature of the level lines. Left: image of a simple shape. Right: difference of this image with itself after it has been smoothed by one iteration of the median �lter. We see, in black, the points which have changed. The width of the difference is proportional to the curvature, as indicated by Theorem 10.4.

148-1 Figure I.65: Fixed point property of the median �lter, showing its griddependence. Left: original image. Right: result of 46 iterations of the median �lter with a radius of 2. The resulting image turns out to be a �xed point of this median �lter. This is not in agreement with Theorem 10.4, which shows that median �lters move images by their curvature: The image on the right clearly has nonzero curvatures! However, the discrete median �lter that we have applied here operating on a discrete image is grid-dependent and blind to small curvatures.

148-2 Figure I.66: Comparing an iterated median �lter and a median �lter. Top-left: original image. Top-middle: 16 iterations of the median �lter with a radius 2, Top-right: one iteration of the same median �lter with a radius 8. Below each image are the level-lines for grey levels equal to multiples of 16. This shows that iterating a small size median �lter provides more accuracy and less shape mixing than applying a large size median �lter. Compare this with the Koenderink– Van Doorn shape smoothing and the Bence–Merriman–Osher iterated �lter in Chapter 2, in particular Figures 2.9, 2.10, and 2.12.

148-3 Figure I.67: Consistency of the median �lter and of the Catté–Dibos. Top row: the sea bird image and its level lines for all levels equal to multiples of 12. Second row: a median �lter on a disk with radius 2 has been iterated twice. Third row: an inf-sup and then a sup-inf �lter based on segments have been applied. On the right: the corresponding level lines of the results, which, according to the theoretical results (Theorems 10.4 and 10.5), must have moved at a speed proportional to their curvature. The results are very close. This yields a cross validation of two very different numerical schemes that implement curvaturemotion based smoothing.

148-4