F Baus

M Nikolova

G Steidl

Fully smoothed ℓ 1 -T V models: Bounds for the minimizers and parameter choice

Keywords: Parameter estimation for smoothed ℓ 1 -TV model, Histogram specification, Quantisation noise, ℓ ∞ error, Convex optimization

. Here we focus on the way the parameters involved in these functionals determine the features of the minimizers û. We give explicit relationships between the minimizers and these parameters.

Given an input digital image f , we prove that the error û-f ∞ obeys b-ε ≤ û-f ∞ ≤ b where b is a constant independent of the input image. Further we can set the parameters so that ε > 0 is arbitrarily close to zero. More precisely, we exhibit explicit formulae relating the model parameters, the input image f and the values b and ε. Conversely, we can fix the parameter values so that the error ûf ∞ meets some prescribed b, ε. All theoretical results are confirmed using numerical tests on natural digital images of different sizes with disparate content and quality.

Introduction

In [START_REF] Nikolova | Exact histogram specifcation for digital images using a variational approach[END_REF] a variational method using C 2 smoothed ℓ 1 -TV functionals were proposed. The goal was to process digital (quantized) images so that the obtained minimizer is quite close to the input digital image but its pixels are real-valued and can be ordered in a strict way. Indeed, the obtained minimizers were shown to enable faithful exact histogram specification outperforming the state-of-the-art methods [START_REF] Coltuc | Exact histogram specification[END_REF][START_REF] Wan | Joint exact histogram specification and image enhancement through the wavelet transform[END_REF]. The intuition behind these functionals was that their minimizer can up to some degree remove some quantization noise and in this way yield an ordering of the pixels close to the unknown original real-valued image. Such an effect can be observed in Fig. 1 where a synthetic real-valued image is quantized and then restored using the proposed variational method. The nonsmooth L 1 -TV model was originally studied in [START_REF] Chan | Aspects of total variation regularized l 1 function approximation[END_REF]. The main feature of its minimizers is that they contain parts that are equal to the data image and parts that are constant (living in a vanishing component of the TV term). Even though the model modification proposed in [START_REF] Nikolova | Exact histogram specifcation for digital images using a variational approach[END_REF] might seem trivial, the minimizers of these C 2 smoothed ℓ 1 -TV functionals exhibit a qualitatively different behavior. Unlike the L 1 -T V (ℓ 1 -TV) minimizers, it was shown in [START_REF] Nikolova | Exact histogram specifcation for digital images using a variational approach[END_REF] that the minimizers of the C 2 smoothed ℓ 1 -TV functionals generically do not have pixels equal to those of the data image and there are no equally valued pixels. Some of the authors of [START_REF] Nikolova | Exact histogram specifcation for digital images using a variational approach[END_REF] observed that once the parameters of the model were fixed, for all kind of realworld digital images f , the residual error obeyed û-f ∞ = b where the constant b typically met b < 0.5. For this reason, they qualified this variational approach as detail preserving. Therefore we were interested in monitoring the error ûf ∞ .

In this paper we consider a wider class of C 1 smoothed ℓ 1 -TV functionals involving also ℓ 2 data fidelity terms. We give explicit relationships between the minimizers and the parameters tuning the model. The observation that û-f ∞ = b, up to a small difference, is independent of Real-valued image f -quantized on {0, • • • , 15}

Restored û

Figure 1: The restored image is obtained by minimizing J(•, f ) of the form [START_REF] Aubert | Mathematical problems in image processing[END_REF] where ψ(t) = √ t 2 + α 1 and ϕ(t) = √ t 2 + α 2 for N 8.

the input image, is confirmed theoretically. Clear indications on the role of the parameter setting and the lower and upper bounds of ûf ∞ enable us to give restrictions on the parameter selection. All theoretical results are confirmed using numerical tests on a set of digital images of different sizes with disparate content and quality.

In spite of the progress in nonsmooth convex optimization [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF], smooth approximations of nonsmooth objectives still remain a common approach in optimization [START_REF] Beck | Smoothing and First Order Methods: A Unified Framework[END_REF]. Our results can help to design smooth approximations of ℓ 1 /ℓ 2 -TV functionals in a proper way.

The outline of this paper is as follows: In the next Section 2 we describe the variational model. Then, in Section 3 we estimate the ℓ ∞ -error between the input image f and the minimizer of the functional. Section 4 provides explicit parameter estimates for the model. In Section 5 we give probability estimates for the behavior of neighboring pixels. Numerical tests demonstrate the quality of our estimates in Section 6. Finally, Section 7 finishes with conclusions and perspectives.

The Fully Smoothed ℓ 1 -TV Model

We consider M × N digital images f with gray values in {0, . . . , L -1}. Let n := M N . To simplify the notation we reorder the image columnwise into a vector of size n and address the pixels by the index set I n := {1, • • • , n}. Further, we denote by I int n ⊂ I n the subset of all inner pixels, i.e., all pixels which are not boundary pixels.

We are interested in the minimizer û of a functional of the form

J(u, f ) := Ψ(u, f ) + βΦ(u), β > 0 (1) with Ψ(u, f ) := i∈In ψ(u[i] -f [i]), Φ(u) := i∈In j∈N i ϕ(γ i,j (u[i] -u[j])) ,
where N i is a neighborhood of pixel i, the γ i,j > 0 are weighting terms for the distance between neighbors, and the functions ψ and ϕ depend on a positive parameter, α 1 and α 2 , respectively.

To emphasize this dependence we use the notation ψ(•, α 1 ) and ϕ(•, α 2 ) when necessary. So ψ : R × (0, +∞) → R and ϕ : R × (0, +∞) → R. The functions ψ and ϕ have to fulfill the properties stated below:

H0 The functions t → ψ(t, α 1 ) and t → ϕ(t, α 2 ) are continuously differentiable and even.

We denote

ψ ′ (t, α 1 ) := d dt ψ(t, α 1 ) and ϕ ′ (t, α 2 ) := d dt ϕ(t, α 2 ) .
When it is clear from the context, we write ψ ′ (t) for ψ ′ (t, α 1 ) and ϕ ′ (t) for ϕ ′ (t, α 2 ). By H0 , ψ ′ (t) and ϕ ′ (t) are continuous and odd functions.

These derivative functions have to satisfy certain conditions given next.

H1 ψ t → ψ ′ (t, α 1 ) : R → (-Y, Y )
, where Y > 0, is a strictly increasing function for any fixed α 1 ∈ (0, +∞) and maps onto (-Y, Y ).

H2 ψ There is a constant T > 0 such that for any fixed t ∈ (0, T ), the function α 1 → ψ ′ (t, α 1 ) is strictly decreasing on (0, +∞).

Here the cases T = +∞ and Y = +∞ are included.

H1 ϕ t → ϕ ′ (t, α 2
) is an increasing function for any fixed α 2 ∈ (0, +∞) satisfying

lim t→∞ ϕ ′ (t, α 2 ) = 1.
H2 ϕ For any fixed t > 0, the function α 2 → ϕ ′ (t, α 2 ) is continuous and decreasing on (0, +∞) and lim

α 2 ց0 ϕ ′ (t, α 2 ) = 1.
These properties imply further useful relations which are collected in the following remark.

Remark 1 i) By H1 ψ we know that ψ is strictly convex and monotone increasing on (0, +∞) and by H1 ϕ that ϕ is convex. Therefore there exists a unique minimizer of (1). This minimizer can be computed, e.g. by using a Weiszfeld-like semi-implicit algorithm, or the nonlinear (preconditioned) conjugate gradient method, see [START_REF] Chan | On the convergence of the lagged diffusivity fixed point method in total variation image restoration[END_REF][START_REF] Nikolova | Exact histogram specifcation for digital images using a variational approach[END_REF][START_REF] Weiszfeld | Sur le point pour lequel la somme des distances de n points donnés est minimum[END_REF], among other viable algorithms.

ii) By H1 ψ there exists the inverse function (ψ ′ ) -1 (•, α 1 ) : (-Y, Y ) → R, and this function is also odd, continuous and strictly increasing.

Some relevant choices of functions θ obeying all properties H0, H1 ψ , H2 ψ , H1 ϕ and H2 ϕ are given in Table 1. For the latter functions, t → θ ′ (t, α) maps onto (-1, 1), i.e., Y = 1 and T = +∞ for any α > 0. A typical graph of such a function, its derivative and inverse derivative is depicted in Fig. 2.

θ θ ′ (θ ′ ) -1 Θ1 t 2 + α t √ t 2 + α y α 1 -y 2 Θ2 |t| -α log 1 + |t| α t α + |t| αy 1 -|y| Θ3 α log cosh t α tanh t α α atanh(y)
Table 1: Options for functions θ obeying all the assumptions stated above. These functions are nearly affine beyond a neighborhood of zero. The size of the latter neighborhood is controlled by the parameter α > 0.

-3 1, where the plots for α = 0.05 are in blue solid line and for α = 0.5 in red dashed line.

0 3 3 -3 0 3 -1 0 1 -1 0 1 -5 0 5 θ(t) = √ t 2 + α θ ′ (t) = t √ t 2 +α (θ ′ ) -1 (y) = y α 1-y 2
Another choice for ψ fulfilling H0, H1 ψ and H2 ψ is the scaled ℓ p -norm for p = α 1 + 1:

ψ(t) := 1 α 1 + 1 | t | α 1 +1 with ψ ′ (t) = | t | α 1 sign(t), (ψ ′ ) -1 (y) = | y | 1 α 1 , α 1 > 0. ( 2
)
Here ψ ′ maps onto R so that Y = +∞. Moreover α 1 → ψ ′ (t, α 1 ) is strictly monotone decreasing for |t| < 1 hence T = 1 in this case. An upper bound for ûf ∞ when α 1 = 1 in (2) was derived in [START_REF] Nikolova | Analytical bounds on the minimizers of least squares[END_REF]. Some general results on the functionals J for α 1 = 1 in (2) can be found in [START_REF] Aubert | Mathematical problems in image processing[END_REF] in a continuous setting.

For ϕ we can also use the scaled Huber function

ϕ(t) :=      t 2 2α 2 if |t| ≤ α 2 , |t| - α 2 2 if |t| > α 2 with ϕ ′ (t) =    t α 2 if |t| ≤ α 2 , sign(t) if |t| > α 2 . (3)
Note that the functions ψ and ϕ in Table 1 and (3) are nearly affine beyond a small neighborhood of the origin.

In this paper, we focus on the neighborhoods N4 and N8 depicted in Fig. 3 top. When taking the gradient of the functional in [START_REF] Aubert | Mathematical problems in image processing[END_REF] we have to take into account that the pixel combination u[i]u[j] appears for j ∈ N 2 i , where N 2 i denotes the "double" neighborhood associated with N i in Fig. 3 bottom. The usual choices are (see e. g. [START_REF] Geman | Constrained restoration and recovery of discontinuities[END_REF]) γ i,j := 1 for vertical and horizontal neighbors,

