The topological gradient method for semi-linear problems and application to edge detection and noise removal.

In this paper we propose a new variationnal method for segmenting/restoring images degraded by diverse noises and blurs. This method is based on the notion of topological gradient. First applied by [11] to restore images degraded by a Gaussian noise, we propose here to extend the segmen-tation/restoration process for possibly blurred images contamined either by an additive Gaussian noise, or a multiplicative noise of gamma law or in presence of Poissonian statistics. We calculate, both for perforated and cracked do-mains, the topological gradient for each noise model. Then we present a segmentation/restoration algorithm based on this notion and we apply it to the three degradation models previously described. Then, we compare our method with the Ambrosio-Tortorelli approximation of the Mumford-Shah functional [23,1]. We also compare our results with those given by a classical TV restoration process (see [4] for a speckle model). Many experimental results showing the ef-ﬁciency, the robustness and the rapidity of the approach are presented.


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 segmentation processes.
The transformation between f and u originates from two phenomena. The first phenomena 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 [16,10], wavelets decomposition [21,13], morphological methods [26]. Here we are interested with variational approaches [6]. In this context, the most famous model is the Mumford-Shah functional [23] (1989) but other works based on variational methods do exist ( [6]). Among more recent papers, we can cite [4] (2008) for the restoration of images contamined by speckle noise and [25] (2013) for the restoration of images degraded by different type of noise.
In this paper we tackle the segmentation problem by using a topological gradient method. First introduced for cracks detection by Sokolowski [27] and Masmoudi [22], 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 calculate the topological gradient, we remove to Ω a small object ω ε of size ε → 0 centered at x 0 ∈ Ω (generally a ball or a curve) 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 the limit : I (x 0 ) = lim ε→0 where d is the dimension of the ambiant space. I (x 0 ) is called the topological gradient at x 0 . It measures the energy contained by a perturbation centered at x 0 and so the structures that we want to detect correspond to the points x 0 where I (x 0 ) is the largest. 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. For detecting edges the usual choice for the cost function is J Ω (u) = 11]). Other imaging tasks such as inpainting or classification have been addressed by using this approach [8,9]. Note also that topological gradient methods have been applied for the detection of fine structures [5] such as filaments or points. In this case the cost function is based on second order derivatives. In [7,11] only Gaussian additive noise are considered and no blur has been introduced. In this paper we propose : (i) to extend the model to another kinds of noise frequent in real images and (ii) to introduce blur in the process. Restoration/segmentation 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 + n where n is Gaussian noise then an anti-log-likelihood estimator amounts to choose as a data fidelity term the L 2 norm u − f L 2 (Ω ) . If the noise follows another statistic, of course this term varies. The regularizing term is in general based on an L p norm of the gradient. The main contribution of this work is to generalize the results given in [7,11] to blurred images contamined by speckle noise and Poissonian process and to give the different expression of the topological gradient associated to the same cost function J Ω (u) = Ω |∇u| 2 . More precisely we will consider variational problems of the following form Speckle and Gaussian noise model : Poissonian model : where I ind (Ω ) is the indices set of the pixels, γ > 0 is a parameter, R j is a regular domain modeling pixel j such that Ω is the disjoint union of (R i ) i∈I ind (Ω ) , K : L 2 (Ω ) −→ L 2 (Ω ) is a convolution operator (generally positive and such that K1 = 0 ) representing the blur. The functions ψ(x, u) and ψ j (X) will be specified in section 5 and section 6. Note that the speckle noise is a multiplicative noise of gamma law. It is present in SAR images, laser images, microscope images [20,18,28]. A Poisson statistic occurs in confocal microscopy [14], emission tomography [29] and single-photon emission computed tomography [17].
In section 2, we recall the rationale for justifying the modelization of the data fidelity term in a Bayesian approach.
In section 3 we set the variational problem taking into account the blurring. Then in section 4 we give the topological gradient for a blurred and Gaussian noisy image. In section 5 we establish the topological gradient for the speckle model and for a more general variational class of problems. In section 6 we treat the particular Poisson model since the form of the energy in this case is not standard. The topological gradient is not only an edge detector but also as a by-product it allows to restore the degraded observed image. We explain how this is possible in section 7. Finally in section 8 we present, for all the models, the numerical analysis for computing the topological gradient and we display some experimental results for each of them.
We conclude this section by giving some notations and assumptions : we suppose to simplify that x 0 = 0 and we denote by u m,Ω the H m (Ω )-norm of the Sobolev space H m (Ω ) = u, D α u ∈ L 2 (Ω ), |α| ≤ m and by u H 1 (Ω )/R the norm on the quotient space H 1 (Ω )/R. We set J Ω (u) = Ω |∇u| 2 and J Ω ε (u) = J ε (u).
Only the proof for a perforated domain is performed since it is more interesting than for a cracked domain in which the explicit dependency on the data is killed by the fact that the crack has a null Lebesgue measure. Hence we just give the topological gradient expression for a cracked domain (b) and develop the full proof for a perforated domain (a).

