Adjoint method in PDE-based image compression

Abstract

We consider a shape optimization based method for finding the best interpolation data in the compression of images with noise. The aim is to reconstruct missing regions by means of minimizing a data fitting term in an $L^{p}$ -norm, for $1 ⩽ p < + \infty$ , between original images and their reconstructed counterparts using linear diffusion PDE-based inpainting. Reformulating the problem as a constrained optimization over sets (shapes), we derive the topological asymptotic expansion of the considered shape functionals with respect to the insertion of small ball (a single pixel) using the adjoint method. Based on the achieved distributed topological shape derivatives, we propose a numerical approach to determine the optimal set and present numerical experiments showing the efficiency of our method. Numerical computations are presented that confirm the usefulness of our theoretical findings for PDE-based image compression.

Keywords

Image compression shape optimization adjoint method image interpolation inpainting PDEs image denoising

1. Introduction and related works

PDE-based methods have attracted growing interest by researchers and engineers in image analysis field during the last decades [8,9,12,17–19,25,27,30,31]. Actually, such methods have reached their maturity both from the point of view of modeling and scientific computing allowing them to be used in modern image technologies and their various applications. Image compression is one of the domain where they appear among the state-of-the-art methods [3,4,11,13,20,24,28]. In fact, the aim for such problems is to store few pixels of a given image (coding phase) and to recover/restore the missing part in an accurate way (decoding). The PDE-based methods use a diffusion differential operator for the inpainting of missed parts from an available data (boundary or small parts of the initial image) therefore their efficiency for decoding is guaranteed/encoded in the operator without any pre- or post-treatment. The question then is how to ensure with these methods a good choice, if it exists, of the “best” pixels to store for high quality reconstruction of the entire image? An answer to this question is given in [6,7] for the harmonic or the heat equation, where its reformulation as a constrained (shape) optimisation problem permitted to exhibit an optimal set of pixels to do the job. In addition, analytic selection criteria using topological asymptotics were derived. Due to the simple structure of the shape functionals considered in these previous works, the topological expansion is easily derived (more or less with formal computations) and gives an analytic criterion to characterize the optimal set in compression. The limitation in obtaining the topological expansion this way is twofold: the criterion gives pointwise information on the importance of the location (pixel) to store which results in hard thresholding selection strategy not robust with respect to the noise. Second, the technique is limited to simple functionals, namely an $L^{2}$ data-fitting term and a linear diffusion operator.

The main contribution of this article is the use of the adjoint method [2–11,13,14] to derive an analytic criterion for PDE-based compression. In fact, we consider the compression problem in the same framework than [6], but we introduce a new approach to the characterization of the set of pixels to select using the adjoint method [2,5,14,15]. This approach to obtain the topological expansion is more general than the one previously studied for the same problem in [6], in the sense that it may be used for other diffusion operators and nonlinear data-fitting term, moreover it allows a better stability with respect to noise. In particular, when the accuracy of the reconstruction (fidelity term) is measured with an $L^{p}$ -norm, $1 ⩽ p < + \infty$ , the adjoint method is still linear and no significant complexity or cost are added. Thus, the main results in the article include the rigorous derivation of the topological expansion based on the adjoint problem in the spirit of [2,14]. We notice that the Dirichlet boundary condition in the inclusion prevents from a direct transposition of the method based on a local perturbation of the material properties by inserting small holes. Therefore, we adapt the sensibility analysis to the problem under consideration and we perform the asymptotic expansion of the proposed shape functional using non-standard perturbation techniques combined with truncation techniques. The asymptotic allows us to deduce a gradient algorithm for the reconstruction that we implement and compare to previous works [6,7].

The article is organized as follows: in Section 2, we introduce the compression problem that takes the form of a constrained optimization problem of finding the best set of pixels to store, denoted K. Section 3 is devoted to describe the adjoint method to compute the topological derivative of the cost functional considered. In Section 4, we perform the computations to obtain the topological expansion and the “shape” derivatives which involve the direct and adjoint states. Finally, in Section 5, we describe the resulting algorithm and we give some numerical results to confirm the usefulness of the theory.

2. Problem formulation

Let $D \subset R^{2}$ and $f : D \to R$ , a given image in some region $K \subset\subset D$ . We consider the mixed elliptic boundary problem for a given $u_{0}$ in $L^{2} (D)$ ,