γ i,j := 1 √ 2 for diagonal neighbors. (4)
In all cases we have γ i,j = γ j,i .

Functionals of the form (1) with functions ψ, ϕ ∈ C s , s ≥ 2 having alike properties (e.g. all functions in Table 1) were successfully used in [START_REF] Nikolova | Exact histogram specifcation for digital images using a variational approach[END_REF] to process digital images f so that the obtained minimizer û is quite close to the input digital image but its pixels can be ordered in a strict way. An analysis of the minimizers û of these functionals has shown that with a probability close to one, û has pixel values that are different from each other and different from the input pixels.

Bounds for the ℓ ∞ -Error

In this section, we give upper and lower estimates for the ℓ ∞ -error between the input image f and the image û obtained by minimizing the functional J(•, f ). The double neighborhoods N 4 2 and N 8 2 appear in the gradient of Φ(u), see [START_REF] Coltuc | Exact histogram specification[END_REF].

If û is a minimizer of u → J(u, f ) we denote by h ∈ R n the vector with components

h[i] := j∈N 2 i γ i,j ϕ ′ (γ i,j (û[i] -û[j])) , i ∈ I n . (5)
First we provide a lemma which gives a useful expression for ûf ∞ .

Lemma 1 Let H0 , H1 ψ and H1 ϕ be satisfied. Let û be the minimizer of u → J(u, f ) and h be given by [START_REF] Chan | Aspects of total variation regularized l 1 function approximation[END_REF]. Then

û -f ∞ = (ψ ′ ) -1 (β h ∞ , α 1 ) . (6)
Proof. In this proof we can omit the parameter α 1 . Using the definition of J and taking into account that ϕ ′ is odd, we have

∂Ψ ∂u[i] = ψ ′ (u[i] -f [i]) and ∂Φ ∂u[i] = j∈N 2 i γ i,j ϕ ′ (γ i,j (u[i] -u[j])) . (7)
The minimizer û of J(•, f ) has to satisfy ∇ u J(û, f ) = 0 which can be rewritten as

∇ u Ψ(û, f ) = -β∇Φ(û) or as ψ ′ (û[i] -f [i]) = -β j∈N 2 i γ i,j ϕ ′ (γ i,j (û[i] -û[j])), i ∈ I n .
Using (5), the latter is equivalent to

ψ ′ (û[i] -f [i]) = -β h[i], i ∈ I n .
Since ψ ′ is by H0 and H1 ψ odd and strictly increasing,

ψ ′ û[i] -f [i] = ψ ′ (û[i] -f [i]) = β h[i] . (8)
Using Remark 1ii), we see that ( 8) is equivalent to

û[i] -f [i] = (ψ ′ ) -1 β h[i] (9)
where (ψ ′ ) -1 is strictly increasing, hence

û -f ∞ = max i∈In (ψ ′ ) -1 β h[i] = (ψ ′ ) -1 (β h ∞ ) .
For inner points i ∈ I int n we define

η := j∈N 2 i γ i,j . (10)
Of course η does not depend on i but just on the choice of the neighborhood. If the weights are defined as in ( 4), we have

η = 4 for N4 , η = 4 + 4 √ 2 = 6
.8284 for N8 .

For i ∈ I n \ I int n we have

j∈N 2 i γ i,j
≤ η whose value depends on the boundary conditions.

In order to extend the obtained result, we shall use a property of (ψ ′ ) -1 which is stated below.

Lemma 2 Let ψ satisfy H0 , H1 ψ and H2 ψ . Set

Y := min{Y, ψ ′ (T )} , where ψ ′ (T ) := lim t→+∞ ψ ′ (t) if T = +∞. Then for any y ∈ (0, Y ), the function α 1 → (ψ ′ ) -1 (y, α 1 )
is strictly increasing on (0, +∞).

Proof. Let 0 < a 1 < a 2 and y ∈ (0, Y ) be arbitrarily fixed. Since t → ψ ′ (t, α 1 ) is one-to-one and odd, there exist t 1 , t 2 ∈ (0, T ) such that

ψ ′ (t 1 , a 1 ) = y = ψ ′ (t 2 , a 2 ) . (11)
Thus we have (ψ ′ ) -1 (y, a 1 ) = t 1 and (ψ ′ ) -1 (y, a 2 ) = t 2 . ¿From H1 ψ , t → ψ ′ (t, α 1 ) is strictly increasing for any fixed α 1 > 0 and from H2 ψ , α 1 → ψ ′ (t, α 1 ) is strictly decreasing for any fixed t ∈ (0, T ). Therefore

t 2 ≤ t 1 ⇒ y = ψ ′ (t 1 , a 1 ) > ψ ′ (t 1 , a 2 ) ≥ ψ ′ (t 2 , a 2 ).
This contradicts [START_REF] Nikolova | Exact histogram specifcation for digital images using a variational approach[END_REF]. Consequently, t 1 < t 2 which implies the assertion.

For all functions in Table 1 and for ψ in (2) we have Y = 1.

The following theorem provides an upper bound for ||ûf || ∞ which is independent of f as well as of the particular shape of ϕ(t, α 2 ) provided that the latter meets the relevant assumptions.

Theorem 1 Assume that H0 , H1 ψ and H1 ϕ are satisfied. Let β η < Y , where η is given in [START_REF] Nikolova | Analytical bounds on the minimizers of least squares[END_REF].

Then the minimizer û of u → J(u, f ) satisfies û -f ∞ ≤ (ψ ′ ) -1 βη, α 1 =: b(β, α 1 ) . ( 12
)
If, in addition, ψ fulfills H2 ψ and β η < Y , where Y = min{Y, ψ ′ (T )}, then α 1 → b(β, α 1 ) is strictly increasing on (0, +∞).

Proof. From H1 ϕ , ϕ ′ is increasing with |ϕ ′ (t)| ≤ 1 for any t ∈ R. Inserting this into the definition of h in (5) yields h ∞ ≤ η. (13)
Since (ψ ′ ) -1 is by Remark 1ii) strictly increasing on (0, Y ), we deduce from ( 6) and ( 13) for

β η < Y that û -f ∞ = (ψ ′ ) -1 (β h ∞ , α 1 ) ≤ (ψ ′ ) -1 (β η, α 1 ) .
If ψ meets H2 ψ and β η < Y we obtain by Lemma 2 that the function α 1 → (ψ ′ ) -1 βη, α 1 is strictly increasing on (0, +∞).

We clarify the statement of Theorem 1 below.

• By Remark 1, the function β → b(β, α 1 ) is strictly increasing since η is a fixed number.

• The equality in [START_REF] Wan | Joint exact histogram specification and image enhancement through the wavelet transform[END_REF] can only be met if ϕ ′ attains the limit in H1 ψ , i.e., if ϕ ′ (t) = 1 for some t ∈ R. This is for example the case for the scaled Huber function in (3).

• The bound in [START_REF] Wan | Joint exact histogram specification and image enhancement through the wavelet transform[END_REF] depends only on ψ(•, α 1 ) and on β but it is independent of the selection of ϕ provided that H1 ϕ holds.

• For all functions ψ listed in Table 1 we have Y = 1 which limits the action of β to less than 1/η . So H2 ψ furnishes a flexible tool to control the upper bound b(β, α 1 ) by using α 1 under the condition that βη < Y , where we remind that Y = 1 for all ψ in Table 1 and in [START_REF] Beck | Smoothing and First Order Methods: A Unified Framework[END_REF].

The lower bound on ûf ∞ exhibited in the next Theorem 2 depends on ϕ(t, α 2 ) and on the input image f as well. In our formula, the reliance on f is expressed via the magnitude ν f defined below:

I := i ∈ I int n sign f [i] -f [j] = σ, ∀ j ∈ N i where σ ∈ {-1, +1} , ν f := max i∈I min j∈N i γ i,j f [i] -f [j] , (14)
where we set ν f := 0 if I = ∅. The values of ν f for some real-world images can be seen in Fig. 7.

Theorem 2 Let H0 , H1 ψ , H2 ψ and H1 ϕ , H2 ϕ be verified. Let β η < Y , where η is given in [START_REF] Nikolova | Analytical bounds on the minimizers of least squares[END_REF]. Assume that ν f > 2b(β, α 1 ). Then the minimizer û of u → J(u, f ) fulfills

û -f ∞ ≥ (ψ ′ ) -1 (c β η, α 1 ) =: ℓ(β, α 1 , α 2 , ν f ) , (15)
where c = c β, α 1 , α 2 , ν f := ϕ ′ ν f -2b(β, α 1 ), α 2 ≤ 1 . The function α 2 → ℓ(β, α 1 , α 2 , ν f ) is decreasing on (0, +∞) and ℓ(β, α 1 , α 2 , ν f ) ր b(β, α 1 ) as α 2 ց 0 . ( 16
)
Moreover, for ε > 0 arbitrarily close to zero, α 2 can be set so that

û -f ∞ ≥ (ψ ′ ) -1 ((1 -ε) β η, α 1 ) . ( 17
)
Proof. ¿From the definition on ν f , there exists i ∈ I int n such that

γ i,j f [i] -f [j] ≥ ν f , ∀ j ∈ N i .
We consider the case

γ i,j (f [i] -f [j]) ≥ ν f > 2b(β, α 1 ), ∀j ∈ N i . (18)
The opposite case, namely

γ i,j (f [j] -f [i]) ≥ ν f > 2b(β, α 1
), ∀ j ∈ N i can be handled in the same way. By Theorem 1, the minimizer û of J(•, f ) meets

-b(β, α 1 ) ≤ û[i] -f [i], -b(β, α 1 ) ≤ f [j] -û[j], ∀j ∈ N i . Thus -2b(β, α 1 ) ≤ û[i] -û[j] -f [i] -f [j] , ∀j ∈ N i , -2b(β, α 1 ) + f [i] -f [j] ≤ û[i] -û[j], ∀j ∈ N i . (19)
Combining ( 18) and ( 19) along with the fact that γ i,j ≤ 1 yields

0 < -2b(β, α 1 ) + ν f ≤ -2b(β, α 1 ) + γ i,j (f [i] -f [j]) ≤ γ i,j (û[i] -û[j]) ∀ j ∈ N i . Since t → ϕ ′ (t, α 2 ) is increasing by H1 ϕ , the value h[i] in (5) satisfies h[i] ≥ j∈N 2 i γ i,j ϕ ′ ν f -2b(β, α 1 ), α 2 = η c β, α 1 , α 2 , ν f .
Using yet again that y → (ψ ′ ) -1 (y, α 1 ) is strictly increasing (Remark 1ii)) we obtain by ( 9) that

û[i] -f [i] ≥ (ψ ′ ) -1 (c βη, α 1 ) = ℓ(β, α 1 , α 2 , ν f ) . Since û -f ∞ ≥ û[i] -f [i] , it follows that û -f ∞ ≥ (ψ ′ ) -1 (c βη, α 1 ) .
Using H2 ϕ , the function α 2 → c(β, α 1 , α 2 , ν f ) is continuous and decreasing on (0, +∞) and lim