A Bayesian approach
In this section we present, according to an a priori modeling of the image, a statistical analysis to deduce the suitable variational model for restoring the observed noisy image . We denote by N the number of pixels in the support of the image Ω . The discrete domain is denoted Ω N . We set by u N (respectively by f N ) the discrete version of the image u to recover (respectively of the observed image f ). For each pixels p x ∈ Ω N , f N (p x ) and u N (p x ) can be viewed as a realization of the random variables U N (p x ) and F N (p x ) where U N and F N stands for the random vector formed by these variables at each pixel. We suppose that they are identically distributed and independent. The reasoning is as follows : we express the a priori density probability g U N |F N and then we apply the classical reasoning which consists of finding u N as the value maximizing this density probability called a maximum a posteriori estimator (MAP estimator). A discrete model associated to a discrete energy is deduced and then passing to the limit when N → ∞ we get the continuous variational model. Let g U N |F N be the a posteriori density probability that we want to maximize with respect to u N . Thanks to the Bayes rule, g U N |F N expresses as : on the noise model and g U N is an a priori density probability. Writing that u n is a minimum of −log(g U N |F N ) we get The a priori density g U N has to be determined, it will play the role of regularizing term . In analogy to statistical mechanics, a priori densities are frequently Gibbs functions [16] of the form : is a discrete version of a non negative energy functional R(u). The choice of the density probability g F N |U N depends on the statistic of the model to be considered. Below we review respectively the Gaussian model, the speckle model and finally the Poisson model.

Gaussian model
A classical modeling of image formations is : F N = U N + G N where U N is the disrete version of the image to recover and G N a Gaussian noise of mean 0 and of standard deviation σ . The density of the Gaussian noise is g G N To simplify we still denote by F N , S N and U N the random variables F N (p x ), S N (p x ) and U N (p x ). Let us express the conditionnal probability density g F N |U N . From the definition of the conditionnal probability we have : The conditionnal probability density g F N |U N ( f |u) is a function of the variable f and depending on a parameter u. From the model F N = U N + G N , the independency of U N and G N and a change of variable we get Hence by identification with (4) we deduce that g F N |U N ( f |u) = g G N ( f −u). Thanks to the independency hypothesis, the density of F N |U N is the product with respect to each pixel p x of the probability density of F N (p x )|U N (p x ). So the energy given in (3) rewrites in this case as with C a constant non depending on u. The constant σ 2 can be neglected in the model because it can be scaled with the regularization parameter γ. By passing to the limit when N → +∞ we get the following continuous energy

Speckle model
For SAR images, the classical modeling of the image is (see [28]) 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 N k for 1 ≤ k ≤ L and for each observations we have is a random variable which follows a negative exponential law with mean 1 and with density g G (x) = e −x 1 {x≥0} . Then, the observed image F N is construct from this L observations as : with Γ (L) = (L − 1)! (the mean of S N is 1 and its variance 1 L ). Now we can express the density g F N |U N . To simplify we still denote by F N , S N and U N the random variables F N (p x ), S N (p x ) and U N (p x ) . We start from the definition of the conditional probability : where g F N |U N ( f |u) is a function of the variable f and depending on a parameter u. Then, from the model F N = S N U N and the independency of U N and S N we have : Thanks to the definition of the probability and by a change of variable we get : Then by identification with (5) we deduce that Thanks to the independency F N (p x ) and U N (p x ) , the density of the conditional variable F N |U N is the product with respect to the pixels p x of the density F N (p x )|u N (p x ). By taking the −log function we deduce that (3) rewrites in this case as 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 Speckle with Log of the image (Speckle-Log model) As we can see in (7), the energy associated with the speckle model is not convex. By taking the logarithm of the speckle model we get : the new model writes as G N = V N + T N . Now the problem is to recover V N from the observation G N . We assume that the random variable V N follows a Gibbs prior. Let us calculate the density function of T N : So the density of T N is g T N (r) = g S N (e r )e r = L L (L−1)! e L(r−e r ) . Concerning the conditional density g G N |V N we have : where C is a constant no depending on v. By scaling the model with the constant γ we get the continuous energy as N → ∞: which is now a convex function of v.The recovered image is then u = e v .