Problem 2.1
Find u in $H^{1} (D)$ such that
$\begin{aligned} {\begin{array}{ll} u - α Δ u = u_{0} & i n D ∖ K, \\ u = f, & i n K, \\ \frac{\partial u}{\partial n} = 0, & o n \partial D . \end{array} \end{aligned}$
(1)
where the available data f is a Dirichlet “boundary” condition and with homogeneous Neumann boundary condition on $\partial D$ . This PDE corresponds to the first term in the time discretization of the homogeneous heat equation, where we assume that the initial condition is $u_{0}$ . For compatibility condition with the “boundary” data on K, we take as $u_{0}$ the image $f \in H^{1} (D)$ , with $Δ f \in L^{2} (D)$ and such that $\frac{\partial f}{\partial n} = 0$ on $\partial D$ . In the compression step (coding phase), the datum f is available in the entire domain D, so we can set the initial condition $u_{0}$ to the function f. The result of this coding step consists of a set K of the pixels to store and the values of f on K. In the decompression step (decoding phase), the data is only available in the subset K of the domain, so we set $u_{0} = 0$ (at least in $D ∖ K$ ). When the reconstruction is performed by solving the heat equation, it means that we start with an initial datum which do not satisfy the compatibility conditions, but this does not influence the dynamic as far as the convergence to an equilibrium, that is a steady state, holds (regularizing effect). Setting $v = u - f$ , we can write equivalently
Problem 2.2.
Find v in $H^{1} (D)$ such that
$\begin{aligned} {\begin{array}{ll} - α Δ υ + υ = α Δ f & i n D ∖ K, \\ υ = 0, & i n K, \\ \frac{\partial υ}{\partial n} = 0, & o n \partial D . \end{array} \end{aligned}$
(2)

Denoting by $v_{K} = u_{K} - f$ the solution of Problem 2.2, the question is to identify the region K which gives the “best” approximation $u_{K}$ , in a suitable sense, that is to say which minimizes some $L^{p}$ -norm. The constrained optimization problem for the compression reads [6], for $1 ⩽ p < + \infty$ ,
$\begin{aligned} min_{K \subseteq D, m (K) ⩽ c} {\frac{1}{p} \int_{D} | u_{K} - f |^{p} d x | u_{K} solution of Problem 2.1}, \end{aligned}$
(3)
where m is a “size measure”. The optimization problem (3) is studied in [6] and the existence of an optimal set is established for $p = 2$ and m is the capacity of sets [32] and the result extends to $1 ⩽ p < + \infty$ as noticed in [7]. The optimal set K is obtained via a relaxation procedure but its regularity is not considered yet, nevertheless the relaxation technique allows us to derive first order optimality conditions via topological derivatives, which was done in the case of the Laplacian as inpainting operator. In this article, we aim to compute the topological gradient [2,10] of the shape functional using the adjoint method which possesses two main advantages on the previous approaches: it is more general and systematic with respect to the inpainting operator and the exponent $1 ⩽ p < + \infty$ , on one side and secondly, it leads to a better characterization of the relevant pixels as it gives a distribution of such pixels taking into account local information from their neighborhood. Loosely speaking, to obtain the topological derivative, let $x_{0} \in D$ and $K_{ε} = K \cup \bar{B (x_{0}, ε)}$ ( $B (x_{0}, ϵ)$ denotes the ball centered at $x_{0}$ with radius $ε$ ), then we look for an expansion of the form
$\begin{aligned} J (u_{K_{ε}}) - J (u_{K}) = ρ (ε) G (x_{0}) + o (ρ (ε)) . \end{aligned}$
where ρ is a positive function going to zero with $ε$ and G is the so-called topological gradient [2,5,14]. Therefore, to minimize the cost functional J, one has to create small holes at the locations $x \in D$ where $G (x)$ is the most negative. For the compression problem this amounts to select the locations where the pixels are the most important to keep.
3. Adjoint method and variations of the cost functional

The adjoint method has been extensively studied and successfully applied to many second order elliptic problems and Helmholtz equation (see [2,14] and references therein). We will recall the main principle (theorem) of the method and apply it to our specific setting. We introduce the following abstract result which describes the adjoint method for the computation of the first variation of a given cost functional (see for instance [2]). Let V be a Hilbert space. For $ε \in [0, ζ]$ , $ζ > 0$ , we consider a symmetric bilinear form $a_{ε} : V \times V \to R$ and a linear form $l_{ε} : V \to R$ such that the following assumptions are fulfilled

$| a_{ε} (v, w) | ⩽ M_{1} ‖ v ‖ ‖ w ‖$ , $\forall (v, w) \in V \times V$ (continuity of the bilinear form),

$a_{ε} (v, v) ⩾ ξ ‖ v ‖^{2}$ , $\forall v \in V$ (uniform coercivity),

$| l_{ε} (w) | ⩽ M_{2} ‖ w ‖$ , $\forall w \in V$ (continuity of the linear form),

with

α, M_{1}, M_{2} > 0

independent of

ε

. Moreover, we suppose that there exists a continuous bilinear form

δ a : V \times V \to R

, a continuous linear form

δ l : V \to R

and a function

ρ : R_{+} \to R_{+}

such that, for all

ε ⩾ 0

$‖ a_{ε} - a_{0} - ρ (ε) δ a ‖_{L_{2} (V)} = o (ρ (ε))$ ,

$‖ l_{ε} - l_{0} - ρ (ε) δ l ‖_{L (V)} = o (ρ (ε))$ ,

$lim_{ε \to 0} ρ (ε) = 0$ .

We emphasize that

δ a

and

δ l

do not depend on

ε

. Finally, for all

ε \in [0, ζ]

, consider a functional

J_{ε} : V \to R

, Fréchet-differentiable at the point

v_{0}

. Assume further that there exists a number

δ J (v_{0})

such that

\begin{aligned} J_{ε} (w) - J_{0} (v) = D J_{0} (v) (w - v) + ρ (ε) δ J (v) + o (‖ w - v ‖ + ρ (ε)), \forall (v, w) \in V \times V . \end{aligned}

Then we have [2]

Theorem 3.1.
Let $v_{ε} \in V$ be the solution of the following problem: find $v \in V$ such that,
$\begin{aligned} a_{ε} (v, φ) = l_{ε} (φ), \forall φ \in V . \end{aligned}$
Let $w_{0}$ be the solution of the so-called adjoint problem: find $w \in V$ such that
$\begin{aligned} a_{0} (w, φ) = - D J_{0} (v_{0}) φ, \forall φ \in V . \end{aligned}$
Then,
$\begin{aligned} J_{ε} (v_{ε}) - J_{0} (v_{0}) = ρ (ε) (δ a (v_{0}, w_{0}) - δ l (w_{0}) + δ J (v_{0})) + o (ρ (ε)) . \end{aligned}$

To be more specific, for $x_{0} \in D$ and $r > 0$ , we denote by $B_{r}$ the open ball centred at $x_{0}$ and of radius r. We set
$\begin{aligned} V_{ε} := {v \in H^{1} (D ∖ B_{ε}) | v = 0 ~on~ \partial B_{ε}} . \end{aligned}$
Then we consider the boundary value problem:
Problem 3.1.
Find ${\tilde{v}}_{ε}$ in $V_{ε}$ such that
$\begin{aligned} {\begin{array}{ll} - α Δ {\tilde{υ}}_{ε} + {\tilde{υ}}_{ε} = h & i n D ∖ B_{ε}, \\ {\tilde{υ}}_{ε} = 0 & i n B_{ε}, \\ \partial_{n} {\tilde{υ}}_{ε} = 0 & o n \partial D . \end{array} \end{aligned}$
(4)
with $h := α Δ f$ , but h can be any $L^{2} (D)$ function. We denote ${\tilde{v}}_{0}$ the solution of the problem
Problem 3.2.
Find ${\tilde{v}}_{0}$ in $H^{1} (D)$ such that
$\begin{aligned} {\begin{cases} - α Δ {\tilde{v}}_{0} + {\tilde{v}}_{0} = h, & i n D, \\ \partial_{n} {\tilde{v}}_{0} = 0, & o n \partial D . \end{cases} \end{aligned}$
(5)

The weak formulation of problems above reads, find ${\tilde{v}}_{ε}$ in $V_{ε}$ such that, for all φ in $V_{ε} \cap H^{1} (D)$ , we have
$\begin{aligned} {\tilde{a}}_{ε} ({\tilde{υ}}_{ε}, φ) = {\tilde{l}}_{ε} (φ), \end{aligned}$
with,
$\begin{aligned} {\tilde{a}}_{ε} ({\tilde{υ}}_{ε}, φ) & := α \int_{D ∖ B_{ε}} \nabla {\tilde{υ}}_{ε} \cdot \nabla φ d x + \int_{D ∖ B_{ε}} {\tilde{υ}}_{ε} φ d x, \\ {\tilde{l}}_{ε} (φ) & := \int_{D ∖ B_{ε}} h φ d x, \end{aligned}$

The dependency of the space $V_{ε}$ on $ε$ prevents us from using Theorem 3.1 directly, therefore, we introduce a truncation technique [15], which consists of inserting a ball $B_{R}$ , for a fixed $R > ε$ and splitting Problem 3.1 into two sub-problems that we glue at their common boundary (see Fig. 1). More precisely, we consider the sub-problems: an internal problem
$\begin{aligned} {\begin{cases} - α Δ v_{ε, R} + v_{ε, R} = h, & in B_{R} ∖ B_{ε}, \\ v_{ε, R} = 0, & on \partial B_{ε}, \\ v_{ε, R} = v_{ε}, & on \partial B_{R}, \end{cases} \end{aligned}$
and an external problem
$\begin{aligned} {\begin{cases} - α Δ v_{ε} + v_{ε} = h, & in D ∖ B_{R}, \\ \partial_{n} v_{ε} = \partial_{n} v_{ε, R}, & on \partial B_{R}, \\ \partial_{n} v_{ε} = 0, & on \partial D . \end{cases} \end{aligned}$

Fig. 1.
Illustration of the splitting.

As the two sub-problems transform the initial one into a transmission problem. Then,
Proposition 3.1.
We have,
$\begin{aligned} {\tilde{υ}}_{ε} = {\begin{array}{ll} v_{ε}, & i n D ∖ B_{R}, \\ v_{ε, R}, & i n B_{R} ∖ B_{ε} . \end{array} \end{aligned}$

Proof.
We set
$\begin{aligned} υ := {\begin{cases} v_{ε}, & in D ∖ B_{R}, \\ v_{ε, R}, & in B_{R} ∖ B_{ε} . \end{cases} \end{aligned}$
(6)
Let φ be in $V_{ε}$ , then,
$\begin{aligned} {\tilde{a}}_{ε} (v, φ) & = α \int_{D ∖ B_{R}} \nabla v \cdot \nabla φ d x + \int_{D ∖ B_{R}} v φ d x + α \int_{B_{R} ∖ B_{ε}} \nabla v \cdot \nabla φ d x + \int_{B_{R} ∖ B_{ε}} v φ d x \end{aligned}$
Replacing v by its expression (6) and integrating by parts yields,
$\begin{aligned} {\tilde{a}}_{ε} (v, φ) & = \int_{D ∖ B_{R}} (- α Δ v_{ε} + v_{ε}) φ d x + \int_{B_{R} ∖ B_{ε}} (- α Δ v_{ε, R} + v_{ε, R}) φ d x \\ + α \int_{\partial B_{R}} \partial_{n_{int}} v_{ε} φ d σ + α \int_{\partial B_{R}} \partial_{n_{ext}} v_{ε, R} φ d σ \\ = \int_{D ∖ B_{ε}} h φ d x = {\tilde{l}}_{ε} (φ) . \end{aligned}$
By the uniqueness of the solution of Problem 3.1, we have $v = {\tilde{v}}_{ε}$ .

For the internal problem, we introduce the notation $v_{ε}^{h, ϕ}$ instead of $v_{ϵ, R}$ , the solution of the more general problem
Problem 3.3.
Find $v_{ε}^{h, ϕ}$ in ${v \in H^{1} (B_{R} ∖ B_{ε}) ∣ v = 0 ~on~ \partial B_{ε}}$ such that
$\begin{aligned} {\begin{array}{ll} - α Δ υ_{ε}^{h, ϕ} + υ_{ε}^{h, ϕ} = h, & i n B_{R} ∖ B_{ε}, \\ υ_{ε}^{h, ϕ} = 0, & o n \partial B_{ε}, \\ υ_{ε}^{h, ϕ} = ϕ, & o n \partial B_{R} . \end{array} \end{aligned}$
(7)

Therefore, $v_{ε, R} = v_{ε}^{h, ϕ}$ , when $ϕ = v_{ε}$ . We also notice that,
$\begin{aligned} v_{ε}^{h, ϕ} = v_{ε}^{h, 0} + v_{ε}^{0, ϕ} . \end{aligned}$
We remind the Dirichlet-to-Neumann operator $T_{ε} : H^{1 / 2} (\partial B_{R}) \to H^{- 1 / 2} (\partial B_{R})$ by
$\begin{aligned} T_{ε} (ϕ) := \nabla v_{ε}^{0, ϕ} \cdot n . \end{aligned}$
and we set
$\begin{aligned} h_{ε} := - \nabla v_{ε}^{h, 0} \cdot n \in H^{- 1 / 2} (\partial B_{R}) . \end{aligned}$
Hence, setting $V_{R} = H^{1} (D ∖ B_{R})$ , we can rewrite the external problem using this operator as following (we still denote by $v_{ε}$ the solution):
Problem 3.4.
Find $v_{ε}$ in $V_{R}$ such that
$\begin{aligned} {\begin{array}{ll} - α Δ v_{ε} + v_{ε} = h, & i n D ∖ B_{R}, \\ - \partial_{n} v_{ε} + T_{ε} v_{ε} = h_{ε}, & o n \partial B_{R}, \\ \partial_{n} v_{ε} = 0, & o n \partial D . \end{array} \end{aligned}$
(8)

For $ε \in [0, ζ]$ , $R > ζ$ , and v, φ in $V_{R} := H^{1} (D ∖ B_{R})$ , we define
$\begin{aligned} a_{ε} (v, φ) & := α \int_{D ∖ B_{R}} \nabla v \cdot \nabla φ d x + α \int_{\partial B_{R}} T_{ε} v φ d σ + \int_{D ∖ B_{R}} v φ d x, \\ l_{ε} (φ) & := \int_{D ∖ B_{R}} h φ d x + α \int_{\partial B_{R}} h_{ε} φ d σ . \end{aligned}$
So that the associated variational formulation reads: find $v \in V_{R}$ , such that
$\begin{aligned} a_{ε} (v, φ) = l_{ε} (φ), \forall φ \in V_{R} . \end{aligned}$
It is easily checked that $a_{ε}$ is symmetric and $l_{ε}$ is continuous.

We take as cost function, for $1 ⩽ p < + \infty$ ,
$\begin{aligned} {\tilde{J}}_{ε} (\tilde{υ}) := \int_{D ∖ B_{ε}} | h - \tilde{υ} |^{p} d x, \forall \tilde{υ} \in V_{ε}, \end{aligned}$

We define now the cost functional on $V_{R}$ as follows: for $v \in V_{R}$ , we set ${\tilde{v}}_{ε} \in V_{ε}$ the extension of v in $D ∖ B_{ε}$ such that,
${\tilde{v}}_{ε} |_{D ∖ B_{R}} = v$ ,

${\tilde{v}}_{ε} |_{B_{R} ∖ B_{ε}} = v_{ε}^{h, ϕ}$ , for $ϕ = v$ on $\partial B_{R}$ .
We notice that ${\tilde{v}}_{ε}$ do not satisfy Problem 3.1 except if v is the solution of Problem 3.4. Then, we may define the restriction of ${\tilde{J}}_{ε}$ to $V_{R}$ by:
$\begin{aligned} J_{ε} (v) := {\tilde{J}}_{ε} ({\tilde{v}}_{ε}), \forall v \in V_{R} . \end{aligned}$

3.1. The adjoint problem and related estimates

We state now the adjoint problem associated to Problem 3.4 when $ε = 0$ : we denote by $w_{0}$ the weak solution in $V_{R}$ of

\begin{aligned} a_{0} (w_{0}, φ) = - D J_{0} (v_{0}) φ, \forall φ \in V_{R}, \end{aligned}

where

v_{0}

is the solution of Problem 3.4. The adjoint state

w_{0}

is then the solution of

Problem 3.5.
Find $w_{0}$ in $V_{R}$ such that
$\begin{aligned} {\begin{array}{ll} - α Δ w_{0} + w_{0} = - v_{0} | v_{0} |^{p - 2}, & i n D ∖ B_{R}, \\ - \partial_{n} w_{0} + T_{0} w_{0} = h_{0}, & o n \partial B_{R}, \\ \partial_{n} w_{0} = 0, & o n \partial D . \end{array} \end{aligned}$
(9)

Note that, formally, when $p = 1$ , the right-hand side of (9) becomes $- sign (v_{0})$ .

We aim to find $δ a$ , $δ l$ and $δ J$ from the adjoint method, Theorem 3.1. Let $h \in L^{2} (D)$ and $ϕ \in H^{1 / 2} (\partial B_{R})$ .
3.2. Variations of the bilinear form

We start by giving an explicit formulation for both $v_{0}^{0, ϕ}$ and $v_{ε}^{0, ϕ}$ , which is analogous to [26], with the following proposition:

Proposition 3.2.
For ϕ in $H^{1 / 2} (\partial B_{R})$ , we have,
$\begin{aligned} v_{0}^{0, ϕ} (r, θ) = \sum_{n \in Z} \frac{I_{n} (α^{- 1 / 2} r)}{I_{n} (α^{- 1 / 2} R)} ϕ_{n} e^{i n θ}, \end{aligned}$
and,
$\begin{aligned} v_{ε}^{0, ϕ} (r, θ) = \sum_{n \in Z} \frac{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} r) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} r)}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} R)} ϕ_{n} e^{i n θ}, \end{aligned}$
where $(r, θ)$ are the polar coordinates in $R^{2}$ , $(ϕ_{n})_{n}$ are the Fourier coefficients of ϕ, $I_{n}$ and $K_{n}$ are the modified Bessel functions of the first and second kind respectively [ 1 , 23 ].
Proof.
Using polar coordinates in $R^{2}$ , we have,
$\begin{aligned} v_{ε}^{h, 0} (r, θ) = \sum_{n \in Z} c_{n} (r) e^{i n θ} and h (r, θ) = \sum_{n \in Z} h_{n} (r) e^{i n θ}, \end{aligned}$
where $c_{n}$ satisfies, for all n in $Z$ , and $0 < r ⩽ R$ ,
$\begin{aligned} - α r^{2} c_{n}^{^{″}} (r) - α r c_{n}^{^{'}} (r) + (r^{2} + α n^{2}) c_{n} (r) = 0. \end{aligned}$
We solve the equation, and we get,
$\begin{aligned} c_{0, n} (r) = A_{0, n} I_{n} (α^{- 1 / 2} r), \end{aligned}$
and
$\begin{aligned} c_{ε, n} (r) = A_{ε, n} I_{n} (α^{- 1 / 2} r) + B_{ε, n} K_{n} (α^{- 1 / 2} r) . \end{aligned}$
By using the boundaries conditions, we have the result.
Remark.
According to [1,23], we have for all n in $Z$ , $I_{- n} = I_{n}$ and $K_{- n} = K_{n}$ .

Next, we state the variation of the solution with respect to hole’s radius:
Proposition 3.3.
For ϕ in $H^{1 / 2} (\partial B_{R})$ , we set,
$\begin{aligned} δ v^{0, ϕ} (r) := - ϕ_{0} \frac{I_{0} (α^{- 1 / 2} R) K_{0} (α^{- 1 / 2} r) - K_{0} (α^{- 1 / 2} R) I_{0} (α^{- 1 / 2} r)}{I_{0} (α^{- 1 / 2} R)^{2}} . \end{aligned}$
Then, for $ε$ sufficiently small, we have the following asymptotic estimation,
$\begin{aligned} v_{ε}^{0, ϕ} (r, θ) - v_{0}^{0, ϕ} (r, θ) - \frac{- 1}{\ln ε} δ v^{0, ϕ} (r) = o (\frac{- 1}{\ln ε}) . \end{aligned}$

Proof.
We have,
$\begin{aligned} v_{ε}^{0, ϕ} (r, θ) - v_{0}^{0, ϕ} (r, θ) & = (\frac{I_{0} (α^{- 1 / 2} ε) K_{0} (α^{- 1 / 2} r) - K_{0} (α^{- 1 / 2} ε) I_{0} (α^{- 1 / 2} r)}{I_{0} (α^{- 1 / 2} ε) K_{0} (α^{- 1 / 2} R) - K_{0} (α^{- 1 / 2} ε) I_{0} (α^{- 1 / 2} R)} - \frac{I_{0} (α^{- 1 / 2} r)}{I_{0} (α^{- 1 / 2} R)}) ϕ_{0} + R_{ε} (r, θ), \end{aligned}$
where,
$\begin{aligned} R_{ε} (r, θ) := \sum_{n \in Z^{}} (\frac{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} r) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} r)}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} R)} - \frac{I_{n} (α^{- 1 / 2} r)}{I_{n} (α^{- 1 / 2} R)}) ϕ_{n} e^{i n θ} . \end{aligned}$
Then,
$\begin{aligned} v_{ε}^{0, ϕ} (r, θ) - v_{0}^{0, ϕ} (r, θ) = - δ v^{0, ϕ} (r) \frac{I_{0} (α^{- 1 / 2} R) I_{0} (α^{- 1 / 2} ε)}{I_{0} (α^{- 1 / 2} ε) K_{0} (α^{- 1 / 2} R) - K_{0} (α^{- 1 / 2} ε) I_{0} (α^{- 1 / 2} R)} + R_{ε} (r, θ) . \end{aligned}$
We use that, [1],
$\begin{aligned} K_{0} (x) = - (γ - \ln 2 + \ln x) I_{0} (x) + x r_{1} (x), \end{aligned}$
and get,
$\begin{aligned} v_{ε}^{0, ϕ} (r, θ) - v_{0}^{0, ϕ} (r, θ) = - δ v^{0, ϕ} (r) (M + I_{0} (α^{- 1 / 2} R) \ln ε + ε r_{2} (ε))^{- 1} + R_{ε} (r, θ), \end{aligned}$
where,
$\begin{aligned} M := \frac{K_{0} (α^{- 1 / 2} R) + I_{0} (α^{- 1 / 2} R) (γ - \ln 2 + \ln (α^{- 1 / 2}))}{I_{0} (α^{- 1 / 2} R)}, \end{aligned}$
is a constant independent of $ε$ . Finally,
$\begin{aligned} v_{ε}^{0, ϕ} (r, θ) - v_{0}^{0, ϕ} (r, θ) = - δ v^{0, ϕ} (r) \frac{- 1}{\ln ε} {(1 + \frac{M}{\ln ε} + ε r_{3} (ε))}^{- 1} + R_{ε} (r, θ) . \end{aligned}$

