{'fraction_non_alphanumeric': 0.05848143920865265, 'fraction_numerical': 0.020408717865101365, 'mean_word_length': 3.685979880300522, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The topological gradient method for semi-linear problems and application to edge detection and noise removal

The goal of this paper is to apply the topological gradient method to edge detection and noise removal for images degraded by various noises and blurs. First applied to edge detection for images degraded by a Gaussian noise, we propose here to extend the method to blurred images contaminated either by a multiplicative noise of gamma law or to blurred Poissonian images. We compute, both for perforated and cracked domains, the topological gradient for each noise model. Then we present an edge detection/restoration algorithm based on this notion and we apply it to the two degradation models previously described. We compare our method with other classical variational approaches such that the Mumford-Shah and TV restoration models and with the classical Canny edge detector. Some experimental results showing the efficiency, the robustness and the rapidity of the method are presented.


1.
Introduction. An important problem in image analysis is the reconstruction of an original image u from an observed image f . In general this includes restoration and edge detection processes. The transformation between f and u originates from two phenomena. The first phenomenon is related to the acquisition process (blur created by a wrong lens adjustement or by a movement, Poissonian photons emission rates ...) and the second is due to the signal transmission. A lot of methods to reconstruct such degraded images exist: stochastic methods [23,14], wavelets decomposition [32,19], morphological methods [38]. Here we are interested with variational approaches [9]. In this context, the most famous model is the Mumford-Shah functional [34] but other works based on variational methods do exist [9]. Among more recent papers, we can cite [6] for restoration of images contaminated by speckle noise, [17] for blind restoration of Poissonian images, and [37] for an overview of image restoration degraded by different type of noise.
In this paper we tackle the edge detection problem by using the topological gradient method. First introduced for cracks detection by Sokolowski et al. [39] and Masmoudi [33], and applied in optimal design and mechanics ( [3,2,5]), this notion consists in the study of the variations of a cost function j(Ω) = J(Ω, u Ω ) with respect to a topological variation, where J(Ω, u) is of the form J(Ω, u) = Ω F (u, ∇u, ∇ 2 u, . . . ) and u Ω is a solution of a PDE defined on the image domain Ω. In order to compute the topological gradient, we remove from Ω a small object ω of size → 0 centered at a point x 0 ∈ Ω (generally a ball or a segment) and we set Ω = Ω\ω . Two typical examples are: for small > 0 (a) Ω = Ω\{x 0 + B} and (b) Ω = Ω\{x 0 + σ(n)}, where B = B(O, 1) is the unit ball of R 2 and σ(n) is a straight segment with normal n (a crack). We compute I(x 0 ) = lim →0 where ρ( ) is a non negative function such that ρ( ) → 0 if → 0 (in our context, ρ( ) = 2 ). I(x 0 ) is called the topological gradient at x 0 . It measures the energy contained by a perturbation centered at x 0 . The type of structure to be detected depends on the choice of the cost function J(Ω, u). Recently this notion has been used in image processing and to the best of our knowledge the first works in this direction are those by [26,31,15]. Then other imaging problems such as inpainting, classification, demosaicing, super resolution, have been adressed using a topological gradient approach [10,16,11,13,12,29,30]. In [15], only Gaussian additive noise is considered and in [30] blur has been introduced. In fact, in [30] more general degradations have been taken into account. The authors consider model of the form f = Lu + b where f is the observed image, L a linear operator and b a Gaussian additive noise. They compute in the case of a crack the topological gradient and illustrate their approach for various imaging problems (edge detection/restoration, super resolution, demosaicing). Note also that topological gradient methods have been also applied for fine structures detection (e.g. points and filaments) [8,7,21]. In this case the cost function is based on second order derivatives. Edge detection/restoration in imaging are in general ill-posed inverse problems and one way to overcome this difficulty is to regularize them. A classical framework to do that is to use a Bayesian formulation which leads to the minimization of an energy consisting in two terms. The first one is a data fidelity term which takes into account both the statistic of the noise and the blur and the second one is an adequate regularizing term. For example if we suppose that the acquisition model is of the form f = u + b where b is Gaussian noise then an anti-log-likelihood estimator amounts to choose as a data fidelity term the L 2 -norm u − f 2 L 2 (Ω) . If the noise follows another statistic, of course this term changes. The regularizing term is often based on an L p norm of the gradient. Our main contribution is to generalize the results given in [30] to blurred images contaminated by Poissonian statistic and images degraded by speckle noise. We give the different expressions of the topological gradient associated to the cost function j . We suppose that Ω is the disjoint union of (R N0 j ) j∈{1,...,N0} . The operator is a convolution (generally positive and such that K1 = 0) representing the blur. We denote by K its adjoint. The functions ψ(x, u) and ψ j (v) will be specified in section 4 and section 5. Note that problems (2) and (3) are semi-linear and one of our contribution is to show they are well-posed and verify some maximum principles. Speckle noise is a multiplicative noise of gamma law, which is present in SAR images, laser images, microscope images [28,25,40]. A Poisson statistic occurs in confocal microscopy [20], emission tomography [41] and single-photon emission computed tomography [24].
In section 2, we recall the classical rationale justifying the modelization of the data fidelity term in a Bayesian approach. In section 3 we set the variational problem taking into account the blur. In section 4 we show that problem (2) is well-posed and give the associated topological gradient both for perforated and cracked domains (in fact, we study a more general class of problem (2)). In section 5 we treat the Poissonian model (3) whose energy is not standard. We summarize in Table 1 (section 6) all the expressions of the topological gradient according to the type of noise and to the infinitesimal perturbation. In section 7 we show how to apply the notion of topological gradient to restore degraded images. Finally in section 8, we present for all the models, the way to numerically implement the computation of the topological gradient and we display various experimental results illustrating each of them. We conclude this section by giving some notations and assumptions.

Notations and assumptions:
• Ω ⊂ R 2 is the image domain.
• u H 1 (Ω)/R the norm on the quotient space H 1 (Ω)/R. • B r (x) the ball centered at x and of radius r, B r the ball centered at 0 and of radius r and B = B 1 denotes the unit ball.
• f : Ω → R is the observed image, u : Ω → R the image to be recovered.
• R N j models pixel j, where N is the number of pixels. We suppose that Ω is the disjoint union of the R j 's.
Only the proof for a perforated domain (a) is performed since for a cracked domain (b) the explicit dependency on the data is killed by the fact that the crack has a null Lebesgue measure.

2.
A Bayesian approach. In this section we recall, in the discrete setting, the classical Bayesian approach allowing to deduce the suitable variational model for restoring noisy images. The reasoning is as follows: we express the a priori conditional probability P (u N = u|f N = f ) and then we look for u N as the value maximizing P (u N = u|f N = f ) (a Maximum A Priori (MAP) estimator). Thanks to the Bayes rule, we have: depends on the noise model and P (u N = u) is an a priori density probability. Writing that u N is a minimum of −log(P (u N = u|f N = f )) leads to looking for u N as the solution of The a priori density P (u N = u) plays the role of a regularizing term. In analogy to statistical mechanics, a priori densities are frequently Gibbs functions [23] of the form: where G N (u) is a discrete version of a non negative energy functional G(u) and C is a constant. Here we will consider that G N (u) = J N (Ω, u) is a discretization of J(Ω, u) defined in (1). The choice of the density probability P (f N = f |u N = u) depends on the statistic of the model to be considered. Below we review respectively speckle and Poisson models.
(i) Speckle model For SAR images, the classical modeling is (see [40]): f N = s N u N where u N is the reflectance of the scene (which is to be recovered) and s N the speckle noise. Let us explicit the law of s N . SAR images are constructed from L ∈ N observations f k,N with 1 ≤ k ≤ L and for each observations we have f k,N = h k,N u N . Generally h k,N is a random variable which follows a negative exponential law with mean 1. Then, the observed image f N is construct from this L observations as: is 1 and its variance 1 L ). Standard computations lead to Classically, it is assumed that f N and u N are independently distributed.
Hence the density of f N |u N is the product of the densities f N j |u N j 1≤j≤N . By applying the −log function we deduce that (4) in this case rewrites as (up to a multiplicative constant) for u ∈ R N and u > 0, where C denotes a constant independent of u. The factor L can be neglected since it can be scaled with the constant γ. Passing to the limit as N → ∞, we deduce the following continuous energy (ii) Speckle with log of the image (Speckle-log model) One drawback of (6) is that the data fidelity term is nonconvex. By setting, v N = log(u N ) and g N = log(f N ), we can still use the Bayesian reasonning with v N = g N + log(s N ). The expression of the data fidelity term is changed and taking the same Gibbs function (G N (v) = J N (Ω, v)) we deduce the continuous energy: which is now a convex function of v. The recovered image is then u = e v .