Poissonian model
This model is classical in astronomical and confocal microscopy images [14]. Poissonian observations originates from the stochastic nature of photons emission. We denote R j , for j ∈ I ind (Ω ), the domain of R 2 modeling the pixel j and such that Ω is the disjoint union of all the (R j ) j∈I ind (Ω ) . We assume that f is a step function constant on each R j and we still denote f N = f the observed image seen as the realization of the random vector F.
More precisely for j ∈ I ind (Ω ), f j is a realization of a Poisson statistic of mean and variance equal to is a discrete version of u N ∈ R N (may be a step function or a bi-linear interpolation). Thanks to the independence of F j and U N j , the conditional probability P(F = f |U N = u) is given by : and by applying the −log function, we have : where C is constant independent of u. We deduce that (3) rewrites in this case as The dependence of E N with respect to u comes from the definition of λ N . Passing to the limit we get the continuous energy : 3 Blurring modeling In most imaging applications the optical material, the motion of the camera or of the target introduce a blur on the observed image (see [24]). Generally spatially invariant blur is modeled as a positive convolution operator u → Ku with K1 = 0. We denote by K N the N ×N matrix associated to the discrete version of K on Ω N . From section 2 we deduce the following models adapted to each kind of noise and taking into account the blur : 1. Gaussian model : the observed image writes as F N = K N U N + G N and by the same reasoning of section 2 we get the following energy : 2. Speckle model : the observed image writes as F N = S N K N U N and the energy is 3. Speckle model with the Log of the image (Speckle-Log model). We recall that the model writes as The deblurring cannot be handled simultaneously with the denoising process. After the denoising step we must solve the problem V N = log K N U N where the unknown U N can be found by a least square formula : but we know that this problem is ill-posed, particularly when K contains small eigenvalues. For this reason the blurring problem is not handled for speckle noise by our method. In this case, if we only want to correctly restore a blurred and speckled image it is preferable to use (11). 4. Poissonian model : the observed image at pixel p x is a realization of a Poisson statistic of mean R px Ku N (x)dx, so the energy is In the sequel we give the topological gradient for the Gaussian and Poisson models with blur and for the Speckle-Log model without blur.

Gaussian noise with blurring
We consider problem (1) with the energy (10) : The main particularity of this model is that it is linear. We do not give the calculus of the topological gradient here because of the similarity with the case without blurring (see [3]). We just give the topological gradient expression and some experimental results. Note that this expression only needs the resolution of two problems: the direct and the adjoint problems (we will see in the next section why an adjoint problem is necessary). By following the notations used in [3], the direct and the adjoint problems u 0 and v 0 are given in this case by: and We can show (see [6] chapter 3) that problems (P 0 ) and (Q 0 ) are well posed in H 1 (Ω ) as soon as K1 = 0 and γ > 0.
The topological gradients at x 0 ∈ Ω for a perforated domain (a) Ω ε = Ω \{x 0 + εB} and for a cracked domain (b) Ω ε = Ω \{x 0 + εσ }, denoted respectively by I b (x 0 ) and I c (x 0 ) can be easily deduced from the case without blur and are given by the following Theorem.

Theorem 1
The topological gradients associated to problems (14) and (15) and to the cost function J ε (u) = Ω ε |∇u| 2 , for a perforated and a cracked domain are respectively : with u 0 and v 0 given by (14) and (15) and γ > 0.

Speckle multiplicative noise
We consider the variational problem (1) with the energy given in (8). More precisely we are going to study the minimization problem : where ψ(x, u) = u − g(x) + e −(u−g(x)) and g = log( f ) is the logarithm of the observed image. We assume that there is no blur i.e. K is the identity operator. To shorten notations we write sometimes ψ(u) instead of ψ(x, u). Remark 1 In [4] the authors propose a speckle denoising model using the total variation model with the data fidelity term associated to (7).
This section is organized as follows. First, we show that (1) admits a unique solution for a more general class of functions ψ (verifying Hypotheses 1) and we prove that the solution verifies some min/max principles. Then we apply the result of the general case to show that problem (17) admits a unique solution for the speckle model ψ(x, u) = u − g(x) + e −(u−g(x)) . Finally, we perform the topological gradient for a general function ψ verifying Hypotheses 1. In the following we denote by I Ω the energy I Ω (u) = γ 2 Ω |∇u| 2 + Ω ψ(x, u) and we recall that we always denote by J Ω (u) = Ω |∇u| 2 the cost function. (1) In this subsection we first establish the well-posedness of (1) for a general class of functions ψ(x, u) and then we check that the function ψ(x, u) associated to the speckle-Log model (17) matches these hypotheses. To simplify we suppose that γ = 1 and sometimes we write ψ(u) for ψ(x, u).