α 2 ց0 c(β, α 1 , α 2 , ν f ) = 1. Combining the latter with Remark 1ii) entails that α 2 → ℓ(β, α 1 , α 2 , ν f )
is decreasing on (0, +∞). Then the definition of b(β, α 1 ) in ( 12) leads to (16).

Finally, H2 ϕ shows that for ε arbitrarily close to zero there is

α 2 > 0 such that c(β, α 1 , α 2 , ν f ) = (1 -ε) and consequently û -f ∞ ≥ (ψ ′ ) -1 ((1 -ε) βη).
Some comments on Theorem 2 may be useful.

• The expression in (17) tells us that by decreasing α 2 , the lower bound ℓ(•) can be adjusted arbitrarily close to the upper bound b(•). The amount of decrease of α 2 needed to reach (1ε) depends on the input image f and can be calculated.

• If t → ϕ ′ (t, α 2
) is nonstrictly increasing on [0, +∞), as the Huber function in (3), it is easy to see that there is α 2 such that c(β, α 1 , α 2 , ν f ) = 1 and hence ℓ(β, α 1 , α 2 , ν f ) = b(β, α 1 ).

Explicit Parameter Estimates

In this section we want to use the error bounds from the previous section to give explicit parameter estimates of β, α 1 and α 2 for the functions ψ, ϕ mentioned in Section 2. More precisely, for a given β satisfying a constraint and for δ fixed, we exhibit the value α 1 = α 1 ensuring that b(β, α 1 ) = δ and then calculate ℓ(β, α 1 , α 2 , ν f ).

For the functions ψ in Table 1 and in [START_REF] Beck | Smoothing and First Order Methods: A Unified Framework[END_REF] we have Y = 1. When the weights γ i,j are chosen as in (4) and H2 ψ holds, the assumption βη < Y = 1 in Theorem 1 reads β < 1 4 = 0.25 for N4, β < 1 6.8284 = 0.1464 for N8.

(20)

In the following we choose β > 0 such that β < 1 η . For δ > 0 fixed, let α 1 solve the equation

b(β, α 1 ) = (ψ ′ ) -1 (βη, α 1 ) = δ . (21)
Then we have by Theorem 1 that ûf ∞ ≤ δ for all α 1 ∈ (0, α 1 ] and there does not exist α 1 > α 1 such that ûf ∞ ≤ δ holds true. In this sense we call α 1 optimal for δ. This claim is ensured thanks to H2 ψ which guarantees that α 1 → b(β, α 1 ) is strictly increasing (see Lemma 2). The value c in Theorem 2 depends on ϕ and on f via ν f . Given the input image f the constant ν f is easy to compute. When

z := ν f -2b(β, α 1 ) > 0, Theorem 2 indicates that the constant c reads c = ϕ ′ (z, α 2 ). (22)
In our experiments on real-world digital images, we always had z ≫ 0 for δ = 0.5. By Theorem 2 a sharper lower bound requires a smaller value for α 2 . According to Theorem 1 and Theorem 2, the upper and lower bounds for fû ∞ and the optimal value for α 1 as defined in (21) for the functions ψ in Table 1 and in (2) are given in Table 2.

ψ(t) b(β, α 1 ) ℓ(β, α 1 , α 2 , ν f ) α 1 t 2 + α 1 α 1 (βη) 2 1 -(βη) 2 α 1 (cβη) 2 1 -(cβη) 2 δ 2 1 β 2 η 2 -1 |t| -α 1 log 1 + |t| α 1 α 1 βη 1 -βη α 1 cβη 1 -cβη δ 1 βη -1 α 1 log cosh t α 1 α 1 atanh(βη) α 1 atanh(cβη) δ atanh(βη) 1 α 1 + 1 | t | α 1 +1 (βη) 1 α 1 (cβη) 1 α 1 ln(βη) ln δ
Table 2: Bounds and parameter α 1 for various functions ψ in Table 1 and in [START_REF] Beck | Smoothing and First Order Methods: A Unified Framework[END_REF]. The parameter c depends on ϕ ′ by (22). The allowed values for β by Theorem 1 are given in (20).

If δ = 0.5 then û has the important property that it preserves the order of the pixel values in a digital image f ∈ {0, . . . , L -1} n . The corresponding values α 1 and β are presented in Table 3 3: Allowed values β < 1/η and the optimal α 1 for δ = b(β, α1 ) = 0.5.

• calculate δ when α 1 and β are given-this can be useful e.g. when ℓ 1 -TV or ℓ 2 -TV are approximated by a fully smooth functional;

• determine the optimal β for fixed α 1 and δ-we remind that from Remark 1, β → b(β, α 1 ) is strictly increasing, hence this value of β is unique.

Probability Estimates for Pixel Neighborhoods

Consider that the assumptions H0, H1 ψ , H1 ϕ and H2 ϕ are met and that the parameters β < Y /η, α 1 and α 2 are fixed. ¿From Theorem 2 we know that the upper bound b(β, α 1 ) in Theorem 1 provides a nearly perfect approximation of the true error û -

f ∞ when c = ϕ ′ (ν f -2b, α 2
) is close to one, which by H1 ϕ means that ν f is large enough. In order to get an intuition-even though very rough-on the behaviour of ν f , we assume in this section that the values of f are realizations of a discrete random variable X taking values in {0, . . . , L -1} whose probability density function (pdf) p X is specialized to real-world digital images. Fig. 4 shows an image together with its histogram which furnishes an empirical estimate of the corresponding pdf. First, we ask for the probability that an inner image pixel i ∈ I int n fulfills

f [i] -f [j] ≥ a and sign(f [i] -f [j]) = σ, ∀j ∈ N i ( 23
)
where σ ∈ {-1, +1} and a > 0 is fixed.

Lemma 3 Let X, X i , i = 1, . . . , k be independent and identically distributed (iid) discrete random variables taking values in {0, . . . , L -1}. Then it holds for a > 0 that q(X, k, a)

:= P (X -X 1 ≥ a, . . . , X -X k ≥ a) = L-1 i=0 (P (X ≤ i -a)) k P (X = i). ( 24
)
Proof. Since the random variables are iid we obtain