(iii) Poissonian model
This model is classical in astronomical and confocal microscopy imaging [20]. Poissonian statistic is due to the stochastic nature of photons counts at sensors. We set N = N 0 the number of sensors which are modelled by the R N0 j 's. To simplify notation, since N 0 is fixed we drop the superscript N 0 in all variables and we assume that |R j | = 1 ∀j ∈ {1, .., N 0 }. We assume that f is a piecewise constant function equal to f j on each R j . The observed image f j is a realization of a Poisson statistic of mean and variance equal to is the continuous image to be recovered. Thanks to the independence of f j and λ j , the conditional probability P (f |λ) is given by: and by applying the −log function, we have: where C is a constant independent of u. We deduce the energy: 3. Blurring modelling. In most imaging applications the optical material, the motion of the camera or of the target introduce a blur on the observed image (see [36]). Generally spatially invariant, the blur is modelled as a positive convolution operator u → Ku with K1 = 0. We denote by K N the N × N matrix associated to a discrete version of K on Ω. From section 2 we can deduce the following models adapted to each kind of noise and taking into account the blur: 1. Speckle model: the observed discrete image writes as f N = s N K N u N and the energy is 2. Speckle-log model. We recall that the model writes as The deblurring cannot be handled simultaneously with the denoising process since in general the equation v N = log K N u N is not invertible. In this case, it is then preferable to use (9). 3. Poissonian model: the observed image f j is a realization of a Poissonian statistic of mean Rj Ku N (x)dx, so the energy is 4. Speckle multiplicative noise. In this section we show that equation (2) admits a unique solution for a more general class of functions ψ including the speckle-log model. For that we consider the minimization problem: 4.1. Well-posedness of problem (2). We begin by stating the hypotheses defining the class of functions ψ(x, u) that we consider, including the speckle-log model (11).
Hypotheses 1. Let I be an interval of R and ψ : Ω × I −→ R such that Lemma 4.1. Let ψ(x, u) a function verifying Hypotheses 1, we assume that K is non decreasing and that K1 = 0, then (11) admits a unique solution u Ω ∈ H 1 (Ω). Moreover, we have a K1 ≤ u Ω ≤ b K1 . Proof. To shorten notations, when it is not ambiguous, we denote ψ(u) for ψ(x, u).
Existence. Let (u n ) a minimizing sequence. There exists a constant C 1 such that E(Ω, u n ) ≤ C 1 . As ψ(x, u) is bounded from below on Ω × I there exists a constant C 2 such that Ω ψ(x, u n ) ≥ C 2 . Therefore: Let be a = a K1 (s.t. Ka = a), v n = max(u n , a ) and Ω − n = Ω ∩ {u n ≤ a }, we have v n ≥ a and By convexity ψ(Ku n ) − ψ(Ka ) ≥ D u ψ(a)(Ku n − Ka ), and from Hypothesis 1-(v) and using that K is not decreasing we get Ω − n ψ(a) − ψ(Ku n ) ≤ 0. We deduce that E(Ω, v n ) ≤ E(Ω, u n ). Thus v n is still a minimizing sequence and v n ≥ a . Similarly by setting b = b K1 and w n = min(v n , b ), we get w n ≤ b and w n is a minimizing sequence. Therefore we can suppose that any minimizing sequence u n verifies a ≤ u n ≤ b . Thus u n is bounded in H 1 (Ω) and up to a subsequence there stands for the weak topology). By using the lower semi-continuity of J(Ω, u) and Fatou's Lemma we obtain that u is a solution of (11). Moreover we have a ≤ u ≤ b a.e. on Ω.
Uniqueness. From the existence, we can work on the set (Ω), we deduce that E(Ω, u) is strictly convex on H(Ω) and admits a unique minimum.
Since problem (11) and equation (2) are equivalent we deduce: Under the same hypotheses as in Lemma 4.1, (2) admits a unique solution We apply below Lemma 4.1 to the speckle-log model for K = I.