Lemma 1 Let ψ(x, u) a function verfying Hypotheses 1, then
Proof Existence : Let (u n ) a minimizing sequence. There exists a constant C 1 such that I Ω (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 . So, we deduce the following inegality Let v n = max(u n , a), and Ω − n = Ω ∩ {u n ≤ a}, we have v n ≥ a and so v n is still a minimizing sequence and v n ≥ a. Similarly by setting w n = min(v n , b), we deduce that w n ≤ b and that w n is still a minimizing sequence. Therefore we can suppose that any minimizing sequence u n verifies a ≤ u n ≤ b. It is easily seen that u n is bounded in H 1 (Ω ). Thus, up to a subsequence there exists u ∈ H 1 (Ω ) such that u n L 2 (Ω ) → u and u n ⇀ stands for the weak topology). By using the lower semi-continuity of J Ω (u) and Fatou's Lemma we get that u is a solution of (1). Moreover we have a ≤ u ≤ b a.e on Ω .
Uniqueness : From the existence, we can work on the set which is a convex set and that I Ω has a unique minimum in H 1 (Ω ). We apply below Lemma 1 to the speckle-Log model.
Let a = log(α) and b = log(β ) the following inequalities hold In the next subsection we compute the topological gradient for a perforated domain. We just give the result for a cracked domain.

Computation of the topological gradient for a perforated domain
Let u ε = u Ω ε be the solution of problem (1) replacing Ω by Ω ε = Ω \{x 0 + εB} and let J ε (u) = J Ω ε (u). In order to establish the topological expansion for a more general class of PDE we assume that ψ(x, u) verifies Hypotheses 1. By writing that DI Ω ε (u ε ).v = 0 for all v ∈ H 1 (Ω ε ), we obtain the following variational formulation : find Then by an integration by parts, we deduce the following Euler equations associated with : We introduce the following functional We denote by u 0 the solution of (18) for ε = 0. We now give the main result of this section. To simplify notations we perform the proof with γ = 1 (but we state the Theorem with γ > 0).

Theorem 2
We denote by I b the topological gradient for a perforated domain. Then I b associated to problem (19) with ψ(x, u) verifying Hypotheses 1 and with the cost function J ε (u) = Ω ε |∇u| 2 is with u 0 and v 0 given by (19) and (25) for ε = 0 and with γ > 0.
Proof The topological gradient is given by the leading term and (21) it is classical [2,27] to introduce an adjoint problem. Due to the non linearity of the direct problem, we first do a Taylor expansion at second order of F ε (u, v) with respect to u at point u 0 : The Euler equations associated with (24) are : Remark 2 -In the proof, we take γ = 1.
-The adjoint problem (Q ε ) is linear and we can notice that the strict convexity of u → ψ(x, u) is necessary to (25) is wellposed and thanks to Lemma 1, we have the following inegality We deduce from (21) and (23) that : with By using an integration by parts, the term F ε (u 0 , v ε ) expresses as : we introduce the following extension on B ε : For with u 0 = u 0 − u 0 (0) and The computation of J ε , E 1 , E 2 and E 4 needs to approximate and that w ε = εQ x ε + r ε where Q is given by (56) with g = −∇v 0 (0).n and r ε 1,Ω ε = O(ε 2 −log(ε)).
Proposition 2 Let I ε , J ε , K ε , E 1 , E 2 , E 3 and E 4 given by (22), (27) and (30), we have the following estimations : Proof The first equality is straightforward. Lemma 4, the regularity of ψ(x, u) stated in Hypotheses 1 and Proposition 1 permit to say that E 1 and E 2 are O( u ε − u 0 2 0,Ω ε ) and then ∼ O(−ε 4 log(ε)). We recall the notation ψ(u) means ψ(x, u(x)). A change of variable (CV) and the continuity of x → D u ψ(x, u 0 (x)) lead to Again with a CV, the equality ∆ v 0 = D 2 u ψ(u 0 )(v 0 − u 0 ) − D u ψ(u 0 ) and the regularity of ψ(u) we get E 3 ∼ O(ε 3 ). By using Lemma 5 and Lemma 2 (see Appendix A) we have : A CV and a Taylor expansion of u 0 and v 0 at 0 lead to F 1 = O(ε 3 ) and F 2 ≤ Cε 2 ∂ n l e ε ε (εX) −1/2,∂ B . For F 2 it suffices to make a CV and use the trace Theorem on B 2 \B: Now from the equivalency of the H 1 (B 2 \B)/R-norm and the semi norm and a CV we get ∂ n l e ε ε (εX) −1/2,∂ B ≤ C|e ε | 1,Ω ε . By using Lemma 5, we obtain F 2 = O(ε 3 −log(ε)). Finally by using a CV, the continuity of ϕ → l ϕ from H 1/2 (∂ B) to H 1 (B), the trace Theorem on B 2 \B, a CV again and Lemma 5 we have : The topological expression is easily deduced from Proposition 2 and Theorem 2 is proven.