P (X -X 1 ≥ a, . . . , X -X k ≥ a) = L-1 i=0 (P (i -X 1 ≥ a, . . . , i -X k ≥ a, X = i) = L-1 i=0 (P (X ≤ i -a)) k P (X = i).
A case relevant to our context is when X is a given inner pixel and X i for i ∈ {1 • • • , k} are the pixels in the "double" neighborhood of X, see Fig. 3. Then the setting of Lemma 3 considers neighborhoods where the central pixel X is bigger than all its neighbors by at least the amount of a. It is clear that the opposite case (when X -X i ≤ -a for all i ∈ 1 • • • , k) is of the same interest and appears with the same probability P (X -X 1 ≤ -a, . . . , X -X k ≤ -a) = q(X, k, a). Of course the "iid" assumption is not realistic for natural images.

For k = 1, the probabilities P (X -X 1 ≥ a) and P (X -X 1 ≤ -a) can be easily exemplified. Let X and X 1 follow independently the same pdf p X . In order to obtain the joint pdf of X and X 1 , one has to compute P (X = i 1 )P (X = i 2 ) for all gray levels i 1 , i 2 obeying |i 1i 2 | ≥ a and then take their sum. Fig. 5 (left) shows for example the joint pdf of X and X 1 when X and X 1 are iid random variables following the pdf p X of the "ducks image" in Fig. 4 left. At position (i 1 , i 2 ) ∈ {0, . . . , 255} 2 the probability P (X = i 1 )P (X 1 = i 2 ) is visualized as a gray value where lighter areas correspond to higher probability.

In Fig. 5 (right) the shaded areas show the points where the pixel difference |i 1i 2 | is larger or equal to a. The sum of the probabilities corresponding to these areas is 2q(X, 1, a). Theorem 3 Assume that the M × N image f is the realization of a discrete iid random vector (X i ) n i=1 with iid components X i as X, where n = M N . Let ν f be defined as in (14) with respect to N4. Then the probability that ν f ≥ a > 0 is not smaller than

1 -(1 -2q(X, 4, a)) m , ( 25
)
where q is defined in (24) and m = ⌊M/3⌋ × ⌊N/3⌋.

For N8 we have to replace q by q(X, 4, a) :=

L-1 i=0 (P (X ≤ i -a)) 4 P (X ≤ i - √ 2a) 4 P (X = i).
Proof. We consider only inner pixels i with non-overlapping neighborhoods as depicted in Fig. 6. Then, by Lemma 3, the probability that one of these pixels does not verify ( 23) is given by 1 -2q(X, 4, a). Hence the probability that all these inner pixels do not fulfill ( 23) is (1 -2q(X, 4, a)) m and the probability that at least one of these pixel satisfies (23

) is 1 -(1 -2q(X, 4, a)) m .
Note that for q(X, 4, a) > 0 the probability in ( 25) is indeed very close to 1 even for moderate sizes of m. For instance, if the random variables are uniformly iid, we have

q(X, 4, a) = 1 L L-1 i=a i -a + 1 L 4 . = (L -a)(L -a + 1)(2(L -a) + 1)(3(L -a) 2 + 3(L -a) -1) 30L 5 .
For a = 137 and L = 256 this formula gives q(X, 4, a) ≈ 0.0044 and for M = N = 128 further 1 -(1q(X, 4, a)) m ≈ 1 -10 -7 . We tested two functionals J(•, f ) as described in Section 2: the first corresponds to ψ = Θ1 and ϕ = Θ1 and the second to ψ = Θ2 and ϕ = Θ1 as given in Table 1. In all tests, N8 was adopted with the weights γ i,j given in [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF]. Two choices for β satisfying (20) were considered along with different values for α 1 and α 2 . The minimizers û were computed using Polak-Ribière conjugated gradients [START_REF] Bonnans | Numerical Optimization (Theoretical and Practical Aspects)[END_REF] with high numerical precision. For each restored image we computed û-f ∞ and present the distance between the theoretical upper bound b(β, α 1 ) and the obtained û

-f ∞ : b(β, α 1 ) -û -f ∞ .
The tables show also the difference between the upper and the lower theoretical bounds on ûf ∞ : bℓ := b(β, α 1 )ℓ(β, α 1 , α 2 , ν f ) , computed using the explicit formulae given in Section 4. Furthermore, we evaluate the amount of pixels that closely approach the ℓ ∞ norm: Over the whole set of these images, for α 2 = 0.02 we have mean 0.5ûf ∞ = 2.968 × 10 -6 and mean 0.5ℓ(β, α 1 , α 2 , ν f ) = 6.0678 × 10 -6 . For α 2 = 100 these values read mean 0.5ûf ∞ = 1.307 × 10 -2 and mean 0.5

q = ♯ i ∈ I n û -f ∞ -| û[i] -f [i] |
-ℓ(β, α 1 , α 2 , ν f ) = 2.491 × 10 -2 .
where ♯ stands for cardinality and ε 0 in order to account for numerical errors. In the experiments, we set ε := 10 -3 .

In all tests, given 0 < β < 1/η, we fixed α 1 = α 1 so that b(β, α 1 ) = δ for δ = 1 2 .

The numerical outcomes confirm the theoretical results on û-f ∞ established in Sections 3 and 4. ¿From Tables 4, 5 and 6 the following observations can be drawn:

• Decreasing α 2 > 0 towards 0 enables to make the difference between the upper and the lower bounds on ûf ∞ arbitrarily small which leads to uf ∞ ≈ b(β, α 1 ).

In this case a large percentage of the pixels i meet |û

[i] -f [i]| ≈ b(β, α 1 ).
• An important increase of α 2 > 0 entails a decrease of the lower bound ℓ(β, α 1 , α 2 , ν f ). Moreover, the number of pixels i verifying |û

[i] -f [i]| ≈ b(β, α 1
) is reduced to a few ones.

Such a situation may be preferable when one wishes that there are not too many pixels close to the upper bound.

Tables 7 and8 show yet again that the gap between the upper bound b(β, α 1 ) and the lower bound ℓ(β, α 1 , α 2 , ν f ) vanishes when α 2 is close to zero and that it increases when α 2 increases. For α 2 fixed, we see that b(β, α 1 )ℓ(β, α 1 , α 2 , ν f ) tends to decrease along with β.

Fig. 8 shows the histograms of the differences {f [i]û[i], i ∈ I n } relevant to "moon", where the upper bound was set to b(β, α 1 ) = 0.5, for an increasing set of values of α 2 . These histograms were plotted for 100 bins equally spaced in [-0.5, +0.5]. For very small values of α 2 , there are many pixels meeting |f For α 2 = 0.02 we have mean 0.5ûf ∞ = 1.777 × 10 -6 and mean 0.5ℓ(β, α 1 , α 2 , ν f ) = 3.666 × 10 -6 . For α 2 = 100, we find mean 0.5ûf ∞ = 8.265 × 10 -3 and mean 0. 6: Results for ψ = Θ2, ϕ = Θ1, β = 0.05 and a small and large value of α 2 , respectively. For α 2 = 0.05 we have mean 0.5ûf ∞ = 5.441 × 10 -6 and mean 0.5ℓ(β, α 1 , α 2 , ν f ) = 10.29 × 10 -6 . For α 2 = 100, we find mean 0.5ûf ∞ = 1.09 × 10 -2 and mean 0.5ℓ(β, α 1 , α 2 , ν f ) = 2.11 × 10 -2 .

[i] -û[i]| ≈ f -û ∞ .
7 Conclusions and Open Questions ℓ 1 -TV and ℓ 2 -TV functionals have been often minimized using a smoothed version of the form we consider in this paper with ad hoc chosen smoothing parameters ("very small"). The results established in our work enable to clearly evaluate the resulting approximation. The functions (ψ, ϕ) studied here have a lot of similarities. However, they produce different image restorations. The question of what couple of functions (ψ, ϕ) would give a better result in the framework of a given application, remains open.

Extension to the rotational-invariant (in a discrete sense) smoothed TV, i.e. Φ(u) = i,j ϕ( ∇ i,j u ), where ∇ i,j u ∈ R 2 stands for a discrete approximation of the gradient of u at pixel (i, j), deserves attention.

Extensions to cases when f are the coefficients of the expansion of the input image using an orthogonal transform as the discrete cosine transform or a frame transform as the curvelet transform, see, e.g., [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] are of interest.

Applications to quantization noise reduction should be envisaged.

Figure 2 :

 2 Figure 2: The function Θ1 in Table1, where the plots for α = 0.05 are in blue solid line and for α = 0.5 in red dashed line.

N 4 N 8 Double N 4 2

 482 Double N 8 2

Figure 3 :

 3 Figure 3: Neighborhoods N 4 (left) and N 8 (right) of a pixel (i, j) are used to formulate Φ(u).The double neighborhoods N 4 2 and N 8 2 appear in the gradient of Φ(u), see[START_REF] Coltuc | Exact histogram specification[END_REF].

Figure 4 :

 4 Figure 4: Left: Duck image. Right: Histogram of "duck image" furnishing an empirical estimate of the corresponding pdf.

Figure 5 :

 5 Figure 5: Left: Joint pdf of two iid random variables X, X 1 where X and X 1 follow the pdf of the "ducks image" in Fig. 4 right. Here light areas correspond to high probability. Right: Areas where |i 1i 2 | ≥ a, i 1 , i 2 ∈ {0, . . . , L -1}. The value 2q(X, 1, a) is the sum of the probabilities in the shaded areas.

Figure 6 :

 6 Figure 6: Disjoint 3 × 3-adjacencies with center pixels "x".

2 + α 2 α 2 =

 222 5ℓ(β, α 1 , α 2 , ν f ) = 1.610 × 10 -2 . N8, ψ(t) = |t|α 1 log 1 + |t| α 1 for α 1 = 0.9645, β = 0.05, hence b = 0.5, ϕ(t) = √ t b-û-f ∞ | ×10 -6 b-ℓ | ×10

.

  Remark 2 Equation (21) offers several other exploits than only fixing the optimal α 1 . For any

			ψ(t)	√	t 2 + α 1	|t| -α 1 log 1 +	|t| α 1
			neighborhood	β	α 1	β	α 1
			N4 N4 N8 N8	0.2 0.1406 0.1 1.3125 0.1 0.2862 0.05 1.8947	0.2 0.1250 0.1 0.7500 0.1 0.2322 0.05 0.9645
			Table
	β <	Y η	one can also

Table 4 :

 4 Figure 7: The set of images used in the tests provided in this section. The values of ν f are computed according to (14) in the case N8 for the weights in (4). Results for ψ = Θ1, ϕ = Θ1, β = 0.1 and a small and large value of α 2 , respectively.

	< ε	and Q% = 100	q n	,

Table 5 :

 5 When α 2 increases, such pixels become more and more rare and the differences |f [i]û[i]| become centered near zero. However they never reach zero: see the value of µ defined in the caption of the figure. Here again, the numerical tests were done with a high precision. Results for ψ = Θ1, ϕ = Θ1, β = 0.05 and a small and large value of α 2 , respectively.

	N8, ψ(t) =	√	t 2 + α 1 for α 1 = 1.895, β = 0.05 hence b = 0.5, ϕ(t) = α 2 = 0.02 α 2 = 100 √ t 2 + α 2
	image chemical	b-û-f ∞ | ×10 -6 b-ℓ | ×10 -6 Q % 2.561 9.055 4.54	b-û-f ∞ | ×10 -3 b-ℓ | ×10 -2 q 14.17 3.993 2
	moon	1.580		3.300	10.2	7.649	1.572	1
	aerial	0.872		2.093	3.92	4.229	1.015	2
	bark	1.673		4.254	6.82	8.239	2.000	1
	couple	1.642		3.432	3.25	7.830	1.632	4
	motioncar	12.39		20.35	0.28	51.43	7.847	1
	stream	0.727		1.240	7.19	3.291	0.608	3
	tank	1.020		1.701	8.31	5.678	0.829	1
	man	0.162		0.374	6.00	0.968	0.186	11
	Pentagon	0.871		1.442	10.2	3.877	0.706	1
	clock	1.013		2.220	2.88	4.193	1.073	1
	boat	0.799		2.795	7.14	5.136	1.342	2
	tree	0.993		2.009	6.06	4.895	0.975	2
	brick wall	0.125		0.329	11.9	1.115	0.164	99
	airplane	0.228		0.403	3.48	1.274	0.200	1

Table

	0.101	0.304	-6 Q % 2.79	b-û-f ∞ | ×10 -3 b-ℓ | ×10 -2 q 18.81 5.236 2
	5.347	11.06	7.03	10.22	2.090	1
	2.670	7.019	2.63	5.663	1.354	2
	5.843	14.26	5.55	11.01	2.653	1
	5.369	11.51	3.25	10.46	2.170	4
	41.36	68.23	0.09	66.99	0.101	1
	1.687	4.155	6.66	4.404	0.813	3
	3.869	5.703	4.45	7.592	1.107	1
	0.673	1.255	3.14	1.298	0.249	10
	2.723	4.837	6.55	5.188	0.943	1
	2.622	7.437	2.88	5.610	1.431	1
	3.879	9.373	3.97	6.874	1.786	2
	4.070	6.737	4.18	6.549	1.301	2
	0.721	1.102	11.3	1.710	0.219	61
	0.682	1.352	0.74	4.983	0.268	1

mean b(β, α 1 )ℓ(β, α 1 , α 2 , ν f ) , b(β, α 1 ) = 0.5, N8 α 2 = 0.01 α 2 = 100 The mean value of the difference b(β, α 1 )ℓ(β, α 1 , α 2 , ν f ) was computed over the selection of images shown in Fig. 7. Here we consider the N8 neighborhood for the weights in [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF].

Table 8: The neighborhood here is N4 with the weights given in (4). The mean is calculated over the set of images in Fig. 7.


Nina Amenta
email: amenta@ucdavis.edu

U C Davis

Dominique Attali
email: dominique.attali@gipsa-lab.inpg.fr.

Olivier Devillers

A tight bound for the Delaunay triangulation of points on a polyhedron

We show that the Delaunay triangulation of a set of n points distributed nearly uni-

), where k = ⌈ d+1 p+1 ⌉. This bound is tight in the worst case, and improves on the prior upper bound for most values of p.

Introduction

Overview. The Delaunay triangulation of a set of points is a fundamental geometric data structure, used in surface reconstruction, mesh generation, molecular modeling, geographic information systems, and many other areas of science and engineering. It is well-known [START_REF] Mcmullen | The maximum number of faces of a convex polytope[END_REF] that the complexity of the Delaunay triangulation of n points in di-

mension d is O(n ⌈ d 2 ⌉
) and that this bound is achieved by distributions of points along one-dimensional curves such as the moment curve. But points distributed uniformly in R d , for instance inside a d-dimensional ball, have Delaunay triangulations of complexity O(n); the constant factor is exponential in the dimension, but the dependence on the number of points is linear. We are interested in filling in the picture in between these two extremes, that is, when the points are distributed on a manifold of dimension

2 ≤ p ≤ d -1.
In this paper we consider the easy case of a p-dimensional polyhedron P .

Main result. We consider a fixed p-dimensional polyhedron P in d-dimensional Euclidean space R d . Our point set S is a sparse ε-sample from P . Sparse ε-sampling requires the sampling to be neither too sparse nor too dense. Our sampling model also contains the restriction that every face of P must be ε-sampled, not just those of the highest dimension. This kind of sampling is used in low dimensions for mesh generation, eg. [START_REF] Cheng | A practical Delaunay meshing algorithm for a large class of domains[END_REF].