Proof. A standard computation leads to
Let a = log(α) and b = log(β), the following inequalities hold From Lemma 4.1, there exists a unique function u ∈ H 1 (Ω) solution of (11). Moreover we have a ≤ u ≤ b.
In the next subsection we give the topological gradients for perforated and cracked domains when K = I and ψ satisfying Hypotheses 1.

Computation of the topological gradient for a perforated domain.
For a perforated domain Ω = Ω\{x 0 + B}, let u be the solution of the problem from the variational formulation of the equation verified by u , computing the topological gradient associated to J(Ω, u) is equivalent to compute the one associated to J(Ω, u). To give the topological gradient expression we need (see [2,22,39]) to introduce an adjoint solution v defined by (13) where u 0 denotes the solution of (12) for = 0 and ∂ J(u 0 ) the gradient of J(Ω, .) at u 0 .
Proof. The proof of Theorem 4.2 is given in [27,22] under slightly different assumptions; see also [22]. Remark 1. The adjoint problem (13) is linear and we can notice that the strict convexity of u → ψ(x, u) is necessary to (Q ) be coercive. Since u → ψ(x, u) is C 2 (I) and thanks to Lemma 4.1 there exist two constants A, B ∈ R such that (13) is well-posed. Moreover, we have the following inequality Expression of the topological gradient for a cracked domain. For the cracked domain Ω = Ω\x 0 + σ(n) the data fidelity term does not affect the topological gradient expression since the crack has a 2d Lebesgue measure equal to 0. The asymptotic expansion of u and v are similar and the computation of the topological gradient is the same as in the linear case and is given in [5,22]. Theorem 4.3. The topological gradient I c Lap associated to problem (12) and to the cost function J (u) = Ω |∇u| 2 for a cracked domain is with u 0 and v 0 respectively solution of (12) and (13) for = 0.