It remains to show that $R_{ε} (r, θ) = o (\frac{- 1}{\ln ε})$ . We have,
$\begin{aligned} R_{ε} (r, θ) := \frac{- 1}{\ln ε} Θ_{ε} (r, θ), \end{aligned}$
where,
$\begin{aligned} Θ_{ε} (r, θ) := \sum_{n \in Z^{}} \ln ε \frac{K_{n} (α^{- 1 / 2} R) I_{n} (α^{- 1 / 2} r) - I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} r)}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} ε)} \frac{I_{n} (α^{- 1 / 2} ε)}{I_{n} (α^{- 1 / 2} R)} ϕ_{n} e^{i n θ} . \end{aligned}$
Moreover,
$\begin{aligned} | Θ_{ε} (r, θ) | ⩽ \sum_{n \in Z^{}} | \ln ε | | \frac{I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} r) - K_{n} (α^{- 1 / 2} R) I_{n} (α^{- 1 / 2} r)}{I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} ε) - I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R)} | | \frac{I_{n} (α^{- 1 / 2} ε)}{I_{n} (α^{- 1 / 2} R)} | | ϕ_{n} | . \end{aligned}$
Since $I_{n}$ is an increasing function and $K_{n}$ is a decreasing function, we have, for $r \in [ε, R]$ ,
$\begin{aligned} | & I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} r) - K_{n} (α^{- 1 / 2} R) I_{n} (α^{- 1 / 2} r) | \\ ⩽ | I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} ε) - K_{n} (α^{- 1 / 2} R) I_{n} (α^{- 1 / 2} ε) | . \end{aligned}$
Thus,
$\begin{aligned} | Θ_{ε} (r, θ) | ⩽ \sum_{n \in Z^{}} | \ln ε | | \frac{I_{n} (α^{- 1 / 2} ε)}{I_{n} (α^{- 1 / 2} R)} | | ϕ_{n} | . \end{aligned}$
Finally, using that
$\begin{aligned} I_{n} (α^{- 1 / 2} ε) = \frac{1}{n!} {(\frac{α^{- 1 / 2}}{2})}^{n} ε^{n} + ε^{n} r (ε), \end{aligned}$
we have that $Θ_{ε}$ tends to 0 when $ε$ goes to 0.