Expression of the topological gradient for a cracked domain
For the cracked domain Ω ε = Ω \x 0 + εσ (n), calculus are similar. The term I ε of (26) is zero and the term The asymptotic expansion of u ε and v ε are similar and then the computation of the topological gradient is the same as in the linear case (see [3,15] for more details).

Theorem 3
We denote by I c the topological gradient for the cracked domain Ω ε = Ω \{x 0 + εσ (n)}. The topological gradient associated to problem (19) and to the cost function J ε (u) = Ω ε |∇u| 2 for a cracked domain is with u 0 and v 0 given by (19) and (25) for ε = 0 and with γ > 0 in front of the laplacian.

Poissonian model with blurring
We recall that we model the observed image f as a step function defined by f j ∈ R on R j for j ∈ I ind (Ω ), where R j is a regular domain of R 2 modeling the pixel j and we denote N = |I ind (Ω )|. To simplify, we suppose that |R j | = 1 and that Ω is the disjoint union of the (R i ) i∈I ind (Ω ) . We assume that min j∈I ind (Ω ) f j > 0. We recall the general minimization problem associated to the Poisson model given in (12) : . First we show that problem (33) is well-posed, then we compute the topological gradient for a perforated domain (a) : Ω ε = Ω \x 0 + εB, and we give the expression for a cracked domain (b) : Ω ε = Ω \x 0 + εσ without proof.

Well-posedness of problem (33)
Proposition 3 Let f a step function such that min i f i > 0 and max i f i < +∞, then problem (33) with ψ j (x) = X − f j log(x) for j ∈ I ind (Ω ) admits a unique solution u ∈ H 1 (Ω ).

Besides this solution is such that
Proof Existence : To simplify the proof we suppose that K is the identity operator and γ = 1. The proof for the general case is quite similar. For more details see chapter 3 of [6]. We must add a priori to (33) the condition R i u > 0, ∀i ∈ I ind (Ω ). We set H = u ∈ H 1 (Ω ), R i u > 0 ∀i ∈ I ind (Ω ) . Then (33) rewrites as : By using the positiveness of Ω |∇u n | 2 we deduce that ∑ j ψ j R j u n ≤ D. By setting K i = ∑ j =i min x ψ j , it is straightforward that ψ i R i u n ≤ D − K i and then Hence the constraint R i u n > 0 is fulfilled. We deduce that ∑ i E i ≤ Ω u n = ∑ i R i u n ≤ ∑ i E i and thanks to Poincaré-Wirtinger Lemma we get that 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 J Ω (u) ≤ lim inf J Ω (u n ) and thanks to (34) and Bolzano-Weierstrass Lemma we can extract a subsequence u n such that R i u n → l i ∈ R and by continuity ψ i R i u n → ψ i (l i ) ∀i ∈ I ind (Ω ).

Finally we have
and so u is a minimizer of I Ω (u) in H 1 (Ω ).
with Dψ j R j u = 1 − f j R j u .
-By taking as function test v = 1 we get the egality N = ∑ j f j R j u . As Uniqueness : From the two previous points we can consider the minimization space for all j ∈ I ind (Ω ), by linearity of the integral we deduce that u → ψ j R j u is strictly convex on H(Ω ). As J Ω is strictly convex we deduce that I Ω (u) is strictly convex on H(Ω ) which is a convex space and that I Ω has a unique minimum in H 1 (Ω ).
Remark 3 -Under the same hypotheses on f, we get the existence and uniqueness of a solution u ε in H 1 (Ω ε ) for (19). For ε small enough, we still have -When K = I is such that K1 = 0 we can show that problem (33) is well-posed in H 1 (Ω ). For more details we refer the reader to [6] chapter 3. -We can show that Proposition 3 holds as soon as ψ j are bounded from below for j ∈ I ind (Ω ) and strictly convex on I ⊂ R. In the general case α and β are implicitly defined in function of the ψ j .