Poissonian model with blurring.
In this section we consider the Poissonian model (3) whose associated energy is given by: . First we show that problem (16) is well-posed, then we compute the topological gradient for a perforated domain and we give the expression for a cracked domain.
it is straightforward that ψ i Ri Ku n ≤ D − K i and then Then we write u n as u n = v n + w n where w n = 1 |Ω| Ω u n . Since ∇v n = ∇u n and Ω v n = 0, from Poincaré-Wirtinger, there exists F ≥ 0 such that v n L 2 (Ω) ≤ F . Moreover, we have Thus, since K1 = 0, we get that Ω u n is bounded and thus u n is bounded in L 2 (Ω). So, there exist a sub-sequence u n k (still denoted u n ) and u ∈ H 1 (Ω) such that u n L 2 (Ω) → u and u n H 1 (Ω) u. We deduce that Ω |∇u| 2 ≤ lim inf Ω |∇u n | 2 and Ri Ku n → Ri Ku. By continuity, we have which proves that u is a minimizer of E(Ω, u).
Bounds. If u is the solution of (16) then Uniqueness. From the two previous bounds we can minimize E(Ω, u) over the .., N 0 }, we get that u → ψ j Rj Ku is strictly convex on H(Ω). As |u| 2 1,Ω is strictly convex, E(Ω, .) is strictly convex on V β α (Ω) and then E(Ω, .) has a unique minimum.
Since problem (16) and equation (3)  and β = Ω f = N0 j=1 f j . Remark 2. We can show that Proposition 2 holds as soon as ψ j is bounded from below for j ∈ {1, ..., N 0 } and strictly convex on I ⊂ R. In the general case α and β are implicitly defined in function of ψ j .

5.2.
Computation of the topological gradient for a perforated domain. We suppose that x 0 ∈ Ω\ ∪ j ∂R j and for sufficiently small we set j 0 the index such that R j0 ⊃ B (x 0 ). For j ∈ {1, ..., N 0 }, let R j be the domain equal to R j0 \B (x 0 ) if j = j 0 and R j otherwise. Finally, let u be the unique solution of problem (3) replacing Ω by Ω : where I j (u) = R j u. As for the speckle-log model, by setting J(Ω, u) = − N0 j=1 ψ j (I j (Ku))I j (Ku), from the variational formulation of the equation verified by u , computing the topological gradient associated to J(Ω, u) is equivalent to compute the one associated to J(Ω, u). To compute the topological gradient we need to introduce an adjoint solution v defined by: We now give the main result of this section.
Theorem 5.1. The topological gradient I b Lap associated to problem (19) and to the cost function J (u) = Ω |∇u| 2 for a perforated domain is with u 0 and v 0 are respectively solution of (19) and (20) for = 0.
Proof. We split the proof in three steps and we assume without loss of generality that x 0 = 0.
Step 1. This step is dedicated to the computation of the difference J (u ) − J 0 (u 0 ). Using the variational formulation of u : find u ∈ H 1 (Ω ) s.t.
we have J (u ) = J (u ) with: Using a Taylor expansion of φ j (s) at I j (Ku 0 ) and by remarking that Step 2. In this step we introduce the adjoint problem in order to estimate the linear term L (u − u 0 ) in (23). For that we need to linearize (22) w.r.t u. A Taylor expansion gives: As a consequence we have