By noticing that [1],
$\begin{aligned} I_{0} (x) K_{0}^{'} (x) - K_{0} (x) I_{0}^{'} (x) = - I_{0} (x) K_{1} (x) - K_{0} (x) I_{1} (x) = - W (K_{0} (x), I_{0} (x)) = - \frac{1}{x}, \end{aligned}$
(10)
with W the Wronskian and by using the previous result, we have the following:
Proposition 3.4.
For ϕ in $H^{1 / 2} (\partial B_{R})$ , we define,
$\begin{aligned} δ T (ϕ) := ϕ_{0} \frac{1}{R I_{0} (α^{- 1 / 2} R)^{2}} . \end{aligned}$
Then, for $ε$ sufficiently small, we have the following asymptotic estimation,
$\begin{aligned} {‖ T_{ε} - T_{0} - \frac{- 1}{\ln ε} δ T ‖}_{L (H^{1 / 2} (\partial B_{R}), H^{- 1 / 2} (\partial B_{R}))} = o (\frac{- 1}{\ln ε}) . \end{aligned}$

Finally, we can derive from the previous proposition the variations of the bilinear form: Proposition 3.5.
For ϕ in $H^{1 / 2} (\partial B_{R})$ , we define,
$\begin{aligned} δ a (v, w) := α \frac{v^{mean}}{I_{0} (α^{- 1 / 2} R)} \frac{w^{mean}}{I_{0} (α^{- 1 / 2} R)}, \end{aligned}$
where $v^{mean}$ and $w^{mean}$ denote the mean value of v and w on $\partial B_{R}$ . Then, for $ε$ sufficiently small, we have the following asymptotic estimation,
$\begin{aligned} | a_{ε} (v, w) - a_{0} (v, w) - \frac{- 2 π}{\ln ε} δ a (v, w) | = o (\frac{- 1}{\ln ε}) ‖ v ‖_{V_{R}} ‖ w ‖_{V_{R}}, \forall v, w \in V_{R} . \end{aligned}$

3.3. Variations of the linear form

We give an explicit formulation for both $v_{0}^{h, 0}$ and $v_{ε}^{h, 0}$ with the following proposition:

Proposition 3.6.
For h in $L^{2} (B_{R})$ , we have,
$\begin{aligned} v_{0}^{h, 0} (r, θ) = \sum_{n \in Z} ((A_{0, n} + A_{n}^{p} (r)) I_{n} (α^{- 1 / 2} r) + B_{n}^{p} (r) K_{n} (α^{- 1 / 2} r)) e^{i n θ}, \end{aligned}$
with,
$\begin{aligned} A_{n}^{p} (r) & := - α^{- 1} \int_{0}^{r} s K_{n} (α^{- 1 / 2} s) h_{n} (s) d s, \\ B_{n}^{p} (r) & := α^{- 1} \int_{0}^{r} s I_{n} (α^{- 1 / 2} s) h_{n} (s) d s, \\ A_{0, n} & := - \frac{A_{n}^{p} (R) I_{n} (α^{- 1 / 2} R) + B_{n}^{p} (R) K_{n} (α^{- 1 / 2} R)}{I_{n} (α^{- 1 / 2} R)} . \end{aligned}$

Moreover,
$\begin{aligned} v_{ε}^{h, 0} (r, θ) = \sum_{n \in Z} ((A_{ε, n} + A_{n}^{p} (r)) I_{n} (α^{- 1 / 2} r) + (B_{ε, n} + B_{n}^{p} (r)) K_{n} (α^{- 1 / 2} r)) e^{i n θ}, \end{aligned}$
with,
$\begin{aligned} A_{ε, n} & = \frac{K_{n} (α^{- 1 / 2} ε) (A_{n}^{p} (R) I_{n} (α^{- 1 / 2} R) + B_{n}^{p} (R) K_{n} (α^{- 1 / 2} R))}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} R)}, \\ B_{ε, n} & = - \frac{I_{n} (α^{- 1 / 2} ε) (A_{n}^{p} (R) I_{n} (α^{- 1 / 2} R) + B_{n}^{p} (R) K_{n} (α^{- 1 / 2} R))}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} R)} . \end{aligned}$