We consider the complexity of the Delaunay triangulation of S, as ε → 0, while P remains fixed, so that n = |S| goes to infinity. The main result in this paper is that the number of simplices of all dimensions is O(n d-k+1 p

) where k = ⌈ d+1 p+1 ⌉. The hidden dimensional manifold is O( √ n). Erickson proved that the Delaunay triangulation of a set of points in R 3 with spread ∆ is O(∆ 3 ). Perhaps even more interestingly, he

showed that this bound is tight for ∆ = √ n, by giving an example of a sparse εsample of points from a cylinder that has a Delaunay triangulation of size Ω(n 3/2 ) [START_REF] Erickson | Nice point sets can have nasty Delaunay triangulations[END_REF].

This surface is not generic and has a degenerate medial axis.

To the best of our knowledge, our earlier abstract [START_REF] Amenta | Complexity of Delaunay triangulation for points on lower-dimensional polyhedra[END_REF] is the only prior result for

d > 3.
Outline of the proof. We begin with a simple lower bound, which illustrates the importance of the quantity k = ⌈ d+1 p+1 ⌉ and motivates the upper bound.

At the coarsest level, the proof of the upper bound is similar to that of [START_REF] Amenta | Complexity of Delaunay triangulation for points on lower-dimensional polyhedra[END_REF]: we map Delaunay simplices to the medial axis in such a way that each medial point can only correspond to a constant number of Delaunay simplices, and then we use a packing argument on the medial axis. But the tight bound seems to require us to consider some additional phenomena which occur in higher dimensions, but not in the familiar settings of dimensions two and three. The most important of these is the observation that when k = ⌈ d+1 p+1 ⌉ > 2, the vertices of any Delaunay simplex, which must span R d , have to be drawn from k > 2 faces of P . This allows us to relate Delaunay simplices to only the lower-dimensional parts of the medial axis, generated by k or more faces of Pthe k-medial axis. This idea is embodied in Corollary 6. Because the k-medial axis is lower-dimensional, we can pack it with fewer balls, giving the tighter bound.

The mapping process, described in Section 6, takes a Delaunay ball Σ to a point z associated with the k-medial axis. We show in Section 7 that all of the original sample points inducing Σ lie near at least one of the projections of z onto the flats spanned by the faces of P ; we say z covers the sample points. This covering allows us to bound the number of Delaunay simplices by the size of the packing on the medial axis.

To motivate some of the complications in the mapping, let us sketch a simpler version which produces a point on the k-medial axis but not necessarily one which covers the sample points. We say that a ball is hollow if its interior does not intersect the polyhedron P . We begin with the center z of Σ, and the largest hollow ball B centered at z; B is tangent to at least one face of P . We choose a hyperplane H containing all of the faces to which B is tangent. Fixing the intersection B ∩ H, we move z in a straight line away from H, increasing its radius as necessary, until B becomes tangent to another face of P . We iterate, stopping when we run out of degrees of freedom; the final z must lie on the k-medial axis.

But notice that as z moves, the part of B on one side of H grows, while the other part shrinks. The difficulty is that B might shrink away from some of the original sample points inducing Σ, and as B recedes the projections of z onto faces of P may move away as well. In our earlier abstract [START_REF] Amenta | Complexity of Delaunay triangulation for points on lower-dimensional polyhedra[END_REF], we started with a single tangent point and used only one iteration (going only to the 2-medial axis), so that the shrinking part was entirely within distance ε of the first point of tangency. But during the additional iterations required for the tight bound the multiple points of tangency can be far apart from each other, and the shrinking part of B can be quite large.

To get around this problem, we define in Section 4 a new structure related to the medial axis, the annular k-medial axis. An annulus in R d is the symmetric difference between two concentric balls. An annulus is tangent to the boundary of P when either the inner bounding sphere or the outer bounding sphere of the annulus is tangent to the boundary of P . It is hollow when the inner bounding sphere is the boundary of a hollow ball. Just as the medial axis consists of the centers of hollow balls tangent to the boundary of P at more than one point, the annular medial axis consists of the centers of hollow annuli tangent to the boundary of P at more than one point.

At the beginning of the mapping process, we start with a hollow annulus with same center z as Σ, which encloses Σ and is tangent to P in at least one point. In each step of the iterative process, now, we move z until either the inner sphere or the outer sphere of the annulus becomes tangent to a new face. This prevents the spheres bounding the annulus from shrinking away from the faces of P containing the sample points inducing the original Delaunay ball Σ. This new mapping procedure produces a point z on the annular k-medial axis.

The annular k-medial axis is typically unbounded, so to apply a packing argument we need to define (in Section 5) a bounded subset of it, the trimmed annular k-medial axis. By definition, this object has dimension at most dk + 1 and we prove that its ). Since our point set S is a sparse ε-sample from a p-dimensional polyhedron, its cardinality is n = Ω(ε -p ).

(d -k + 1)-dimensional
Eliminating ε gives the O(n

d-k+1 p
) upper bound.

We carefully structure the arguments so that we can avoid making any non-degeneracy assumptions on either P or the vertex set S of the Delaunay triangulation. This means that parts of the annular k-medial axis of a degenerate polyhedron may have dimension greater than dk + 1 (which is what one expects in the generic case). Nonetheless we can show that any simplex of the Delaunay triangulation is mapped to a part of the annular k-medial axis which does have dimension at most dk + 1, by mapping the simplex to an annulus tangent to at least k independent faces (Definition 4).

Lemma 5 shows that k independent faces generate a piece of the correct dimension, and

Lemma 14 shows that a simplex can be mapped to an annulus tangent to k independent faces.

Statement of Theorem

To formally state the theorem, we must first define the sampling condition.

Sampling

Let P ⊆ R denote our input polyhedron, not necessarily convex or connected 1 The density condition implies that all faces, of all dimensions, are nearly uniformly sampled, not just the faces with maximal dimension. We treat λ as a constant. As ε goes to zero, the number n of points in a λ-sparse ε-sample of a p-dimensional polyhedron is related to ε by n = Θ(ε -p ), i.e. n tends to infinity (to see this, consider any packing of every face of P by balls of radius ε, and then extend this packing to a covering by increasing the radius of each ball to 2ε).

Main Theorem

We are now ready to state our main result:

Theorem 1 Let S be a λ-sparse ε-sample of a p-dimensional polyhedron P in R d ,
with n = |S|. In the worst case, the Delaunay triangulation of S has size Θ(n

d-k+1 p )
where k = ⌈ d+1 p+1 ⌉.

The lower bound

We begin with the lower bound, which is comparatively simple and conveys the intuition as to why Θ(n

d-k+1 p
) is the correct bound.

Proof of Lower Bound for Theorem 1: We first construct a polyhedron P , and then examine the complexity of the Delaunay triangulation of a λ-sparse ε-sample on P .

We select a set of d + 1 affinely independent points, and partition them into groups

Q 1 , . . . , Q k so that groups Q 1 . . . Q k-1 contain p + 1
i , for 1 ≤ i ≤ k, producing a set X.
The point set X is linearly independent, since the C i span independent subspaces, so the convex hull of X is a (k -1)-simplex σ.

Also, for any X that we can construct this way by picking one sample point s i on each C i ,, there is a unique (d -1)-sphere Σ tangent to P at the s i ∈ X. To see this, notice that being tangent to C i at s i is roughly speaking equivalent to passing through

c i + 1 coincident points; and k i=1 (c i + 1) = d + 1. (1)
Σ encloses no sample point of S in its interior, since it is tangent to each of the C i at a single point. Thus σ must be a Delaunay simplex.

Since C i contains Ω(ε -ci ) points of S, the number of distinct Delaunay (k -1)simplices that we can construct is at least

Ω(ε -c1 × • • • × ε -c k ) = Ω(ε -(d-k+1) ) = Ω(n d-k+1 p ).
and this is a lower bound on the overall complexity of the Delaunay triangulation.

The annular medial axis

The annular medial axis AM is the key geometric object in our proof of the upper bound. The annular medial axis is formed by tangent annuli rather than tangent spheres.

Setting k = ⌈ d+1 p+1 ⌉, we are particularly interested in AM k , the part of the annular medial axis generated by at least k faces of P . If the faces of P were in general position, it would be fairly simple to establish that AM k had dimension at most dk + 1

everywhere, but when P is degenerate this is not the case; see Figure 2. Instead, we identify the subset of AM k which is guaranteed to have dimension at most dk + 1.

Definitions

It is clear what it means for a sphere to be tangent in the interior of an open face F of P ; let us introduce a definition which also handles the boundaries. We denote the closure of F by cl(F ) and write Aff(F ) for the affine space spanned by F . We say that a (d -1)-sphere Σ is tangent to F at point x if both cl(F ) and Aff(F ) intersect Σ in a unique point x. Since faces are relatively open, this means that x might be a limit point of F , so that Σ ∩ F = ∅. Note also that a sphere can be tangent to several faces at x, only one of which is the face F x containing x (the one of smallest dimension).

An annulus with center z, inner radius r and outer radius R is the set of points x whose distance to the center satisfies r ≤ xz ≤ R. The boundary of an annulus consists of two (d -1)-spheres and we call the smaller one the inner sphere and the larger one the outer sphere. We say that an annulus A is tangent to F at x if one of the two spheres bounding A is tangent to F at x (see Figure 1, or Figure 6 for a more intricate example). Point x is called a tangency point of A. An annulus is hollow if its inner sphere bounds a d-ball whose interior does not intersect P ; notice that P might intersect the annulus itself. The width of an annulus is the difference between the outer and inner radii Rr. We define the penetration of an annulus to be the maximum radius of the intersection of the annulus with a flat, such that the open ball bounded by the inner sphere remains empty; see Figure 1. Specifically, the penetration

is √ R 2 -r 2 .
We are mostly concerned with annuli of penetration ε; such an annulus has width ε 2 R+r .

Definition 2

The annular k-medial axis AM k of P is the set of points z ∈ R d for which the largest hollow annulus of penetration ε centered at z is tangent to at least k faces of P .

We write A(z, ε) for the largest hollow annulus of penetration ε centered at z.

Dimension and degeneracy

To describe the annular medial axis, we use the language of stratifications [START_REF] Goresky | Stratified Morse theory[END_REF][START_REF] Gunnells | Stratified spaces twigs[END_REF].

Given a subset X ⊆ R d , a stratification of X is a filtration ∅ = X -1 ⊆ X 0 ⊆ • • • ⊆ X j = X
by subspaces such that the set difference X i \ X i-1 is an open i-dimensional manifold (possibly not connected), called the i-dimensional stratum S i of X. For example, the Voronoi diagram of a point set in R 2 admits a stratification into its cells, edges and vertices. Semi-algebraic sets admit a stratification, and since the annular k-medial axes of polyhedra are piecewise semi-algebraic, they also admit a stratification. In this section, we give conditions under which a point z ∈ AM k belongs to a stratum of dimension dk + 1 or less.