Kv
Let us introduce now the adjoint problem: It is easily seen that (26) is the variational formulation of (20). From (23) and (25), we obtain By using an integration by parts and equation (19) for = 0, A (u 0 , v ) expresses as In the sequel, we need to extend v inside the ball B . To do that, for ϕ ∈ H 1/2 (∂B ), we denote by l ϕ the harmonic function defined by Step 3. In this step we give the asymptotic expansion of all the previous terms. and Proof. The first equality is straightforward. A Taylor expansion of u 0 at 0 gives the first estimation. Again a Taylor expansion of u 0 at 0, Lemma 9.4 (see Appendix A) and the fact that R j0 →0 −→ R j0 give the second estimation. For K , from Lemmas 9.1 and 9.4 we have: For F 2 it suffices to make a change of variable and use the trace theorem on B 2 \B: Now from the equivalence of the H 1 (B 2 \B)/R-norm with the semi-norm and a change of variable, we get ∂ n l e ( X) −1/2,∂B ≤ C|e | 1,Ω . By using Lemma 9.4, we obtain F 2 = O( 3 −log( )). For E 1 and E 2 , we use Lemma 9.3, the boundedness of Ω Kv independently from and the fact that η and ξ j are in ]α, β[ : By using that ∆v 0 = ψ j0 (I j0 (Ku 0 ))I j0 (Kv 0 ) − K1 and a Taylor expansion of u 0 at 0 we get . For E 4 , by using a change of variable, the continuity of ϕ → l ϕ from H 1/2 (∂B) to H 1 (B), the trace theorem on B 2 \B, again a change of variable and Lemma 9.4 we obtain: Remark 3. The topological gradient is the same in the general case of functions ψ j ∈ C 3 (I), strictly convex on an interval I and bounded from below on I. Just in the right hand-side of (20), K 1 must be replaced by φ j (I j (Ku 0 ))K 1 with φ j (s) = ψ j (s)s. 6. Summary table of the topological gradient expressions. We summarize in Table 1, all the expressions of the topological gradient according to the type of noise and to the infinitesimal perturbation.
Ball Crack Table 1. Summary of the topological gradient expressions 7. Restoration using the topological gradient for a cracked domain. As a by product the computation of the topological gradient I c Lap for a cracked domain allows to restore images degraded by blur or/and various noise statistics. Once I c Lap is computed, we define for a fixed threshold δ > 0, the set E δ = x ∈ Ω; |I c Lap (x)| ≥ δ and the approximated characteristic function where η is a small positive parameter. Typically, we take η = 10 −5 in the numerical experiments. From the computation of I c Lap we also get the normalized direction τ = n ⊥ of the edge. If n = (cos(ϕ), sin(ϕ)) is the normal to the crack given by I c Lap , we have τ = (sin(ϕ), −cos(ϕ)). Then, if f is the degraded observed image, we want to find a restored version u ∈ H 1 (Ω) of f as the solution of the following anisotropic PDE: and where P ϕ η (x) is a tensor constructed from ϕ(x) and χ η (x) and γ is a parameter to tune. More precisely, we choose P ϕ η (x)∇u(x) = (∇u.τ )τ + χ η (x)(∇u.n)n. A simple identification shows that P ϕ η (x) is the matrix where n 1 = cos(ϕ(x)) and n 2 = sin(ϕ(x)). The interpretation of this matrix P ϕ η (x) is as follows: (i) if x belongs to the background, thanks to the definition of χ η (x), P ϕ η (x) writes as P ϕ η (x) = I, so div(P ϕ η (x)∇u) = ∆u and the smoothing is isotropic.
(ii) if x belongs to an edge (i.e. x ∈ E δ ), then χ η (x) is close to zero and P ϕ η (x)∇u(x) ≈ (∇u.τ )τ and the diffusion is in the direction of the edge. As we will see in section 8 on numerical examples, the restoration results obtained when applying equation (32) are very good. To compute a numerical solution of equation (32), we consider the associated strictly convex energy : E(u) = Ω P ϕ η ∇u.∇u + ψ(x, u) that we discretize by using a classical finite differences scheme. Then we compute the minimizer of the discrete energy by applying a Scale Gradient Projection (SGP) type algorithm described in Algorithm 2 below. 8. Numerical application to 2D imaging. In this section we illustrate the theory of the topological gradient by giving various experimental results for models (2) and (3).
The topological gradient expressions for the two models are stated in sections 4 and 5 and are summarized in section 6 ( Table 1).
For each model, to compute the topological gradient (TG) we use Algorithm 1. The computation of the direct and adjoint solutions is specific to each model.
We first perform the discretization of problems (11) and (16) and then we give the experimental results. As the adjoint problems (13) and (20) are linear with non constant coefficients we discretize them by a finite difference scheme and we compute the discrete solution by using a sparse solver. 8.1. Solving numerically Poisson and Speckle-log models. In this section, we assume that f and u are vectors of R N and min f > 0. By Proposition 1 and Proposition 2, problems (11) and (16) can be discretized as: where α s = min f and α p = min f N ; J p (u) and J s (u) are respectively the discrete versions of energies (10) and (7).
During the construction of a minimizing sequence u (k) , the condition u (k) ≥ α p for the Poisson model (respectively u (k) ≥ log(α s ) for the speckle model) must be fulfilled at each step. Hence a projection ensures this condition. To solve these problems we use an iterative algorithm based on the SGP algorithm [18]. Let us write the discrete cost functions: where A is the Neumann Laplacian matrix, K is a discretization of the blurring operator (circulant block matrix if we assume that the image is periodic) and we recall that g i = log(f i ). Let us give the main ideas of the SGP algorithm. The discrete energies J s and J p are denoted by J as soon as we do not use their expression and δ will be the number equal to α for the Poisson model and equal to log(α) for the Speckle-log model. We set Λ = u ∈ R N , u ≥ δ . We want to find u ∈ Λ such that ∇J(u ) = 0. At step k, a first order Taylor expansion at point u = u (k) leads to the following equation . We deduce by this reasoning that the direction of the descent algorithm can be given by ∇ 2 J(u (k) ) −1 ∇J(u (k) ), but the computation of this direction is very costly. We denote by D L the compact set of symmetric positive definite N × N matrices such that D ≤ L and D −1 ≤ 1 L . The main idea of the SGP algorithm is to construct two sequences α k and D k ∈ D L such that α k D k approximates ∇ 2 J(u (k) ) and to project each iterate on Λ with respect to the norm u D = √ u T Du. We set P Λ,D −1 for D ∈ D L the projector on Λ with respect to the norm . D . For more details on Algorithm 2, we refer the reader to [18] and we only give some explanations on the construction of the sequences D k and α k . Choose the parameter α k ∈ [α min , α max ] and the scaling matrix D k ∈ D L 4: if y (k) = u (k) then 6: Stop, u (k) is a stationary point. λ k ← 1 and J max ← max 0≤j≤min(k,M −1) J(u (k−j) ) 10: λ k fixed by backtracking: 11: while f (u (k) + λ k d (k) ) ≤ J max + βλ k ∇J(u (k) ) T d (k) do 12: λ k ← θλ k 13: end while 14: u (k+1) ← u (k) + λ k d (k) 15: end for The approximation of the Hessian matrix ∇ 2 J(u (k) ) is B(α k ) = α k D k . By using a first order Taylor expansion of ∇J(u) at point u (k−1) we get that Hence two possible choices of α k can be made: where s (k−1) = u (k) − u (k−1) and z (k−1) = ∇J(u (k) ) − ∇J(u (k−1) ). In [18] the choice of α k is the output of an algorithm called SGP-SS Algorithm (SGP step length selection) which uses two thresholds 0 < α min < α max . The derivative of the discrete cost functions J p and J s are: where 1 ∈ R N denotes the vector whose each coordinate is equal to 1, diag(x) for x ∈ R N is the diagonal matrix with diagonal entries equal to x. For x ∈ R N and ϕ : R −→ R, ϕ(x) stands for the vector (ϕ(x i )) 1≤i≤N . The choice of the parameters in Algorithm 2 is the following: β = 10 −4 , θ = 0.4, k max = 600, M = 1 and for the Poisson model (16) we set α min = 10 −10 , α max = 10 5 while for the Speckle-log model (11) we set α min = 10 −5 , α max = 10 15 . The initial value of u (0) is either the observed image for the Poisson model or its logarithm for the Speckle-log model. Let us note that in the case of equation (32), the matrix A is the finite differences discretisation of the operator div(P ϕ η ∇.) and for the speckle model of restoration we use the function ψ(u) = log(u) + f u . 8.2. Comparison of our method with some classical models. We compare the results given by the topological gradient approach with the ones performed by the Mumford-Shah and TV models. We will also compare our detection process with the well-known Canny edge detector.

Mumford-Shah model of segmentation/restoration and its approximation
Let u be the image of support Ω, the functional introduced by Mumford and Shah in 1989 (see [34]) is: where f is the observed image, u is a function defined on Ω\γ (the restored version of f ) and γ ⊂ Ω is the set of discontinuity of u. H 1 is the Hausdorff measure of γ, λ and α are positive parameters. The difficulty is that the unknown are not of same nature: u is a function and γ is a set. Ambrosio and Tortorelli [1] proposed an approximation of this functional as follows: We will change the data fidelity term |u−f | 2 according to the a priori model (speckle and Poissonian models) i.e. we will compare (32) with the solution of (see [37]): where K is the blur operator, u(x) is the restored image, 1 − b(x) ≈ 0 is the characteristic function of the edges and ψ(x, u) is given in (33). We will call this model the Mumford-Shah model (MS).

TV model of restoration
The TV model is well-known: we search for a restored version u minimizing an energy functional which is the sum of the total variation ||Du|| and a data fidelity term which depends on the a priori model (see [6] for the speckle model and [17] for the Poisson one). Thus, we will compare model (32) to: (37) min where λ is a parameter, K the blur operator, and ψ(u) is given in (33). In the sequel we call this model the TV model.

The Canny edge detector
We will also compare the edge detection performed by our approach with the Canny edge detector which is the norm of the gradient of a regularization of the image by a Gaussian convolution at scale σ.
For more details on these models we refer the reader to [17,6,37,34] and for the numerical approximations to [9] and the references therein.
For restoration comparisons on synthetic images we use two indicators: • the PSNR which is defined for a noisy observation I of an image I 0 by P SN R(I) = 10Log 255 2 I − I 0 2 2 /N where N is the number of pixels in the image.
• the SSIM defined as where µ x stands for the mean of x, cov is the covariance operator and c 1 and c 2 are constants given for RGB images by c 1 = (255k 1 ) 2 and c 2 = (255k 2 ) 2 with k 1 = 0.01 and k 2 = 0.03.

8.3.
Numerical results for the speckle-log model. Comparisons. In this section we illustrate the detection process given by I c Lap (15), I b Lap (14), b(x) (computed with the MS model (36)) and the Canny edge detector. All edge detectors are displayed up to a normalisation in order to have edges in black and the background of the image in white. We also display the restoration performed by (32), the MS model (36) and the TV model (37).
On Fig. 1, we compare I c Lap (15) and I b Lap (14) for different values of γ with b(x) and with the Canny edge detector for a synthetic speckled image. I b Lap seems more adapted to detect isotropic structures and we deduce that γ must be tuned with respect to the noise and to the size of structures to detect. Comparing I c Lap , MS and the Canny edge detector, similar results are obtained for the cheetah.  (14) and (15) for different values of γ with the MS edge detection (36) and the Canny edge detector for an initial speckled image (L = 6) containing mainly isotropic small structures.
Remark 5. Let us note that for large γ an edge doubling phenomena appears for I c Lap on Fig. 1-(e). As shown in [4], in the linear case (Gaussian noise), the edge doubling can be avoided by using as cost function the opposite of the square of the L 2 -norm. More generally we have remarked that the edge doubling can be avoided for the speckle-log model by adding in the cost function the term Ω ψ(u). Let us recall that the direct and adjoint models given in (12) and (13) are associated to the auxiliary cost function J(Ω, u) = − Ω D u ψ(x, u)u. Fig. 2 shows profiles of the direct and adjoint models and I c Lap (15) associated to J(Ω, u) and to J 2 (Ω, u) = J(Ω, u) + Ω ψ(u).
The result given in Fig. 3 for a real SAR image is similar to the one of  (14) and (15) with the MS edge detection (36) and the Canny edge detector for a real SAR image. and the other indicators detect its edges. Let us notice that I c Lap detects only edges while the Canny edge detectors is also sensitive to texture. Fig. 4 and Fig. 6 compare the restoration performed by (32), (36) and (37) respectively on a real SAR image and on a synthetic speckled image. Restoration given by the MS degrades contours while (32) and the TV model (37) are nearly equivalent. However, on Fig. 6, we can notice that the restoration performed by (32) is better than the TV one and computation times are equivalent (see section 8.5).
On Fig. 5 we compare I c Lap , b(x) computed with MS (36) and the Canny edge detector for a very noisy synthetic image (L = 1 i.e. the worst case for this model).
Here I c Lap gives a quite good result with respect to MS where edges are spread out and the Canny edge detector which detects a lot of noise.
Finally Fig. 7 shows the 1D profiles of the image to recover, its noisy version, the restored version (32) and I c Lap (15) across an edge. The restored version matches very well the image to recover and edges are not degraded except on strong angular point with low contrast where the direction of the diffusion in (32) is not well defined and the topological gradient I c Lap is low in magnitude. This shows that (32) is a good restoration process.   (32), (36) and (37) for a real SAR image (Zoom on Fig. 3-(a)).   (32), (36) and (37) for a speckled synthetic image (L = 1).
Remark 6. We can still do a similar remark as Remark 5 for Fig. 8-(e). Here also we have remarked that the edge doubling can be avoided for Poissonian model by adding in the cost function the term N0 j=1 ψ j (I j (Ku)). Let us recall that the direct and adjoint models given in (19) and (20) are associated to the auxiliary cost function J(Ω, u) = − N0 j=1 Ω ψ (I j (Ku)I j (Ku). Fig. 9 shows profiles of the direct and adjoint models and I c Lap (31) associated to J(Ω, u) and to J 2 (Ω, u) = J(Ω, u) + N0 j=1 ψ j (I j (Ku)).  Fig. 11 and Fig. 13 display the restoration computed by (32), (36) and (37) on respectively a real confocal image and a synthetic Poissonian image blurred by Gaussian convolution. We notice that (32) and (37) restore very well the image preserving edges unlike to the MS model (36) which degrades contours and which does not annihilate the blur effect.
Finally, Fig. 14 shows the 1D profiles of the image to recover, its degraded versions (blurred, blurred + Poissonian process), the restored version (32) and I c Lap (31) across an edge. We see that (32) allows to recover the initial image and that I c Lap detects very well the edge. 8.5. Computational time comparisons for the three methods. On Fig. 15 are compared the three methods for K = I (no blur): Mumford-Shah (MS), TV and the topological gradient method (TG). The experiments are performed on a computer equipped with a processor Intel Core 1.9 GHz and all algorithms are implemented in Matlab. The computation time for the Canny edge detector is not displayed because it just consists in a convolution which is very fast and so it would be very close to the y-axis compared to other methods.
For Poisson and speckle models, the computation of the topological gradient still remains the fastest but the restoration step performed by (32) takes approximately the same time as for a restoration given by a TV model. Let us notice that the TV model is implemented by using an iterative algorithm (explicit schema with fixed step length with at maximum 1000 iterations) and that the Mumford-Shah model solution is computed by minimizing the approximate functional alternatively with respect to u and b. The minimization with respect to u (at b fixed) is made by using Algorithm 2 and we perform 15 iterations. Finally let us precise that we present on Fig. 15 the computation time only for the Poissonian model because similar results are obtained from the speckle one.

9.
Appendices. In these appendices we give the asymptotic expansion of the differences u − u 0 and v − v 0 for the non linear problems (Poisson and Speckle-log models). Some proofs are similar to the linear case and so we will refer the reader to [5]. To establish these asymptotic expansions we need the following exterior (a) Initial image where g ∈ H −1/2 (∂B) and ∂B gdσ = 0. For the computation of the topological gradient we will use the two following lemma. We omit the proofs and we refer the reader to [22] for more details  and l P expresses as l P (x) = ∂B λ(y)E(x − y)dσ.
The following asymptotic estimations holds.

1
|Ω| Ω e 3 and using Poincaré-Wirtinger inequality and the fact that K1 = 0 (as in Proposition 2) we get e 3 1,Ω ≤ C 2 . From the inequality e 1,Ω ≤ e 1 H 1 (Ω )/R + e 2 H 1 (Ω )/R + e 3 1,Ω , we obtain the estimation. For the L 2 (Ω )-norm estimation of X , it suffices to take the L 2 (Ω )-norm of its asymptotic expansion and to use the first point of Lemma 9.3 and Lemma 9.2.  (20) for > 0 and for = 0, the Euler equations associated to w are: This problem is linear and from Proposition 2 we have: Then the topological expansion of w can be deduced from the proof of Lemma 9.3 or from the linear case with constant coefficient [5]. See also Nedelec [35] for details on the analysis. The two last estimations are straightforward by using the topological expansion of w and Lemma 9.2.