Proof.
Using polar coordinates in $R^{2}$ , we have,
$\begin{aligned} v_{ε}^{h, 0} (r, θ) = \sum_{n \in Z} c_{n} (r) e^{i n θ} and h (r, θ) = \sum_{n \in Z} h_{n} (r) e^{i n θ}, \end{aligned}$
where $c_{n}$ satisfies, for all n in $Z$ and $0 < r ⩽ R$ ,
$\begin{aligned} - α r^{2} c_{n}^{″} (r) - α r c_{n}^{'} (r) + (r^{2} + α n^{2}) c_{n} (r) = r^{2} h_{n} (r) . \end{aligned}$
(11)
Firstly, we solve the homogeneous equation, and we get,
$\begin{aligned} c_{0, n}^{h} (r) = A_{0, n} I_{n} (α^{- 1 / 2} r), \end{aligned}$
and
$\begin{aligned} c_{ε, n}^{h} (r) = A_{ε, n} I_{n} (α^{- 1 / 2} r) + B_{ε, n} K_{n} (α^{- 1 / 2} r) . \end{aligned}$
Secondly, we use the variation of parameters method to get the particular solution,
$\begin{aligned} c_{n}^{p} (r) = A_{n}^{p} (r) I_{n} (α^{- 1 / 2} r) + B_{n}^{p} (r) K_{n} (α^{- 1 / 2} r) . \end{aligned}$
(12)
By replacing (12) into (11), and by supposing that,
$\begin{aligned} (A_{n}^{p})^{^{'}} (r) I_{n} (α^{- 1 / 2} r) + (B_{n}^{p})^{^{'}} (r) K_{n} (α^{- 1 / 2} r) = 0, \end{aligned}$
we get,
$\begin{aligned} (A_{n}^{p})^{^{'}} (r) I_{n}^{^{'}} (α^{- 1 / 2} r) + (B_{n}^{p})^{^{'}} (r) K_{n}^{^{'}} (α^{- 1 / 2} r) = - α^{- 1 / 2} h_{n} (r) . \end{aligned}$
Solving the last two equations, we get the value of $A_{n}^{p} (r)$ and $B_{n}^{p}) (r)$ as stated in the theorem. Finally, by the superposition principle and by using the boundaries conditions, we have the result.
Remark.
Since h is in $L^{2} (D)$ , using the Parseval’s equality we have that the Fourier coefficients $h_{n}$ are in $L^{2} (D)$ as well. Similarly for the Fourier coefficients $c_{n}$ . As a result, we have that the integral in the $A_{n}^{p}$ is convergent.
Remark.
We have,
$\begin{aligned} v_{0}^{h, 0} (x_{0}) = \sum_{n \in Z} A_{0, n} I_{n} (0) e^{i n θ}, \forall θ \in [0, 2 π [, \end{aligned}$
and in particular, for $θ = 0$ , we have
$\begin{aligned} v_{0}^{h, 0} (x_{0}) = \sum_{n \in Z} A_{0, n} I_{n} (0) . \end{aligned}$
Moreover, [1,23], $I_{n} (0)$ vanishes for $n \in N^{}$ and $I_{0} (0) = 1$ . Then,
$\begin{aligned} A_{0, 0} = v_{0}^{h, 0} (x_{0}) . \end{aligned}$
(13)

Now, we can state the variation of the solution with respect to hole’s radius:
Proposition 3.7.
For h in* $L^{2} (B_{R})$ , we set,
$\begin{aligned} δ v^{h, 0} (r) := - A_{0, 0} \frac{I_{0} (α^{- 1 / 2} R) K_{0} (α^{- 1 / 2} r) - K_{0} (α^{- 1 / 2} R) I_{0} (α^{- 1 / 2} r)}{I_{0} (α^{- 1 / 2} R)}, \end{aligned}$
where $A_{0, 0}$ is defined in Proposition3.6. Then, for $ε$ sufficiently small, we have the following asymptotic estimation,
$\begin{aligned} v_{ε}^{h, 0} (r, θ) - v_{0}^{h, 0} (r, θ) - \frac{- 1}{\ln ε} δ v^{h, 0} (r) = o (\frac{- 1}{\ln ε}) . \end{aligned}$

Proof.
We have,
$\begin{aligned} v_{ε}^{h, 0} (r, θ) - v_{0}^{h, 0} (r, θ) = (A_{ε, 0} - A_{0, 0}) I_{0} (α^{- 1 / 2} r) + B_{ε, 0} K_{0} (α^{- 1 / 2} r) + R_{ε} (r, θ), \end{aligned}$
and,
$\begin{aligned} R_{ε} (r, θ) := \sum_{n \in Z^{}} ((A_{ε, n} - A_{0, n}) I_{n} (α^{- 1 / 2} r) + B_{ε, n} K_{n} (α^{- 1 / 2} r)) e^{i n θ} . \end{aligned}$
We have, for all n in $Z$ ,
$\begin{aligned} A_{ε, n} - A_{0, n} = - A_{0, n} \frac{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R)}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} R)}, \end{aligned}$
and,
$\begin{aligned} B_{ε, n} = A_{0, n} \frac{I_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} R)}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} R)} . \end{aligned}$
Thus,
$\begin{aligned} v_{ε}^{h, 0} (r, θ) - v_{0}^{h, 0} (r, θ) = - δ v^{h, 0} (r) \frac{I_{0} (α^{- 1 / 2} R) I_{0} (α^{- 1 / 2} ε)}{I_{0} (α^{- 1 / 2} ε) K_{0} (α^{- 1 / 2} R) - K_{0} (α^{- 1 / 2} ε) I_{0} (α^{- 1 / 2} R)} + R_{ε} (r), \end{aligned}$
and,
$\begin{aligned} R_{ε} (r, θ) := - \sum_{n \in Z^{}} A_{0, n} \frac{I_{n} (α^{- 1 / 2} ε) (I_{n} (α^{- 1 / 2} r) K_{n} (α^{- 1 / 2} R) - K_{n} (α^{- 1 / 2} r) I_{n} (α^{- 1 / 2} R))}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - K_{n} (α^{- 1 / 2} ε) I_{n} (α^{- 1 / 2} R)} e^{i n θ} . \end{aligned}$
We use that, [1],
$\begin{aligned} K_{0} (x) = - (γ - \ln 2 + \ln x) I_{0} (x) + x r_{1} (x), \end{aligned}$
and get,
$\begin{aligned} v_{ε}^{h, 0} (r, θ) - v_{0}^{h, 0} (r, θ) = - δ v^{h, 0} (r) (M + \ln ε + ε r_{2} (ε))^{- 1} + R_{ε} (r) . \end{aligned}$
where,
$\begin{aligned} M := \frac{K_{0} (α^{- 1 / 2} R) + I_{0} (α^{- 1 / 2} R) (γ - \ln 2 + \ln (α^{- 1 / 2}))}{I_{0} (α^{- 1 / 2} R)}, \end{aligned}$
is a constant independent of $ε$ . It remains to show that $R_{ε} (r, θ) = o (\frac{- 1}{\ln ε})$ . We have,
$\begin{aligned} R_{ε} (r, θ) := \frac{- 1}{\ln ε} Θ_{ε} (r, θ), \end{aligned}$
where,
$\begin{aligned} Θ_{ε} (r, θ) := \sum_{n \in Z^{}} \ln ε \frac{K_{n} (α^{- 1 / 2} R) I_{n} (α^{- 1 / 2} r) - I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} r)}{I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R) - I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} ε)} I_{n} (α^{- 1 / 2} ε) A_{0, n} e^{i n θ} . \end{aligned}$
Moreover,
$\begin{aligned} | Θ_{ε} (r, θ) | ⩽ \sum_{n \in Z^{}} | \ln ε | | \frac{I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} r) - K_{n} (α^{- 1 / 2} R) I_{n} (α^{- 1 / 2} r)}{I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} ε) - I_{n} (α^{- 1 / 2} ε) K_{n} (α^{- 1 / 2} R)} | | I_{n} (α^{- 1 / 2} ε) | | A_{0, n} | . \end{aligned}$
Since $I_{n}$ is an increasing function and $K_{n}$ is a decreasing function, we have, for $r \in [ε, R]$ ,
$\begin{aligned} | I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} r) - K_{n} (α^{- 1 / 2} R) I_{n} (α^{- 1 / 2} r) | \\ ⩽ | I_{n} (α^{- 1 / 2} R) K_{n} (α^{- 1 / 2} ε) - K_{n} (α^{- 1 / 2} R) I_{n} (α^{- 1 / 2} ε) | . \end{aligned}$
Thus,
$\begin{aligned} | Θ_{ε} (r, θ) | ⩽ \sum_{n \in Z^{}} | \ln ε | | I_{n} (α^{- 1 / 2} ε) | | A_{0, n} | . \end{aligned}$
Finally, using that,
$\begin{aligned} I_{n} (α^{- 1 / 2} ε) = \frac{1}{n!} {(\frac{α^{- 1 / 2}}{2})}^{n} ε^{n} + ε^{n} r (ε), \end{aligned}$
we have that $Θ_{ε}$ tends to 0 when $ε$ goes to 0.

Using the previous result and the property of the Wronskian [1], we have the following:
Proposition 3.8.
For h in* $L^{2} (B_{R})$ , we define,
$\begin{aligned} δ h (h) := - A_{0, 0} \frac{1}{R I_{0} (α^{- 1 / 2} R)} . \end{aligned}$
where $A_{0, 0}$ is defined in Proposition 3.6. Then, for $ε$ sufficiently small, we have the following asymptotic estimation,
$\begin{aligned} {‖ h_{ε} - h_{0} - \frac{- 1}{\ln ε} δ h ‖}_{- 1 / 2, \partial B_{R}} = o (\frac{- 1}{\ln ε}) . \end{aligned}$

Finally, we can derive from the previous proposition the variations of the linear form: Proposition 3.9.
For h in $L^{2} (B_{R})$ , we define,
$\begin{aligned} δ l (w) := - α A_{0, 0} \frac{w^{mean}}{I_{0} (α^{- 1 / 2} R)} . \end{aligned}$
Then, for $ε$ sufficiently small, we have the following asymptotic estimation,
$\begin{aligned} | l_{ε} (w) - l_{0} (w) - \frac{- 2 π}{\ln ε} δ l (w) | = o (\frac{- 1}{\ln ε}) ‖ w ‖_{V_{R}}, \forall w \in V_{R} . \end{aligned}$

4. Computation of the topological derivative