If the faces of P were in general position, this would be easy: any z belonging to the stratum of dimension d-k +1 would be the center of an annulus tangent to k faces.

But in degenerate situations an annulus tangent to k faces might belong to a stratum of dimension greater than dk + 1, as illustrated in Figure 2. So we need to make a more careful argument.

Recall that AM k is defined as the set of points for which A(z, ε) is tangent to at least k faces. Let us concentrate instead in this section on the subset C k , defined as the set of points z such that annulus A(z, ε) is tangent to exactly k faces of P . When does

C k form a (d -k + 1)-dimensional stratum?
We start by writing down some equations that determine C k locally around z. Since A(z, ε) is tangent to exactly k faces F 1 , . . . , F k , there exists δ > 0 such that every face of the polyhedron not in {F 1 , . . . , F k } is at distance at least δ to the boundary of A(z, ε). Using a compactness argument as in [START_REF] Lieutier | Any open bounded subset of has the same homotopy type as its medial axis[END_REF], it follows that for a point y close enough to z, the only faces that might possibly be tangent to annulus A(y, ε) are F 1 , . . . , F k . We set e i = -ε 2 if F i is tangent to the outer sphere of A(z, ε) and e i = 0 if F i is tangent to the inner sphere of A(z, ε). Since any annulus A(y, ε) with y ∈ C k in the neighborhood of z is tangent to both F i and F k , C k thus obeys the following k -1 equations:

g i (y) = d(y, F i ) 2 -d(y, F k ) 2 + e i -e k = 0, for 0 < i < k.
The equations expressing the fact that A(y, ε) is tangent to any other pair F i , F j are redundant. Each g i (y) is a polynomial of degree two, so every neighborhood in C k is a subset of the intersection of k -1 quadrics. In general, k -1 hypersurfaces meet at point z in a (dk + 1)-manifold, but this might not hold in the presence of degeneracies.

Lemma 3 Suppose z ∈ C k is the center of an annulus A(z, ε) tangent to the polyhedron at k affinely independent points x 1 , . . . , x k . Then the neighborhood of z in C k is a (dk + 1)-manifold. Furthermore, the tangent space to C k at z is the set of vectors orthogonal to the affine space spanned by x 1 , . . . , x k .

PROOF. In a small neighborhood of z, C k coincides with the zero-set of the functions g i (y). Let F i = F xi be the face to which x i belongs. One can check that the gradient of g i (y) is 2(x kx i ). Since the x i are affinely independent for 1 ≤ i ≤ k, we have that the Jacobian of the map g = (g 1 (y), . . . , g k-1 (y)) has rank k -1, so the implicit function theorem implies that its zero-set g -1 (0) is a (dk + 1)-dimensional manifold in a neighborhood of z. Furthermore, the tangent space of g -1 (0) at z is the null space of the derivative Dg(z), which is the set of vectors orthogonal to the affine space spanned by x 1 , . . . , x k .

While the assumption that the tangency points of A(z, ε) are independent is sufficient to show that z belongs to the stratum of dimension at most d-k+1, this condition is not necessary. Again, this is a phenomenon that occurs only in high dimension. As an example of an annulus A(z, ε) for which the tangency points are not independent but the dimension of the stratum containing z remains small, let's consider the following example in R 6 . The polyhedron P consists of segments of four lines, represented parametrically by (0, 0, s, 0, 0, 0), (1, 0, 0, t, 0, 0), (0, 1, 0, 0, u, 0), and (1, 1, 0, 0, 0, v), with -1 ≤ s, t, u, v ≤ 1. Let us take ε = 0, and consider the part of C 4 consisting of points z with A(z, ε) tangent to the interiors of all four line segments. This is a three-dimensional manifold, as one would expect, and at most points z in this portion of AM 3 the four points of tangency span a three-dimensional subspace. But the point z = (1/2, 1/2, 0, 0, 0, 0) also belongs to this part of C 4 and its points of tangency all lie in a common 2-flat. Nonetheless z is a regular point of C k , and C k has dimension three in the neighborhood of z. This example illustrates why we need to consider not just the local points of tangency in the neighborhood of z but also the tangent faces.

Definition 4

We say that k faces F 1 , . . . , F k are independent if none of them is contained in the affine space spanned by the union of the others, that is for

1 ≤ i ≤ k, F i ⊆ Aff(F 1 ∪ • • • ∪ F i ∪ • • • ∪ F k ),
where the symbol over F i indicates that it is omitted in the union.

For example, the four segments in R 6 which we just considered are independent. PROOF. Recall that C k is the set of points z ∈ R d such that A(z, ε) is tangent to exactly k faces of P . We partition C k into k pieces (some of which might be empty), as follows. We write S i for the set of points y ∈ C k whose tangency points span a space of dimension i. Thus we have C k = i S i , for 0 ≤ i ≤ k -1. All we need to prove is that in the neighborhood of any z ∈ S i , S i has dimension at most dk + 1, for all 0 ≤ i ≤ k -1. By Lemma 3, we already know that S k-1 is a (dk + 1)dimensional manifold, so let us consider i < k -1. The idea of this proof is to define S i , in the neighborhood of z, as the intersection of a family of surfaces, and to compute the dimension of the intersection. Let F 1 , . . . , F k be the k faces tangent to A(z, ε). Using the same compactness argument as before, there exists a small neighborhood U of z in R d such that for every point in w ∈ S i ∩ U , the only faces possibly tangent to the annulus A(w) are F 1 , . . . , F k .

For y ∈ S i ∩U , let x i (y) be the orthogonal projection of y onto Aff(F i ) (see Figure 3). We can always reduce the size of U such that the tangency points x 1 (y), . . . , x k (y)

span an affine space of dimension i; without loss of generality, let us assume that the first i tangency points x 1 (y), . . . , x i (y) are affinely independent. The first surface we construct is the locus of annulus centers tangent just to the faces F 1 , . . . , F i . Specifically, let

P ′ = cl(F 1 ) ∪ • • • ∪ cl(F i )
and write S ′ for the set of points which are the center of a hollow annulus (with respect to P ′ ) of penetration ε tangent to the i faces F 1 , . . . , F i . By Lemma 3, S ′ is a (di + 1)-manifold in a neighborhood of y. Clearly

S i ⊆ S ′ .
Now, for each of the other F j , we construct another surface containing y. For i < j ≤ k, the point x j (y) is an affine combination of x 1 (y), . . . , x i (y) and therefore belongs to

X j = Aff(F 1 ∪ • • • ∪ F j ∪ • • • F k-1 ).
Since F j ∩ X j contains the tangency point x j (y), it is not empty. We define H j as the set of points w ∈ R d such that the nearest point to w on Aff(F j ) lies in Aff(F j ∩ X j ), that is, the set of points that, like y, fall into Aff(F j ∩ X j ) when projected to Aff(F j ).

H j is an affine space, orthogonal to F j , of dimension ddim

F j + dim(F j ∩ X j ). So all y ∈ C k in the neighborhood of z lie in S ′ ∩H i+1 ∩H i+2 • • •∩H k . Let us prove
that this intersection has dimension at most d-k+1. In Lemma 3, we observed that the normal space to S ′ at y is spanned by the i -1 vectors v 2 = x 1 (y)x i (y), . . . , v i =

x i-1 (y)x i (y). For i + 1 ≤ j ≤ k, we can always find a vector v j in the normal space to H j , by choosing v j in the tangent plane to F j and orthogonal to F j ∩ X j .

Because the F j are independent faces, the vectors v 2 , . . . , v k are linearly independent and all belong to the normal space of the intersection

S ′ ∩ H i+1 ∩ H i+2 ∩ • • • ∩ H k . It
follows that the intersection is a manifold of dimension at most dk + 1, and that S i is a stratified space of dimension at most dk + 1.

We deduce immediately the following corollary:

Corollary 6 Let z ∈ AM k and suppose that A(z, ε) is tangent to j faces amongst which k faces are independent. Then, z lies on a stratum of AM k with dimension at most dk + 1.

18

We have established that the subset of AM k induced by independent faces has the "right" dimension dk + 1, which is necessary to our plan of counting Delaunay simplices via a packing argument. It may, however, extend to infinity and hence have infinite (dk + 1)-dimensional volume. In this section we define a bounded (dk + 1)-dimensional subset of AM k , which we call the trimmed annular k-medial axis T AM k . We end this section by proving that the volume of T AM k is bounded by a constant that does not depend on ε.

We need some definitions. We say that a hyperplane supports X ⊆ R d if it has nonempty intersection with the boundary of X and empty intersection with the interior of X. We begin our proof that T AM k is bounded and that this bound does not depend on ε with two simple, and similar, technical lemmata. Lemma 9 Let A be an annulus tangent to a hyperplane H at point x and whose inner sphere does not enclose point q, with q and the center of A on the same side of H. Let R and r be respectively the outer and inner radii of A. Suppose that there exist two scalars D and h > 0 such that d(q, H) ≥ h, xq ≤ D and Rr ≤ h 2 . Then, the inner radius of A satisfies r ≤ D 2 h .

PROOF. We first consider what happens when the inner sphere of A passes through q and point x lies on the outer sphere of A (see Figure 5). Let y be the intersection of the inner sphere of A with the segment connecting x to the center z of A. Let c be the midpoint of the segment yq. Since the angle between the two vectors cz and qz is equal to the angle between the vector qy and the hyperplane H, we have

q -y 2r = d(q, H) -(R -r) q -y
The bound on r follows immediately.

For the case in which x lies on the inner sphere of A, apply the same argument, assuming that r = R and x = y.

Lemma 10 Let A(z, ε) be tangent to two hyperplanes H i , H j , at points x i , x j such that their normal vectors v i , v j , pointing towards z, have angle difference ∠(v i , v j ) > α. Let D be the diameter of P . Then the inner radius r of A is at most D sin α/2 .

PROOF. If x i is tangent to the larger sphere of A(z, ε), we let y i be the intersection point of segment x i , z with the inner sphere; similarly we define y j if x j is tangent to the outer sphere. In any case, we consider the distance D ′ between the two points on the inner sphere, and note that D ′ ≤ D. We have sin α/2 = D ′ /2, so that r =

d(y, z) ≤ D 2 sin α/2 .

21

We plan to show that the radius of any essential annulus can be bounded using one or the other of these technical lemmata. Let F = {F i } i∈I be the set of faces of P tangent to A(z, ε). Each F i is tangent to A(z, ε) at a point x i , and we define the normal v i to A(z, ε) at x i to be the unit vector with direction of zx i , i.e.

v i = (z -x i )/||z -x i ||.
We will choose a constant α later; it will depend only on P and not on z or ε. If there is a pair of normals v i , v j such that ∠(v i , v j ) > α, we bound the radius of A(z, ε)

using Lemma 10, with D taken to be the diameter of P . So it only remains to bound the size of an essential annulus, for which every pair of normals differs by at most α, which is done in Lemma 12.