Computation of the topological gradient for a perforated domain
In this section we compute the topological gradient for a perforated domain Ω ε = Ω \{x 0 + εB}.
is the ball centered at x 0 and of radius ε. For j ∈ I ind (Ω ), let R ε j be the domain equal to R j 0 \B ε (x 0 ) if j = j 0 and R j otherwise. Now let us denote I ε j (u) = R ε j u, I j (u) = R j u, J ε (u) = J Ω ε (u) and u ε the solution of (33) replacing Ω by Ω ε .
By writing that DI Ω ε .v = 0 ∀v ∈ H 1 (Ω ε ), we deduce the following variational formulation of (33) : By taking v ∈ D(R ε j ) the space of C ∞ (R ε j ) functions with compact support in R ε j , we get −γ∆ u ε +Dψ j (I ε j (Ku ε ))K ⋆ 1 = 0 on R ε j , ∀ j ∈ I ind (Ω ). Then if v is any test function v ∈ H 1 (Ω ), we deduce ∂ n u ε = 0 on ∂ Ω ε therefore we have : where [u ε ] ∂ R j denotes the jump of u ε through ∂ R j . We now to give the main result of this section. The proof is performed in the case γ = 1 and K = I but the proof is the same in the general case.

Theorem 4
The topological gradient I b associated to problem (37) and to the cost function J ε (u) = Ω ε |∇u| 2 for a perforated domain is with u 0 and v 0 given by (37) and (43) for ε = 0.
Proof We set Then in order to introduce the adjoint problem, we make a second order Taylor expansion with respect to u for F ε (u, v): By taking u ∈ D(R ε j ) and by integrating by parts, we deduce the following Euler equations associated with (42) Remark 4 For γ = 1 and K = I the adjoint problem will be defined by From (39) and (41) we deduce that with By using an integration by parts, G ε (u 0 , v ε ) expresses as With a similar calculus as the one made in (29) for the speckle model, G ε (u 0 , v ε ) rewrites as : In the next section we give the asymptotic expansion of the previous quantities. (45) and (46), then we have the following estimations : 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 7 (see Appendix B) and the fact that R ε j 0 ε→0 −→ R j 0 give the second estimation. For K ε we refer the reader to the proof of Proposition 2. For E 1 , we use Lemma 6, the regularity of ψ(x) and that ξ ε ∈ [α, β ] : By using that ∆ v 0 = D 2 ψ j 0 (I j 0 (u 0 ))I j 0 (v 0 ) and a Taylor expansion of u 0 at 0 we get E 2 = O(ε 3 ). For E 3 , by a change of variable and Lemma 7, we get Denoting by j 0 the integer such that R j 0 ∋ x 0 , we deduce the expression given in Theorem 4.

Expressions of the topological gradient for a cracked domain
For a cracked domain Ω ε = Ω \{x 0 + εσ (n)}, the calculus are similar. The term I ε of (44) is zero and the term F ε (u 0 , v ε ) expresses as F ε (u 0 , v ε ) = − σ ε ∂ n u 0 [w ε ]. The topological expansion of u ε and v ε are similar with the perforated domain and the computation of the topological gradient is the same as the linear case (see [3,15] for more details).

Theorem 5
The topological gradient I c associated to problem (37) and to the cost function J ε (u) = Ω ε |∇u| 2 for a cracked domain is with u 0 and v 0 given by (37) and (43) for ε = 0.

Remark 5
The topological gradient is the same in the general case of functions ψ j ∈ C 3 (I), strictly convex on I and bounded from below on I. Just in the right hand-side of (43), K ⋆ 1 must be replaced by D 2 ψ j (I j (Ku 0 ))I j (Ku 0 ) + Dψ j (I j (Ku 0 )) K ⋆ 1.

Restoration using the topological gradient for a cracked domain
As a by product the calculus of the topological gradient for a cracked domain (TGC) allows to restore images degraded by blur or/and various noise statistics. Once the TGC is computed we define for a fixed threshold δ > 0, the set E δ = {x ∈ Ω ; |T GC(x)| ≥ δ } and the characteristic function where η is a small positive parameter. From the computation of the TGC we also get the normalized direction τ = n ⊥ of the edge. If n = (cos(ϕ), sin(ϕ)) is the normal to the crack given by the TGC, we have τ = (sin(ϕ), −cos(ϕ)). Then, if f is the degraded observed image, we want to find a restored version u 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 : 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.
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 (49) are very good.

Numerical application to 2D imaging
In this section we illustrate the theory of topological gradient by giving various experimental results for models (13), (17) and (33).
The topological gradient expressions for the three models are stated in section 4, 5 and 6.
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. If the crack model is used, the topological expression is the same for these three models. First, since equations (14) and (15)  develop the numerical analysis in a specific subsection. Then we perform the discretization of problems (17) and (33) and finally we give the experimental results. As the adjoint problems (25) and (43) 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.