We now gather the previous section results to derive the topological derivative. We consider the adjoint problem of Problem 3.2:

Problem 4.1.
Find ${\tilde{w}}_{0}$ in $H^{1} (D)$ such that
$\begin{aligned} {\begin{cases} - α Δ {\tilde{w}}_{0} + {\tilde{w}}_{0} = - {\tilde{v}}_{0} | {\tilde{v}}_{0} |^{p - 2}, & i n D, \\ \partial_{n} {\tilde{w}}_{0} = 0, & o n \partial D . \end{cases} \end{aligned}$
(14)

Then, we have the following proposition,
Proposition 4.1.
$w_{0}$ , solution of Problem 3.5, is the restriction of ${\tilde{w}}_{0}$ to $D ∖ B_{R}$ .
Proof.
We set $w_{R} := {\tilde{w}}_{0} |_{D ∖ B_{R}}$ . We have to show that $w_{R} = w_{0}$ i.e. $a_{0} (w_{R}, φ_{R}) = - D J_{0} (v_{0}) φ_{R}$ , $\forall φ_{R} \in V_{R}$ . Let $φ_{R} \in V_{R}$ . We denote $\tilde{φ} \in V_{0}$ the extension of $φ_{R}$ to $V_{0}$ such that $- α Δ \tilde{φ} + \tilde{φ} = 0$ in $B_{R}$ . Thus,
$\begin{aligned} a_{0} (w_{R}, φ_{R}) & = α \int_{D ∖ B_{R}} \nabla w_{R} \cdot \nabla φ_{R} d x + α \int_{\partial B_{R}} T_{0} w_{R} φ_{R} d σ + \int_{D ∖ B_{R}} w_{R} φ_{R} d x \\ = α \int_{D ∖ B_{R}} \nabla w_{R} \cdot \nabla φ_{R} d x + α \int_{\partial B_{R}} T_{0} w_{R} φ_{R} d σ + \int_{D ∖ B_{R}} w_{R} φ_{R} d x \\ + \int_{B_{R}} \underset{= 0}{\underset{⏟}{(- α Δ \tilde{φ} + \tilde{φ})}} {\tilde{w}}_{0} d x . \end{aligned}$
And after integrating by parts,
$\begin{aligned} a_{0} (w_{R}, φ_{R}) & = α \int_{D ∖ B_{R}} \nabla w_{R} \cdot \nabla φ_{R} d x + \int_{D ∖ B_{R}} w_{R} φ_{R} d x + α \int_{B_{R}} \nabla \tilde{φ} \cdot \nabla {\tilde{w}}_{0} d x + \int_{B_{R}} \tilde{φ} {\tilde{w}}_{0} d x \\ = α \int_{D} \nabla {\tilde{w}}_{0} \cdot \nabla \tilde{φ} d x + \int_{D} {\tilde{w}}_{0} \tilde{φ} d x \\ = {\tilde{a}}_{0} ({\tilde{w}}_{0}, \tilde{φ}) = - D {\tilde{J}}_{0} (v_{D}) \tilde{φ} . \end{aligned}$
Moreover, by definition ${\tilde{J}}_{0} (v_{D}) = J_{0} (v_{R})$ , thus,
$\begin{aligned} D {\tilde{J}}_{0} (v_{D}) \tilde{φ} = D J_{0} (v_{0}) φ_{R} . \end{aligned}$
By uniqueness of the solution, $w_{R} = w_{0}$ .

Using that, if u and v are solutions of the linear diffusion equation on $B_{R}$ ,
$\begin{aligned} v (x_{0}) = \frac{v^{mean}}{I_{0} (α^{- 1 / 2} R)} and w (x_{0}) = \frac{w^{mean}}{I_{0} (α^{- 1 / 2} R)}, \end{aligned}$
it follows the topological gradient based on the adjoint method given by:
Proposition 4.2.
For $ε$ small enough, we have,
$\begin{aligned} j (K_{ε}) - j (K) = 2 α {\tilde{v}}_{0} (x_{0}) {\tilde{w}}_{0} (x_{0}) \frac{- 2 π}{\ln ε} + o (\frac{- 1}{\ln ε}), \end{aligned}$
with ${\tilde{v}}_{0}$ solution of Problem 3.2 and ${\tilde{w}}_{0}$ solution of Problem 4.1.

We notice that with this expansion, we get the main theoretical result of the paper which might be summarized as follows: to minimize the $L^{p}$ -error between an image and its reconstruction from linear diffusion inpainting, we have to keep in the mask the pixels $x_{0}$ which minimize the product ${\tilde{v}}_{0} (x_{0}) {\tilde{w}}_{0} (x_{0})$ .
5. Numerical results

In this section, we present numerical results when the cost functional is the $L^{1}$ -error and the $L^{2}$ -error, respectively, as they are the most representative for noise in practice. In fact, we take these specific values of p, depending on the nature of the noise considered: $p = 2$ is well suited for Gaussian noise while $p = 1$ is better for impulse noise such as salt and pepper noise.

Let us denote by f the original image, $f_{δ}$ the noisy one, K the inpainting mask, and u the reconstructed image. Proposition 4.2 gives a hard-threshold criterion for the selection of K, namely, we must keep the pixels $x \in D$ which minimize the quantity ${\tilde{v}}_{0} (x) {\tilde{w}}_{0} (x)$ (or equivalently, which maximize $- {\tilde{v}}_{0} (x) {\tilde{w}}_{0} (x)$ ). In fact, the adjoint state is obtained by solving a linear PDE which is, by the elliptic regularity, $H^{2}$ (the right-hand side is $L^{2}$ even when $p = 1$ with the regularization adopted in the article). The distribution of this function measures the influence of a considered single pixel and its neighborhood in the cost variations. This is referred to as soft-thresholding (in contrast to hard-thresholding with pixel-wise asymptotics) and is implemented using the Floyd-Steinberg halftoning technique [29]. Digital halftoning is a technique used to simulate continuous-tone images using a limited palette of colors or shades, often just black and white. The compression step states as follows:

solve Problem 3.2 with $h = α Δ f_{δ}$ , and save the result in ${\tilde{v}}_{0}$ ,

solve Problem 4.1 and save the result in ${\tilde{w}}_{0}$ ,

apply hard/soft-thresholding to $- {\tilde{v}}_{0} {\tilde{w}}_{0}$ and store the selected pixels in K.

We emphasize that the inpainting mask is built from

f_{δ}

, which is available in D during the mask selection step, while it is only available in K for the reconstruction.

The reconstruction step is performed by solving the heat equation with the semi-implicit discrete scheme of time step $d t > 0$ and with initial condition $u_{0}$ . In fact, Problem 2.1 is formally the first iteration of the semi-implicit discrete scheme of the heat equation, and directly solving Problem 2.1 will result to a small diffusion of the data, leading to a poor reconstruction’s quality. We choose the time step $d t = 10$ and the initial condition $u_{0} = 1_{K} f_{δ}$ for boundary conditions compatibility. The diffusion process is stopped as soon as the error $e_{n}$ between the reconstruction $u_{n}$ and the noisy data $f_{δ}$ at the iteration $n \in N^{*}$ starts to increase i.e. when $e_{n} ⩾ e_{n - 1}$ where $e_{n} := ‖ u_{n} - f_{δ} ‖_{L^{p} (D)}$ .

We denote by Lp-ADJ-T the algorithm using the adjoint method by selecting the pixels using hard-thresholding, and by Lp-ADJ-H the algorithm combining with the Floyd-Steinberg halftoning. For comparison purposes, we consider the methods H1-T and H1-H, which correspond to mask selection based on the asymptotic expansion provided in [6]. In these methods, the authors use a hard/soft-thresholding of $| Δ f_{δ} |$ as pixel selection criterion, combined with homogeneous diffusion inpainting for the reconstruction. Note that the homogeneous diffusion equation can be viewed as a limiting case of the heat equation as t approaches infinity. In their paper, Belhachmi et al stated that they smoothed the data prior to computing its Laplacian. However, since we aim to demonstrate the effectiveness of our methods on noisy data, we believe that pre-processing the data introduces an unfair advantage. For this reason, we do not smooth $f_{δ}$ when calculating its Laplacian, although we do apply smoothing for the H1- methods as described by the authors.

To solve the equations, we discretize the problems using the finite difference method where every degree of freedom correspond to a pixel of the image, and we manually select the parameter α to achieve the best reconstruction. The implementation of the algorithms is done in Python and can be found in [16].

5.1. Impulse noise

We begin the experiments by evaluating the robustness of our methods against impulse noise. Impulse noise, such as salt-and-pepper noise, is often addressed by minimizing the $L^{1}$ -error, as discussed in [21,22]. Here, we focus on the case $p = 1$ , and our proposed methods are refereed by L1-ADJ-. Figure 2 shows the input images in the following order, from left to right: the noiseless image, the image with salt noise, the one with pepper noise, and the image with a combination of both salt and pepper noise.

Fig. 2.

Input images.

We give in Table 1 (and in Appendix A) the $L^{1}$ -error for the H1- and L1-ADJ- methods for several amounts of salt and/or pepper noise, and hand-picked regularization parameter α ( $‖ \cdot ‖_{1}$ denotes the $L^{1}$ norm). The first notable observation is that L1-ADJ-H gives the lowest $L^{1}$ -error with and without presence of noise. Additionally, consistent with findings from [6,7], the soft-thresholding rule consistently results in a more effective inpainting mask for image reconstruction.

Table 1

$L^{1}$ -error between the original image f and the reconstruction u (built from $f_{δ}$ ) with $10 %$ of total pixels saved.

Noise		L1-ADJ-T		L1-ADJ-H		H1-T	H1-H
Salt	Pepper	α	$‖ f - u ‖_{1}$	α	$‖ f - u ‖_{1}$	$‖ f - u ‖_{1}$	$‖ f - u ‖_{1}$
0	0	0.01	4248.14	2.27	904.21	4648.99	961.01
0.02	0	0.41	3203.43	0.66	1194.91	6531.56	3875.63
0	0.02	0.36	2505.63	0.71	1204.80	5461.22	4019.73
0.01	0.01	0.56	2252.71	0.76	1321.55	4573.62	3097.18
0.04	0	0.36	4132.62	0.56	1453.11	12381.64	9611.86
0	0.04	2.07	2938.43	0.56	1454.21	15171.93	11195.98
0.02	0.02	0.46	2912.86	0.61	2098.27	6654.28	5901.86
0.1	0	0.26	8414.41	0.51	1845.45	28163.84	22580.17
0	0.1	2.42	4320.62	0.51	1840.34	34035.09	27595.06
0.05	0.05	0.36	6191.32	0.51	4623.07	17716.39	15806.10

Figure 3 (and Appendix A) shows the resulting masks and reconstructions corresponding to the data in the tables. We note that most of the corrupted pixels are not selected in K for the L1-ADJ- methods, while they are in the case of H1- ones. It is explained as follows: the asymptotic provided for the H1- methods leads to a pixel selection criterion which is proportional to $| Δ f_{δ} |^{2}$ . In the presence of impulse noise, the Laplacian is high at the edges, with even greater intensities at noisy pixels. Thus, noisy pixels are selected as a priority on pixels close to edges. On the other hand, the L1-ADJ- methods keep pixel in the inpainting mask that are located near the edges and around noisy pixels. In fact, for well-chosen α, if the Laplacian is too large at an isolated pixel, we are in presence of the noise and if it is too small, we are in homogeneous parts, and in both case, the L1-ADJ- methods do not select the pixel.

Fig. 3.

Masks and reconstructions from image with $2 %$ of salt and pepper noise and with $10 %$ of total pixels saved.

To provide further clarification on this last statement, we present the following case study: a binary image with two corrupted pixels (see Fig. 4 (a)). On Fig. 4 (a), (b), we plot the value of the criterion evaluated at several positions (noise, edge and homogeneous) with respect to the value of the regularization parameter α. Formally, the criterion reads, for a well-chosen $γ > 0$ , depending on α,

\begin{aligned} α (K_{γ} ⋆ Δ f_{δ}) (K_{γ} ⋆ sign (Δ f_{δ})), \end{aligned}

where

⋆

denote the convolution product and

K_{γ}

is the heat kernel, fundamental solution of the linear diffusion equation. When α is small, the criterion is close to

α sign (Δ f_{δ}) Δ f_{δ} = α | Δ f_{δ} |

which behaves similarly to the H1- method, and, in particular, the noisy pixel is kept. For large α, the criterion in the homogeneous part is dominant, and the edges of the image are lost. There is a window for the value of α in which the edge pixels are sectioned. We give in Fig. 5 illustrations for under/over-regularization.

Fig. 4.

Influence of α on the criterion evaluated at a noisy pixel, an edge pixel and a pixel in the homogeneous part.

Fig. 5.

Example of under-regularization on the left and of over-regularization on the right. For both case, L1-ADJ-H is used on the salt and pepper corrupted input (Fig. 2 (d)).

5.2. Gaussian noise

Now, we consider images with Gaussian noise and we take $p = 2$ . Figure 6 presents the input images without noise on the left, followed by the noisy images.

Fig. 6.

Input images $f_{δ}$ with Gaussian noise of deviation σ.

Table 2 (and Appendix B) gives the $L^{2}$ -error for the L2-ADJ- methods and the H1-T methods with respect to the deviation $σ > 0$ of Gaussian noise and the value of α (we write $‖ \cdot ‖_{2}$ the $L^{2}$ norm). For a reasonable level of noise, the L2-ADJ-H gives the lowest $L^{2}$ -error. Moreover, and as for the L1-ADJ- methods, using digital halftoning enhance the reconstruction by giving a better distribution of the pixel in the mask.

Table 2

$L^{2}$ -error between the original image f and the reconstruction u (built from $f_{δ}$ ) with $10 %$ of total pixels saved.

Noise	L2-ADJ-T		L2-ADJ-H		H1-T	H1-H
σ	α	$‖ f - u ‖_{2}$	α	$‖ f - u ‖_{2}$	$‖ f - u ‖_{2}$	$‖ f - u ‖_{2}$
0	0.01	23.08	0.01	9.70	25.57	4.99
0.03	0.71	9.38	0.96	7.91	9.23	8.57
0.05	0.86	13.25	0.76	12.39	13.78	12.55
0.1	0.71	26.95	0.66	24.47	27.15	22.95
0.2	0.01	56.08	2.27	46.46	61.94	42.62

In Fig. 7 (and in Appendix B), we present the resulting masks and reconstruction for various noise level. In the presence of Gaussian noise, the Laplacian amplifies the noise, and we cannot distinguish the edges from the noise if the noise level is sufficiently large. This explains why the H1- methods result in poor reconstructions. Conversely, the edges of the image are clearly discernible in the L2-ADJ- masks, and the reconstructed image appears to contain less noise compared to the input.

Fig. 7.

Masks and reconstructions from image with Gaussian noise of deviation $σ = 0.1$ and with $10 %$ of total pixels saved.

For small α, the criterion $- {\tilde{v}}_{0} {\tilde{w}}_{0}$ is formally close to $α | Δ f |^{2}$ which is similar to the result found in [7], but the fact that the adjoint state ${\tilde{w}}_{0}$ and primal variable ${\tilde{v}}_{0}$ are computed by solving linear PDEs improves distribution of the topological derivative. When α is large, we lose the information given by the edges, as showed in Fig. 8.

Fig. 8.

Example of under-regularization on the left and of over-regularization on the right. For both case, L2-ADJ-H is used on the salt and pepper corrupted input (Fig. 6 (c)).

6. Conclusion and perspectives

In this article, we formulated a (PDE-)compression problem based on the linear diffusion inpainting as a shape optimization one, and we estimated the topological expansion for the optimality condition by the adjoint method. More precisely, and following the approaches of [2,14], we computed the asymptotic development with respect to insertion of a circular hole for variations of the $L^{p}$ -error with general exponents $p ⩾ 1$ , which leads to an analytic criterion to select relevant pixels of the mask using a hard-thresholding rule. This criterion is determined by solving the linear diffusion equation and its adjoint equation which make the method easy to implement with a reasonable cost as both equations are linear. Finally, we presented numerical results in the two following settings: when the data is corrupted by salt and pepper noise ( $p = 1$ ) and when the data is affected by Gaussian noise ( $p = 2$ ). In these experiments, we relaxed the hard-thresholding from the topological gradient by introducing a soft-thresholding, enforced with the Floyd-Steinberg digital halftoning [29]. This new approach delivers better results than hard-thresholding by giving a better distribution of the selected pixels. Then, we compared our methods with those presented in [6], which is also a PDE-based compression problem, but rely on homogeneous diffusion inpainting and use the $H^{1}$ semi-norm as the cost function. Our mask selection methods outperform the other approaches when compressing images affected by Gaussian or impulse noise. Specifically, our methods effectively denoise the data by excluding corrupted pixels from the mask, resulting in a cleaner reconstruction. In contrast, the homogeneous diffusion inpainting methods include all noisy pixels, as the Laplacian of the noisy image is high around the corrupted pixels.

For future research, we aim to explore additional inpainting operators and alternative cost functionals to address various types of noise or combinations of different noise patterns. In the numerical section, the selection of α proves to be critical: if it is too low, there is no denoising, preserving the corrupted pixels; if it is too high, we end up losing the edges. However, determining the appropriate value of α remains unclear, which opens up potential research direction. Another goal would be to extend the method to colored images. Finally, further numerical studies, along with the use of advanced scientific computing tools – such as parallelization and more precise algorithms – will undoubtedly enhance this approach.

Footnotes

Additional numerical results for images affected by impulse noise

Table 4

$L^{1}$ -error between the original image f and the reconstruction u (built from $f_{δ}$ ) with $15 %$ of total pixels saved.

Noise		L1-ADJ-T		L1-ADJ-H		H1-T	H1-H
Salt	Pepper	α	$‖ f - u ‖_{1}$	α	$‖ f - u ‖_{1}$	$‖ f - u ‖_{1}$	$‖ f - u ‖_{1}$
0	0	0.01	2004.52	1.16	631.49	3069.92	629.18
0.02	0	0.46	1997.46	0.61	901.53	3681.41	2169.59
0	0.02	0.41	1606.48	0.71	928.76	3453.97	2263.76
0.01	0.01	0.61	1629.86	0.61	993.84	2936.01	2057.86
0.04	0	0.41	2623.90	0.61	1132.25	7601.93	5195.02
0	0.04	0.36	2188.72	0.71	1157.62	7299.16	5652.56
0.02	0.02	0.51	2159.01	0.56	1599.72	5063.83	4101.42
0.1	0	0.31	6111.74	0.56	1455.41	21270.53	16461.51
0	0.1	0.36	4011.69	0.51	1472.71	26678.66	20173.49
0.05	0.05	0.41	4962.55	0.51	4018.63	11612.23	10171.57

Additional numerical results for images affected by Gaussian noise

Table 6

$L^{2}$ -error between the original image f and the reconstruction u (built from $f_{δ}$ ) with $15 %$ of total pixels saved.

Noise	L2-ADJ-T		L2-ADJ-H		H1-T	H1-H
σ	α	$‖ f - u ‖_{2}$	α	$‖ f - u ‖_{2}$	$‖ f - u ‖_{2}$	$‖ f - u ‖_{2}$
0	0.01	11.62	0.01	6.58	18.21	3.35
0.03	0.71	8.25	0.56	7.69	8.14	7.64
0.05	0.51	12.79	0.66	12.23	13.00	11.91
0.1	0.31	25.95	0.76	24.32	25.89	22.98
0.2	0.01	51.20	1.11	46.93	54.71	43.29

References

Abramowitz

Stegun

I.A.

, Handbook of Mathematical Functions: With Formulas, Graphs and Mathematical Tables, Dover Books on Advanced Mathematics, Dover, New York, 1972, unabridged, unaltered and corr. republ. of the 1964 ed edition.

Amstutz

, Sensitivity analysis with respect to a local perturbation of the material property, Asymptotic Analysis 49 (2006).

Andris

Peter

Weickert

, A proof-of-concept framework for PDE-based video compression, in: 2016 Picture Coding Symposium (PCS), 2016, pp. 1–5, ISSN: 2472-7822. doi:10.1109/PCS.2016.7906362.

Bae

Weickert

, Partial differential equations for interpolation and compression of surfaces, in: Mathematical Methods for Curves and Surfaces, Dæhlen

Floater

Lyche

Merrien

J.-L.

Mørken

Schumaker

L.L.

, eds, Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 2010, pp. 1–14. doi:10.1007/978-3-642-11620-9_1.

Belhachmi

Ben Abda

Meftahi

, Topology optimization method with respect to the insertion of small coated inclusion, Asymptot. Anal. (2018). doi:10.3233/ASY-171441.

Belhachmi

Bucur

Burgeth

Weickert

, How to choose interpolation data in images, SIAM Journal of Applied Mathematics 70 (2009), 333–352. doi:10.1137/080716396.

Belhachmi

Jacumin

, Optimal interpolation data for PDE-based compression of images with noise, Communications in Nonlinear Science and Numerical Simulation 109 (2022), 106278, https://www-sciencedirect-com-443.web.bisu.edu.cn/science/article/pii/S1007570422000193 . doi:10.1016/j.cnsns.2022.106278.

Bertalmio

Sapiro

Caselles

Ballester

, Image inpainting, in: Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH’00, ACM Press/Addison-Wesley Publishing Co., USA, 2000, pp. 417–424. doi:10.1145/344779.344972.

Catté

Lions

P.-L.

Morel

J.-M.

Coll

, Image selective smoothing and edge detection by nonlinear diffusion, SIAM Journal on Numerical Analysis 29(1) (1992), 182–193. doi:10.1137/0729012.

10.

Cea

Gioan

Michel

, Quelques resultats sur l’identification de domaines, CALCOLO 10(3) (1973), 207–232. doi:10.1007/BF02575843.

11.

Chan

Shen

, Nontexture inpainting by curvature-driven diffusions, Journal of Visual Communication and Image Representation 12 (2001), 436–449. doi:10.1006/jvci.2001.0487.

12.

Dirk Adam

Peter

Weickert

, Denoising by inpainting, in: Scale Space and Variational Methods in Computer Vision, Lauze

Dong

Bjorholm Dahl

, eds, Lecture Notes in Computer Science, Springer International Publishing, Cham, 2017, pp. 121–132. doi:10.1007/978-3-319-58771-4_10.

13.

Galić

Weickert

Welk

Bruhn

Belyaev

Seidel

H.-P.

, Image compression with anisotropic diffusion, Journal of Mathematical Imaging and Vision 31 (2008), 255–269. doi:10.1007/s10851-008-0087-0.

14.

Garreau

Guillaume

Masmoudi

, The topological asymptotic for PDE systems: The elasticity case, SIAM J. Control and Optimization 39 (2001), 1756–1778. doi:10.1137/S0363012900369538.

15.

Guillaume

Idris

, The topological asymptotic expansion for the Dirichlet problem, SIAM J. Control and Optimization 41 (2002), 1042–1072. doi:10.1137/S0363012901384193.

16.

Jacumin

, thomasjacumin/vcodec-lp: v2, September 2024. doi:10.5281/ZENODO.13825644.

17.

Larnier

Fehrenbach

Masmoudi

, The Topological Gradient Method: From Optimal Design to Image Processing, Milan Journal of Mathematics 80 (2012). doi:10.1007/s00032-012-0196-5.

18.

Lenzen

Scherzer

, Partial differential equations for zooming, deinterlacing and dejittering, International Journal of Computer Vision 92 (2011), 162–176. doi:10.1007/s11263-010-0326-x.

19.

Masnou

Morel

J.-M.

, Level lines based disocclusion, in: Proceedings 1998 International Conference on Image Processing. ICIP98 (Cat. No. 98CB36269), Vol. 3, 1998, pp. 259–263. doi:10.1109/ICIP.1998.999016.

20.

Mohideen Kaja Mohideen

Peter

Alt

Weickert

Scheer

, Compressing Colour Images with Joint Inpainting and Prediction, 2020. arXiv:2010.09866.

21.

Nikolova

, Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers, SIAM J. Numerical Analysis 40 (2002), 965–994. doi:10.1137/S0036142901389165.

22.

Nikolova

, A Variational Approach to Remove Outliers and Impulse Noise. Journal of Mathematical Imaging and Vision, 20, January 2004 (2004). doi:10.1023/B:JMIV.0000011920.58935.9C.

23.

Oldham

Myland

Spanier

, An Atlas of Functions, Springer, US, 2009. doi:10.1007/978-0-387-48807-3.

24.

Pascal

Weickert

, Colour image compression with anisotropic diffusion, in: 2014 IEEE International Conference on Image Processing, ICIP 2014, 2015, pp. 4822–4826. doi:10.1109/ICIP.2014.7025977.

25.

Perona

Malik

, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence 12(7) (1990), 629–639. doi:10.1109/34.56205.

26.

Samet

Amstutz

Masmoudi

, The topological asymptotic for the Helmholtz equation, SIAM Journal on Control and Optimization 42(5) (2003), 1523–1544. http://epubs.siam.org/doi/10.1137/S0363012902406801 . doi:10.1137/S0363012902406801.

27.

Scherzer

Grasmair

Grossauer

Haltmeier

Lenzen

, Variational Methods in Imaging, Springer Science & Business Media, 2008.

28.

Schmaltz

Weickert

Bruhn

, Beating the quality of JPEG 2000 with anisotropic diffusion, in: Pattern Recognition, Denzler

Notni

Süße

, eds, Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 2009, pp. 452–461. doi:10.1007/978-3-642-03798-6_46.

29.

Ulichney

, Digital Halftoning, MIT Press, Cambridge, MA, USA, 1987.

30.

Weickert

, Theoretical foundations of anisotropic diffusion in image processing, Computing, Suppl 11 (1996), 221–236. doi:10.1007/978-3-7091-6586-7_13.

31.

Weickert

Ishikawa

Imiya

, Linear scale-space has first been proposed in Japan, Journal of Mathematical Imaging and Vision 10(3) (1999), 237–252. doi:10.1023/A:1008344623873.

32.

Ziemer

W.P.

, Extremal length and p-capacity, Michigan Mathematical Journal 16(1) (1969), 43–51, https://projecteuclid.org/journals/michigan-mathematical-journal/volume-16/issue-1/Extremal-length-and-p-capacity/10.1307/mmj/1029000164.full . doi:10.1307/mmj/1029000164.