We begin by characterizing an annulus as essential or non-essential using only its set F of tangent faces, irrespective of z and ε.

Let S d-1 = {v ∈ R d | v = 1}
δ Fi (v) = max{ q -x, v | ∀x ∈ cl(F i ), ∀q ∈ P }.
Roughly speaking, δ Fi (v) represents the distance between an extreme point in direction v on P and an extreme point in direction -v on the closure of F i . Notice that because

F i ⊂ P , δ Fi (v)
≥ 0 for all v, since zero is achieved for any x ∈ F i by choosing q = x.

In addition, we have δ Fi (v) = 0 exactly when every point x ∈ F i is an extreme point of P in direction v, that is, when there is a supporting plane of the convex hull of P containing F i . Finally, we observe that δ Fi (v) is continuous, and, since it is defined on the compact space S d-1 , it is uniformly continuous (by the Heine -Cantor Theorem).

22

We now consider any set F = {F i } i∈I of faces. Let

δ F (v) = max i∈I δ Fi (v).
As the maximum of a set of continuous non-negative functions, this is again a continuous non-negative function.

Observation 11

The annulus A(z, ε) with tangent face set F is non-essential if and only if there exists a unit vector v such that δ F (v) = 0.

PROOF. We have δ F (v) = 0 if and only if δ Fi (v) = 0 for all i. This happens when the tangent plane to P in direction v contains all of the F i ∈ F. This means, by definition, that A(z, ε) is non-essential.

Thus, since A(z, ε) is essential, we must have δ F (v) > 0 for all unit vectors v.

Since the map δ F is continuous and defined on a compact set, it attains a global minimum and this minimum is positive. We define

2h = min F min v δ F (v),
where v ranges over all unit vectors and F ranges over all subsets of faces such that F is not contained in a hyperplane supporting the convex hull of P . We have h > 0. We use this constant, dependent on P , in the following theorem.

Lemma 12 There exists a constant α, dependent on the geometry of P but not on ε, such that for any essential annulus A(z, ε) for which every pair v i , v j of its normal vectors has angle distance ∠(v i , v j ) < α, the distance d(z, P ) < D 2 h , where D is the diameter of P .

PROOF. From the fact that A(z, ε) is essential, we can conclude that for any of the

F i ∈ F, there is some v j such that δ Fi (v j ) > 2h. If δ Fi (v i ) > 2h, we are done, but
this is not guaranteed to be the case. So now we use the fact that the lemma assumes that the directions of the v j are all similar. Since δ Fi is uniformly continuous, we can select an α Fi (again depending on F and not on z or ε), such that

∠(v i , v j ) < α Fi =⇒ |δ Fi (v i ) -δ Fi (v j )| < h.
We define α = min F α F over all faces F of the polyhedron.

With this choice of α, the fact that

δ Fi (v j ) > 2h implies that δ Fi (v i ) > 2h-h = h.
This in turn implies that there exists a point q ∈ P such that d(q, H i ) ≥ h, where H i is the hyperplane containing F i with normal v i . So we can apply Lemma 9 to bound the inner radius of A(z, ε).

Combining Lemma 10 and Lemma 12, we conclude that the set of essential points is bounded and at distance at most max D sin α/2 , D 2 h from P . Let B(P ) be the smallest ball containing the set of essential points. The T AM k is a stratified set, of dimension

at most d -k + 1, in B(P ). Its i-dimensional stratum, for 1 ≤ i ≤ d -k + 1, is the
union of pieces of semi-algebraic sets, each of dimension i, formed by the intersection of at least k polynomials, each of degree at most two. The number of such pieces is bounded by a function of P (corresponding to the choices of at least k faces producing the polynomials), that is independent of ε. Each piece is also trimmed by B(P ). The idimensional volume of the intersection of a ball of radius R with an algebraic variety of dimension i formed by m polynomials of bounded degree is bounded by a function of i, m and R (see [START_REF] Demmel | The probability that a numerical analysis problem is difficult[END_REF], Lemma 7.4). So the overall i-dimensional volume of the stratum is at most the sum of these trimmed varieties, and hence is itself bounded by some function of the geometry of P , independent of ε.

Mapping Delaunay spheres to T AM k

The goal of this section is to assign every Delaunay sphere Σ to some point z on the trimmed annular k-medial axis. We begin with a quick geometric lemma.

Lemma 13 A Delaunay sphere Σ with center z 1 is contained in the annulus A(z 1 , ε) .

PROOF. Let x be a point in P with minimum distance to z 1 and let s be the sample

point in S ∩ cl(F x ) closest to x. The distance x -s ≤ ε and therefore s ∈ A(z 1 , ε). Because Σ encloses no sample point, Σ ⊆ A(z 1 , ε).
It is the preceding simple lemma which requires the restrictive sampling requirement that the faces of all dimensions must be sampled; if F x was a subface of some larger-dimensional face F , and not itself sampled, a hollow annulus with center z 1 could contain a large portion of F x , extending far away from the projection x of z 1 onto F x and F . While we could find a sample s on F within distance ε of x, s would be contained in an annulus of width ε, but not however in the annulus A(z 1 , ε) of penetration ε (recall Section 4.1).

We use the following incremental construction to associate an annulus, called Expansion(Σ),

with each Delaunay sphere Σ; we will then prove that the center z of Expansion(Σ) lies in T AM k , and that Σ ∩ P is contained in a set that we call the cover of z, the definition of which depends on d, p and ε.

Recall that we assume that P spans R d (Section 2.1), and that p is the dimension of P . Initially, we let z 1 be the center of Σ. The inner sphere B(z 1 , r) of the annulus A(z 1 , ε) is tangent to P in at least one point x, contained in a face F 1 = F x . At each step j of the construction we perturb the annulus to find another tangent face, and we increase the dimension of the subspace spanned by {F 1 . . . F j } by at least one.

At an arbitrary step of the construction, we have a set {F 1 , . . . , F j } of j faces tangent to A(z j , ε). If the dimension of Aff({F 1 , . . . , F j }) < d, we can find a hyperplane H that contains their union j i=1 F i (possibly there are many such H; we can pick one arbitrarily). Let L be the line orthogonal to H and containing the current z j . Consider any annulus A(y, ε) centered on L in the neighborhood of z j . Since A(y, ε) also has penetration ε, its intersection A(y, ε) ∩ H = A(z j , ε) ∩ H, and A(y, ε) remains tangent to F 1 , . . . , F j . Consider the y, nearest to z j , such that A(y, ε) is tangent to an additional face F j+1 , not contained in H. Such a y must exist, since P spans R d . We set z j+1 = y.

When there is a set of faces F 1 , . . . , F j tangent to A(z j , ε) that spans R d , we stop.

Clearly this occurs after at most d steps. We associate Σ with the final annulus: define Expansion(Σ) = A(z j , ε), for the final z j .

Lemma 14 For every Delaunay sphere Σ, the center z of the annulus Expansion(Σ) belongs to the trimmed annular k-medial axis T AM k .

PROOF. Since we know that there is a set of faces {F 1 . . . F j } tangent to A(z, ε) = Expansion(Σ) that span R d , we know that they cannot all be contained in a hyperplane supporting the convex hull, and this implies that z is essential. The crux of the proof is showing that A(z, ε) is tangent to at least k independent faces. This will establish that z belongs to the (dk + 1)-dimensional stratum of AM k , and hence to T AM k .

Initially, since any single face is independent, the set {F 1 } is independent. The dimension of Aff({F 1 . . . F j+1 }) is at least one greater than the dimension of Aff({F 1 . . . F j }),

since Aff(F 1 . . . F j ) ⊂ H, and F j+1 is not contained in H. This also implies that at each step F j+1 is independent of F 1 . . . F j . However, it is possible that some F i in {F 1 . . . F j } may be spanned by some subset of {F 1 . . . F j+1 }, for instance if

F i ⊆ F j+1 . Any such F i can be removed, with Aff({F 1 . . . F j+1 }-F i ) = Aff({F 1 . . . F j+1 }).
Let F be the remaining set of independent faces.

We claim that at the end of the construction, there are at least k independent faces in F, for k = ⌈ d+1 p+1 ⌉. Consider adding each F i to the F in turn. Each face has dimension at most p, so that adding F i increases the dimension of the union {F 1 . . . F i-1 } by at most p + 1 (adding to the subspace basis a vector from the previous subspace Aff(F 1 . . . F i-1 ) to some point in F i , and then a set of at most p vectors spanning F i ).

So the number of independent faces required to span all of d-dimensional space is at

least k = ⌈ d+1 p+1 ⌉.
This establishes that the final A(z, ε) is tangent to at least k independent faces.

Covering Delaunay simplices

Now we would like to relate the intersection Expansion(Σ) ∩ P back to the original Delaunay sphere Σ. We do this by defining the cover of a point z (roughly its region of influence on P ), and then showing that the intersection Σ ∩ P is contained in the cover of the annulus Expansion(Σ) centered at a point z ∈ T AM k .

So lets begin by defining the cover of a point z. Write π x (z) for the orthogonal projection of z onto the flat supporting the face F x containing x in P and let Π(z) = {π x (z) | x ∈ P } be the set of orthogonal projections of z onto the flats supporting all faces of P . For any non-negative number w, called the radius of the cover, we define the w-cover as a union of balls:

Cover(z, w) = x∈Π(z) B(x, w).
We say that x is an anchor point of z if x ∈ A(z, ε) and π x (z) = x. Thus, anchor points of z form a subset of Π(z) that contains the tangency points of the annulus with the polyhedron and possibly other points of P ∩ A(z, ε). Now we establish that any point in P ∩ A(z, ε) must be close to an anchor point of z. Lemma 15 For every point x ∈ P ∩ A(z, ε) there exists an anchor point y of z in the closure of F x such that

x -y ≤ (dim F x -dim F y + 1) ε.
PROOF. The proof is by induction over the dimension d x = dim F x of the face F x containing x. If d x = 0, the result holds for y = x. Suppose d x > 0 and let q = π x (z) be the orthogonal projection of z onto Aff(F x ). We distinguish two cases:

(1) if q ∈ F x , the segment xq lies inside A(z, ε) and therefore xq ≤ ε; and the result holds for y = q. (2) If q ∈ F x , we consider the point y ∈ P on the segment xq, which is closest to x on the boundary of F x (as in Figure 6). Since the segment xy is contained in A(z, ε), this implies xy ≤ ε. Furthermore, since y belongs to the boundary of F x , the dimension d y < d x . We apply our induction hypothesis: cl(F y ) contains an anchor point w of z, with yw ≤ (dim F ydim F w + 1) ε, so

xw ≤ ε + (dim F ydim F w + 1) ε, which gives the result.

The preceding Lemma shows that any point in P ∩ A(z, ε) is at distance at most (p + 1)ε to an anchor point of z. Thus, P ∩ A(z, ε) ⊆ Cover(z, (p + 1)ε). From there, we show how to cover a Delaunay sphere using an annulus centered on the trimmed medial axis.