Numerical analysis for Gaussian model with blurring
To discretize (14) and (15) we use a DCT1 (discrete cosine transform of type 1) because of the symmetry properties and the fact that the DCT1 of a convolution product of two vectors is the product of the DCT1 of each vector. We choose this discretization because of the symmetry properties guaranteed by the algorithm. A DCT1 of N points is equivalent to a DFT (discrete Fourier transform) of 2N-2 points.
We use the FFT (fast Fourier transform) to perform the DFT. The computation time of a FFT is a O (Nlog(N)). A numerical study shows that for N ≤ 10 10 , FFT is faster than a finite difference scheme. Algorithm 2 gives the different steps to compute the discrete solutions (14) and (15). It consists in the following steps : -Symmetric extension of the initial N y × N x image in an 2(N y − 1) × 2(N x − 1) image and extension of the 2n y + 1 × 2n x + 1 kernel in a 2(N y − 1) × 2(N x − 1) kernel. To fix ideas in 1D and for n x = 2, the extension of the dis- }, but it is not a good choice. This fact is explained in [19]. Let us give the reasoning in 1D for a 1-periodic function. The trigonometric function associated with the vector of DFT coefficients is where y N ∈ R N is the vector such that y N = DFT (y N ), y N = u N (l/N) 0≤l≤N−1 and m k ∈ Z are coefficients which does not change the function u at points x l = l N , but it greatly modifies u between these points (aliasing appears). If we compute the L 2 (0, 1)-norm of the first derivative we get that From (51), we see that m k changes considerably the L 2 (0, 1)norm of the first derivative and the good choice for m k is the value minimizing (k + m k N) 2 . If 0 ≤ k < N/2 then (k + m k N) 2 is minimized for m k = 0 and if N/2 ≤ k < N then (k + m k N) 2 is minimized for m k = −1. A special consideration is made for k = N/2 when N is even because of the two possible choices (m k = −1 or m k = 0, see [19] for more details). By following these considerations we define the fre- − 2), ..., −1} and E y = 0, ..., N y − 1, −(N y − 2), ..., −1 . We denote by Λ = (Λ x ,Λ y ) the 2(N y − 1) × 2(N x − 1) the mesh grids associated with this discrete space and the vector of Fourier coefficients associated to a discrete signal x ∈ R N is denoted by X.

Numerical analysis for Poisson and Speckle models
If in f Ω f > 0, by Proposition 1 and Proposition 3, problems (17) and (33) are well-posed and can be discretized as : . f N is a discretization of f ; J p (x) and J s (x) are respectively the discrete versions of the energy functions (12) and (8). We choose a simple discretization : u N (x) is the step function equals to u( j) on pixel j, and we represent u N by a vector of R N .

Given a blurring kernel of convolution k
, use the procedure described in section 8.1 to calculate its symmetric extension k mn for 0 ≤ m < 2(Ny − 1) and 0 ≤ n < 2(N x − 1).

Use an FFT to compute
4. Given Λ = (Λ x ,Λ y ) the meshgrid associated to the frequencies domain described in section 8.1, compute

Use an inverse FFT to compute u i j and
During the construction of the sequence x (k) , the condition x (k) ≥ α p for the Poisson model and x (k) ≥ log(α s ) must be fulfilled at each step. Hence a projection ensures this condition. To solve these problems we use an iterative algorithm based on the descent method called the SGP algorithm [12] (scaled gradient projection). 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 by assuming 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 p and J s are denoted by J as soon as we do not use their expression and by δ the number equal to α for the Poisson model and equal to log(α) for the Speckle-Log model. We set by Λ = x ∈ R N , x ≥ δ . We want to find x ⋆ ∈ Λ such that ∇J(x ⋆ ) = 0. At step k, a first order Taylor expansion at point x = x (k) leads to the following equation We deduce by this reasoning that the direction of the descent algorithm can be given by ∇ 2 J(x (k) ) −1 ∇J(x (k) ), but we see that the computation of this direction is very costly. We denote by D L the compact set of the 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(x (k) ) and to project each iterate on Λ with respect to the norm x D = √ x T Dx. We set P Λ ,D −1 for D ∈ D L the projector on Λ with respect to the norm . D .