Lemma 16 For every Delaunay sphere Σ, we have Σ ∩ P ⊆ Cover(z, 4dε) where z is the center of Expansion(Σ).

PROOF. Again, this follows from the construction of A(z j , ε) = Expansion(Σ).

Initially, z 1 is the center of Σ. By Lemma 13, Σ ∩ P ⊂ A(z 1 , ε), and by Lemma 15, this implies that for every point x 0 ∈ Σ ∩ P , there exists an anchor point x 1 of z such that x 0 -x 1 ≤ (dim F x0 -dim F x1 +1)ε. In particular, Σ∩P ⊆ Cover(z 1 , (p+1)ε). Now we need to show that as we proceed from z i to z i+1 , Σ ∩ P remains covered for 1 ≤ i < j. So we construct a sequence of j points x 1 , x 2 , . . . , x j such that x i+1 is an anchor point of z i+1 and

D D ′ H - H + H - H + A(y, ε) L + H H xi A(zi, ε) A(zi+1, ε) Aff Fi zi y zi+1 xi
x i -x i+1 ≤ (dim F xi -dim F xi+1 + 2)ε, (2)
for all 1 ≤ i < j. For simplicity, we write F i = F xi and d i = dim F i . Let H be the hyperplane used during the construction of z i+1 and which contains the union i l=1 F l .

Let H + be the closed half-space bounded by H in the direction z i+1z i , with H - being the open half space on the other side. We consider two cases:

1. The original anchor point x i ∈ H + , including the case in which x i ∈ H. Let R i be the outer radius of A(z i , ε). The intersection H + ∩ B(z i , R i ) ⊆ H + ∩ B(z i+1 , R i+1 ), and x i is of course outside both inner balls, so x i remains in the annulus A(z i+1 , ε). It may no longer be an anchor point of z i+1 , but by Lemma 15, there exists an anchor point x i+1 of z i+1 on the closure of F i such that x ix i+1 ≤ (d id i+1 + 1) ε. To summarize, for every point x 0 ∈ Σ ∩ P , we have constructed a sequence of points

x 1 , x 2 , . . . , x j such that x i is an anchor point of z i for all 1 ≤ i ≤ j, x 0 -x 1 ≤ (d 0 -PROOF. Recalling that π x (z) is the orthogonal projection of z onto the tangent plane to P at x, we have π x (z)π x (z ′ ) ≤ zz ′ ≤ ε. The claim follows immediately.

Using Lemma 17 and Lemma16, we get that for d ≥ 1 there exists a point z ′ ∈ M such that

Σ ∩ P ⊆ Cover(z ′ , 5dε)
The cover of z ′ is a union of d-balls, each with radius 5dε, at most one for each face of the polyhedron, and therefore it contains a constant number of points of S. So the number of simplices that we can form by picking points in the cover of z ′ is also a constant. This means that only a constant number of Delaunay simplices can be charged to each point z ′ ∈ M , and the size of the Delaunay triangulation is proportional to the size m of M .

Recall that the number of points in a λ-sparse ε-sample S of a p-dimensional polyhedron P is n = Θ(ǫ -p ) and that the i-faces of P have Θ(ε -i ) points of S. Using n = Ω(ε -p ), we get that the number of Delaunay simplices is

O(m) = O(ε -(d-k+1) ) = O(n d-k+1 p ),
where k = ⌈ d+1 p+1 ⌉.

Conclusion

This paper answers only the first of many possible questions about the complexity of the Delaunay triangulations of points distributed nearly uniformly on manifolds. It would be interesting to establish a bound which does not require all faces of P to be sampled. Similar bounds for smooth surfaces rather than polyhedra would be of more practical interest. The proof in this paper seems to rely critically on some properties specific to polyhedra, particularly that sample points on k faces are needed to form a simplex, so other techniques will be needed for the cases of more general manifolds of dimension 1 < p < d.

  volume is bounded from above by a constant that does not depend on ε. It follows that we can construct an ε-sample M of the trimmed annulark-medial axis with m = O(ε -(d-k+1) ) points.The mapping process in Section 6 produces a point z on trimmed annular k-medial axis. We define the cover of z in Section 7 asx∈Π(z) B(x, r),where Π(z) is the set of all orthogonal projections of z onto the flats spanned by faces of P , and B(x, r) is the ball centered at x with radius r. We then assign z to the nearest sample z ′ ∈ M , and argue that the cover of a point z ′ ∈ M , with r = 5dε, contains only a constant number of points in S. Since each z ′ is charged for a constant number of Delaunay simplices. It follows that the size of the Delaunay triangulation is bounded from above by the size of M , which is m = O(ε -(d-k+1)

Figure 1 Figure 1 :

 11 Figure1pictures an example of the annular 2-medial axis in R 2 . Observe that this is a superset of the medial axis, even for ε = 0: the medial axis of the polyhedron is the set of points z ∈ R d for which A(z, 0) touches the polyhedron in two points or more, while the annular 2-medial axis with ε = 0 is the set of points z for which A(z, 0) is tangent to two faces of P or more (possibly at the same point).

Figure 2 :

 2 Figure 2: The box is P , and we take ε = 0. The set C 4 ⊂ AM 4 , consisting of the centers of annuli tangent to four faces is the segment connecting the hollow dots. It has dimension one, rather than dimension zero as we would expect in the absence of degeneracy.

Lemma 5

 5 Suppose that A(z, ε) is tangent to exactly k faces. If those k faces are independent, then C k is a manifold of dimension at most dk + 1 in a neighborhood of z.

Figure 3 :

 3 Figure 3: Notation for the proof of Lemma 5. A(y, 0) is tangent to F1, F2 and Fj.

Definition 7 ADefinition 8

 78 point z ∈ AM k , and the corresponding annulus A(z, ε) are nonessential if there exists a hyperplane supporting the convex hull of P and containing all faces tangent to A(z, ε). Otherwise, the point z and the corresponding annulusA(z, ε) are essential.An annulus A(z, ε) may be essential because it is tangent to a face F not on the boundary of the convex hull, and/or because the union of faces tangent to A(z, ε) spans R d ; see Figure4. Since we assume that P spans R d (Section 2.1), the set of essential points is non-empty. The trimmed annular k-medial axis, T AM k , is the set of essential points lying on the i-dimensional strata of the annular k-medial axis for i ≤ dk + 1.

Figure 4 :

 4 Figure 4: A polyhedron formed of four faces and its annular 2-medial axis, taking ε to be zero. The set of essential points forms the closed piece of parabola consisting of points equidistant to F0 and F3. The two endpoints are essential because the spheres centered at the endpoints are also tangent to F3.

Figure 5 :

 5 Figure 5: Notation for the proof of Lemma 9.

  be the space of direction vectors in R d . To each face F i of the polyhedron P , we associate the function δ Fi : S d-1 → R which maps every unit vector v ∈ S d-1 to

F3Figure 6 :

 6 Figure 6: Notation for the proof of Lemma 15. The annulus is tangent to the four faces F0, F1, F2 and F3.

Figure 7 :

 7 Figure 7: On the upper left, three annuli centered at zi, y and zi+1 that share the same intersection with a hyperplane H. On the lower right, intersection of the three annuli with the tangent space passing through xi. The restriction of those intersections to H -are nested.

2 .

 2 The case in which x i ∈ H -is illustrated in Figure7. Let us begin by consideringthe intersections D = Aff(F i ) ∩ A(z i , ε) and D ′ = Aff(F i ) ∩ A(z i+1 , ε). Since x i is an anchor point contained in A(z i , ε), the intersection D is a d i -dimensionalball of radius at most ε, with center x i . The annulus is shrinking on this side of H, implying D ′ ∩ H -⊆ D ∩ H -. It is not possible for D ′ ∩ cl(F i ) to be empty, since that would require there has been some last point at which the annulus A(y, ε), with y on the segment connecting z i and z i+1 , was in contact with cl(F i ), contradicting the choice of z i+1 as the annulus with the first new tangent point. Since D ′ ∩ cl(F i ) is non-empty, it must contain some point p and by Lemma 15, there exists an anchor point x i+1 of z i+1 on the closure ofF i such that px i+1 ≤ (d id i+1 + 1) ε. Since p ∈ D ′ ⊆ Dand D is a ball centered at x i with radius at most ε, we have x ip ≤ ε. The triangle inequality then gives us Inequality 2.

  points each, and Q k contains between 1 and p + 1 points. Thus it must be that k = ⌈ d+1 p+1 ⌉.Let C i be the convex hull of Q i , and let c i be the dimension of C i . For all 1 ≤ i < k, this dimension c i = p, and c k ≤ p. We let the polyhedron P =

	consider choosing one sample point from each of the C
	k i=1 C i ; note that
	P indeed has dimension p. Also note that each C i spans a subspace independent of the
	other C j .
	Now let us consider a λ-sparse ε-sample of P ; let us call this point set S. Now

More formally, a p-dimensional polyhedron can be defined as the underlying space of any geometric simplicial complex of dimension p.

Acknowledgments

We would like to dedicate this paper to the late Hazel Everett. We thank Hazel, Sue Whitesides, and Sylvain Lazard for inviting us to their 2006 Bellairs Workshop on Computational Geometry, where we started to work on this problem. We gratefully acknowledge André Lieutier for helpful suggestions, and the insights and suggestions of our patient anonymous reviewers. The work of the first author was supported by NSF awards CCF-0331736 and CCF-0635250.

Supported by the EU under STREP contract FET-255827 (CGLearning) 1

Since the number of steps in the construction is bounded by d, i.e. j -1 ≤ d, we get x 0x j ≤ (3d + 1)ε yielding the result.

Size of Delaunay triangulation

To establish the upper bound on the number of Delaunay simplices, it remains only to combine what we know about the size of T AM k with our method of mapping Delaunay simplices to T AM k .

We first consider a sample M of the trimmed annular k-medial axis T AM k , such that every point x ∈ T AM k has a sample within distance ε, and such that the number of samples in a ball of radius O(ε) centered at x is at most λ (unlike our original sample S of P , it is not necessary to sample the lower-dimensional strata of T AM k ).

p+1 ⌉. This follows from the fact that the dimension of T AM k is dk + 1, and the results in Section 5, which established that the volume of T AM k is bounded by a constant that does not depend on ε.

Next, we map each Delaunay simplex σ ∈ Del(S) to a point z ∈ M . Consider the Delaunay sphere Σ passing through the vertices of σ. Lemma 14 and Lemma 16 tell us that the Delaunay sphere Σ associated with σ belongs to Cover(z, 4dε) for some point z ∈ T AM k , where z is an arbitrary point not belonging to M . But there must be some z ′ ∈ M at distance at most ε from z.

Lemma 17 For every pair of points z and z ′ with zz ′ ≤ ε:

Cover(z, wε) ⊆ Cover(z ′ , (w + 1)ε).