5:
if y (k) = x (k) then 6: Stop, x (k) is a stationary point. 7: end if 8: Descent direction d (k) = y (k) − x (k) ; 9: λ k ← 1 and J max ← max 0≤ j≤min(k,M−1) J(x (k− j) ) 10: λ k fixed by backtracking : end while 14: x The construction of the sequences D k and α k needs some explanations. We choose The approximation of the Hessian matrix ∇ 2 J(x (k) ) is B(α k ) = α k D k . By using a first order Taylor expansion of ∇J(x) at point x (k−1) we get that Hence there two possible choices of α k can be made : In [12] 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 . Let us give the derivative of the discrete cost functions J p and J s : where 1 ∈ R N denotes the vector with each coefficient equal to 1, diag(x) for x ∈ R N is the diagonal matrix with the diagonal equal to x and for x ∈ R N and ϕ : R −→ R a function, ϕ(x) stands for the vector (ϕ(x i )) 1≤i≤N . The choice of the parameters in Algorithm 3 is the following : β = 10 −4 , θ = 0.4, k max = 600, M = 1 and for the Poisson model (33) we set α min = 10 −10 , α max = 10 5 while for the Speckle-Log model (17) we set α min = 10 −5 , α max = 10 15 . The initial value of x (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 the speckle model, (49) is performed with ψ(u) = log(u) + f u and this problem is discretized for both the speckle and the Poisson models by the SGP algorithm.

Comparison of our method with some classical models
As said in the introduction other variational methods exist for segmenting/restorating images.
We will compare the topological gradient segmentation process with the one performed by the Mumford-Shah model. We will also compare the restoration proposed in (49) with the ones given by the Mumford-Shah restoration and by the TV restoration.

Mumford-Shah model of segmentation/restoration and its approximation
Let u the image of support Ω and K ⊂ Ω , the functional introduced by Mumford and Shah in 1989 (see [23]) 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 was 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 (Gaussian, Poissonian and speckle model)i.e. the model that we will compare with (49) is (see [25]) : where K is the blur operator, u(x) is the restored image, 1 − b(x) ≈ 0 is the characteristic function of the edges and ψ(u) depends on the model as follows : We will call this model the Mumford-Shah model.

TV model of restoration
In most recent variational models we search for a restored version minimizing an energy functional composed of the total variation of the function and a data fidelity term which depends on the a priori model (see [4] for the speckle model).
The model studied for comparisons with model (49) is the following : where λ is a parameter, K the blur operator, and ψ(u) is given in (54). In the sequel we call this model the TV model. For more details on these models we refer the reader to [4], [25] and [23].

Numerical results for the Gaussian model
An interesting study is the comparison of formula (16a) and (16b) giving the topological gradient for respectively a perforated and a cracked domain. A priori formula (16a) (TGB), associated to a perforated domain, would be adapted for the detection of isotropic small structures (about 10 pixels of aera) and formula (16b) (TGC), associated to the cracked domain, for the detection of edges of big objects (more than 20 pixels of aera). Fig. 1 displays the TGC and the TGB for different values of γ compared to the function b(x) performed by minimizing the Mumford-Shah functional (53). For small γ the TGC seems more robust than the TGB. For γ = 1 the TGC better detects the black spots contours of the cheetah compared to those given by the Mumford-Shah model. We remark by increasing γ that the TGB becomes singular on the black spots contours while the TGC better detects the border of the cheetah.
We deduce that γ must be tuned with respect to the noise but also by taking into account the size of structures to detect. Moreover we must notice that if the TGB seems to be less robust with respect to noise than the TGC, it is easier and faster to compute the TGB than the TGC which is performed by minimizing an expression over the normal n of the crack.    On Fig. 5, we compare the TGC (32) and the TGB (20) for different values of γ with the function b(x) given by the Mumford-Shah model for a synthetic speckled image. Similarly to the Gaussian case, the TGC seems better compared to the TGB and the Mumford-Shah segmentation. We still deduce that γ must be tuned with respect to the noise and to the size of structures to detect. The result given in Fig. 6 for a real SAR image is comparable to the one of Fig. 5. Here we see that the TGB and the TGC can be used for different objectives : particularly on small structures we see that the TGB detects the entire object while the TGC detects its edges.

Numerical results for the Poisson model. Comparisons
In this section we compare the segmentation performed by the TGC (48), the TGB (38) and the Mumford-Shah model (53). We also compare the restoration computed with (49) the Mumford-Shah model and the TV model. Fig. 11 and Fig. 12 show respectively the segmentation results in the case of a synthetic Poissonian image and of a real confocal image of a rat's neuron. TGC (48) detects edges quite well compared to the Mumford-Shah model. We see that TGB fills small structures (the size of these structures is related to γ).  Finally, Fig. 16 shows the 1D profiles of the image to recover, its degraded versions (blurred, blurred + Poissonian process), the restored version (49) and the TGC (48) across an edge. We see that (49) allows to recover the initial image and that the TGC detects very well the edge.

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 [3]. To establish these asymptotic expansions we need the following exterior problem where g ∈ H −1/2 (∂ B) and ∂ B gdσ = 0. For the computation of the topological gradient we wiil need the two following usefull lemma. We omit the proofs and we refer the reader to [15] for more details The following asymptotic estimations holds.

Lemma 3
Let H the solution of (56), then :