Image Restoration via Topological Derivative

The problem of image restoration has deserved considerable attention in resent years. For the visual analysis of images, clarity and visibility of details are important factors, but for advanced processing, a high signal-to-noise ratio (SNR) is essential, as further processing steps (such as segmentation and classification) are sensitive to noise. Though the years, different techniques have been studied to improve the SNR or a degraded image. Techniques based on post-processing have the advantage of not affecting the acquisition process (Gonzalez & Woods (2001); Jain (1989); Weickert (1995)). More recently, the work of Tschumperle (2006) has explored the more extensive use of curve-preserving PDE’s for restoration of images. The calculation of mean intensities over neighboring pixels, equivalent to isotropic diffusion, considerably increases the SNR, but degrades the quality of image features (edges, lines and dots). This effect can be reduced with non-linear filters. The median filter has the characteristic of maintaining these features, but details are lost, degrading the image resolution. Perhaps the most popular technique introduced in the last couple of years is anisotropic diffusion, initially proposed by Perona and Malik (Black et al. (1998); Perona & Malik (1990)).


Introduction
The problem of image restoration has deserved considerable attention in resent years.For the visual analysis of images, clarity and visibility of details are important factors, but for advanced processing, a high signal-to-noise ratio (SNR) is essential, as further processing steps (such as segmentation and classification) are sensitive to noise.Though the years, different techniques have been studied to improve the SNR or a degraded image.Techniques based on post-processing have the advantage of not affecting the acquisition process (Gonzalez & Woods (2001); Jain (1989); Weickert (1995)).More recently, the work of Tschumperle (2006) has explored the more extensive use of curve-preserving PDE's for restoration of images.The calculation of mean intensities over neighboring pixels, equivalent to isotropic diffusion, considerably increases the SNR, but degrades the quality of image features (edges, lines and dots).This effect can be reduced with non-linear filters.The median filter has the characteristic of maintaining these features, but details are lost, degrading the image resolution.Perhaps the most popular technique introduced in the last couple of years is anisotropic diffusion, initially proposed by Perona and Malik (Black et al. (1998); Perona & Malik (1990)).This problem has motivated interdisciplinary research and the use of techniques actually born in other areas, as is the case of the Topological Derivative (TD).The TD has been originally conceived for the study of topology optimization and property identification problems.Since 1994 different works proposed this new paradigm for the study of such problems.The pioneer works of Eschenauer in 1994 (Eschenauer et al. (1994)) and Schumacher in 1995 (Schumacher (1995)) introduced a way to obtain the optimal shape and topology using Topological Sensitivity Analysis.In summary, this new concept called asymptotical topological expansion is posed as follows.Let J (Ω)=F (u(Ω)) be an arbitrary cost function that measures the "quality" associated to a given topology characterized by the state function u(Ω), which is restricted to be the solution of a variational equation defined in Ω.Given a sufficiently small positive number ǫ, a positive function f (ǫ) that goes to zero with ǫ, we call Ω ǫ the perturbed domain after the inclusion at x of a hole of infinitesimal size governed by ǫ.Therefore, the asymptotical topological expansion J (Ω ǫ )=J (Ω)+ f (ǫ)D T ( x)+O( f (ǫ)) (1) provides an estimate of the cost function value in the perturbed domain for a sufficiently small ǫ, where D T is known as the TD.The Topological Derivative can then be defined as (Novotny (2003)) A scalar function, defined in Ω, indicating in each point x ∈ Ω the sensitivity of a cost function when a hole of infinitesimal size ǫ is introduced at that point.

Motivation for the use of Topological Derivative in image restoration
In 2005 appeared the first papers using the TD in image processing: restoration by Belaid et al. (2007) and Larrabide et al. (Larrabide, Novotny, Feijóo & Taroco (2005;2006)), where the objective was to recover an image that suffered some kind of degradation; segmentation by Larrabide et al. (Larrabide, Feijóo, Novotny, Taroco & Masmoudi (2005); Larrabide, Feijóo, Taroco & Novotny (2006)) and Hintermüler (2005), in medical images where the interest is in the identification of different organs for posterior reconstruction; and image classification by Auroux et al. (2006).He & Osher (2006) established a relation between the TD and other techniques broadly used in image processing as is the case of level sets.A remarkable feature of the TD is that it allows to compute the variation of a cost function with respect to a parameter that changes non-smoothly (e.g., characteristic function of the domain, material properties or non-continuous change of the forces acting on the problem).This derivative can be used to identify, according to some criterion, the characteristic function of an optimal domain, the material properties and their distribution in a given domain, or even the forces acting on a given domain and how they are distributed.This problem appears frequently in the context of image processing.As examples we can mention identification of edges, identification of objects, object tracking, decomposition of texture and geometry and reconstruction from projections, where the use of the TD appears as a natural way to solve these problems.
Regarding image restoration, the stationary heat equation has been used for this purpose (Kornprobst et al. (1997)).In this approach, the diffusion coefficient is usually given by The problem consists in determining the function φ that allowed to remove noise preserving edges in the image.One property that characterized restoration methods based on non-linear isotropic diffusion was removing noise along the edges in the image.This unwanted property of non-linear diffusion, but still isotropic, was partially reduced when a non-linear anisotropic diffusion tensor was introduced (Frangakis & Hegerl (2001)).In this case, a diffusivity tensor c(|∇u|) was constructed from two eigenvectors and eigenvalues of tensor J = ∇u ⊗∇u.Still, as diffusion across edges is not completely stopped, heuristic criterion must be introduced to avoid the loss of image details.
Stoping the diffusion in the direction orthogonal to the intensity iso-lines of u is somehow equivalent to introducing a crack in the same direction.In this way, the use of the TD in image restoration appears naturally as it allows to analyze abrupt variations in the material

98
Image Restoration -Recent Advances and Applications www.intechopen.comproperties.For instance, the TD can be used to determine the diffusion tensor, which in the following we recall as K.To formalize this, we consider the image characterized by its intensity u 0 ∈ L 2 (Ω) defined in a limited open domain Ω ∈ R 2 (the extension to R 3 is straightforward).For each restored image, characterized by the intensity u ∈ H 1 (Ω),w e can associate the cost function measuring the "quality" of the restoration given by the solution of the following variational problem: Determine u ∈ H 1 (Ω) such that Different methods exist to remove noise and enhance edges of a degraded image.We can distinguish two types: based on the solution of a stationary problem and based on the solution of an evolutionary problem.Both types of methods are based on the application of non-linear/anisotropic diffusion on the image.In both cases, the diffusion coefficient or tensor is computed as a function of the local image gradient.This coefficient takes small values for high gradients (edges) thus stoping diffusion, and higher values when the gradient is small (homogeneous regions) promoting a higher diffusion.Both methods have two parameters.
The stationary method has a parameter determining which gradients will be considered as edges and which ones not, and a second one characterizing the intensity of the diffusion to be applied.In the case of evolutionary methods, the first parameter is in some way similar to the one used by the stationary method to determine the threshold for gradients that are considered edges, and the number of iterations.For both methods, the parameter determining the gradient threshold can be estimated.But this does not happen for the other parameters, which need to be set depending on the noise type and intensity.In both cases, the the selection of this parameter will determine the quality of the result.

Topological Derivative
Topological sensitivity analysis allows to characterize the sensitivity of a problem when the domain Ω where the problem is defined is perturbed in some way.This perturbation can be a: • change in topology: the domain Ω is perturbed introducing a hole of an arbitrary shape ω ǫ at the point x ∈ Ω, and the TD provides the sensitivity of a cost function when ǫ → 0 (see Eq. ( 1)).• change in material properties: a perturbation in the material properties at point x ∈ Ω of an arbitrary shape ω ǫ is introduced and the TD provides the sensitivity of the cost function when ǫ → 0 (see Eq. ( 1)).By "material properties" is meant the coefficients that define the variational problem associated to the problem.• change in the forces/sources acting on Ω: similar to the previous case, but perturbing the forces/sources.
In the following, and for the sake of simplicity, only the first case will be considered, namely the perturbation of the domain by the introduction of a hole.The extension to the other two cases is straightforward.
where U characterizes a set (usually a linear manifold of V) of admissible functions defined in Ω e V = V (Ω) the vector space of admissible variations.Also, a(., .): U×V → R is a symmetrical bilinear form and l(.) : V → R a linear form.These forms also satisfy the properties of continuity and coercivity to warrant existence and uniqueness of solution u.
Let ω be a open domain arbitrarily shaped and of regular boundary ∂ω containing the origin.Given ǫ > 0 sufficiently small it can be defined for any point x ∈ Ω the domain ω ǫ given by ω ǫ = x + ǫω.In this way, the introduction of a hole ω ǫ centered in x ∈ Ω allows to characterize the perturbed domain Ω ǫ (Fig. 1) given by From Eq. ( 1), D T in x ∈ Ω can be defined as where , being u ǫ the solution of the same state equation now defined in the perturbed domain, namely in In the work of Novotny et al. (2003) a relation between the TD and classical shape sensitivity analysis (Haug et al. (1986); Murat & Simon (1976)) is established.This result permits to use tools developed in classical sensitivity analysis for the computation of the TD.This new approach can be stated in the form of theorem as: Theorem 1.Let f (ǫ) be a function chosen such that 0 < |D T ( x)| < ∞, then the limit when ǫ → 0 appearing in (4), can be written as where dJ (Ω τ ) dτ is the classical shape sensitivity.
Proof 1.The proof of this theorem can be found in the work of Novotny et al. (2003).
In the previous expression is implicit the domain transformation (deformation) χ τ : where v is the velocity field characterizing the shape change and For further information on this type of transformation see the work of Haug et al. (1986) and Haug & Céa (1981), Pironneau (1984), Sokolowski & Zolésio (1992) and Zolésio (1981).
From this theorem is naturally deduced the Topological-Shape Sensitivity Method, which will be explored in the following.The shape change derivative of the cost function in relation to the parameter τ can be written as where and where the notation J τ (u τ ) evidences the dependency of the cost function on u τ and on Ω τ .
To compute the derivative to change in shape considering the state equation as a restriction, the Lagrangian method is used (i.e., relaxing the restriction by the introduction of a Lagrange multiplier).The Lagrangian of this problem is written as We verify, for v = u τ , that We compute the derivative with respect to τ in Eq. ( 10), then We then work term by term on Eq. ( 12).Starting by the third term on the right-hand side we have 101 Image Restoration via Topological Derivative www.intechopen.com Then, for the particular case v = u τ , Eq. ( 13) is zero.Considering the second term of Eq. ( 12), we have where the symmetry of a τ (•, •) was used.In this expression, η can be chosen arbitrarily.In articular, it will be selected η = q τ , being q τ ∈V τ the solution to the adjoint equation given by The previous equation is known as "(variational) adjoint equation", and its solution q τ (or q ǫ and q if the adjoint equation is defined in the domain Ω ǫ and Ω respectively) as "adjoint solution".We note that, because of the properties of a(•, •), the adjoint equation is of the same kind as the state equation (Eq.( 2), or its counterpart in the perturbed domain Eq. ( 5)).From the computational point of view, the former means that the same computational system used to compute the solution of the state equation u (or u ǫ ) can be used to compute q (or q ǫ ).
The total derivative with respect to parameter τ of the Lagrangian is given by We notice that u τ | τ=0 = u ǫ and q τ | τ=0 = q ǫ .Therefore, the former expression results in a function of u ǫ and q ǫ and its derivatives.As we noted before, only the boundary ∂ω ǫ is perturbed by a uniform expansion (Eq.( 8)).Then, the derivative to shape change results in an integral only defined on the boundary ∂ω ǫ .Therefore, the topological derivative is given by an expression of the form where Σ ǫ depends on u ǫ and q ǫ , and it can be interpreted as a generalization of the Eshelby momentum energy tensor proposed by Eshelby (1985).The tensor Σ ǫ must be identified for each particular problem, which depends on the cost function adopted and the state equation associated to u ǫ .
Finally, its necessary to compute the limit when ǫ → 0 to obtain the D T expression.For this, it is necessary to study the behavior of solutions u ǫ and q ǫ when ǫ → 0. This behavior can be obtained with asymptotical analysis on the neighborhood of the hole.At this point, different alternatives can be used depending on the problem under consideration (boundary conditions, type of perturbation, etc.).In all these cases, asymptotical analysis allows to express u ǫ and q ǫ as a function of ǫ, u( x), q( x) and there derivatives in x, respectively.Namely, as a function of the solutions of the state and adjoint equation defined in the domain without perturbation providing as well the function f (ǫ).In this way, substituting these equation in Eq. ( 17) and from the computation of the limit for ǫ → 0 the final expression of the D T at point 102 Image Restoration -Recent Advances and Applications www.intechopen.com x is obtained.As mentioned, this expression will only depend on the values of u and q and its derivatives at point x.The former has consequences from the computational point of view and in fact, once u and q are obtained, the computation of the TD is just a post processing.
Then, for a given cost function F (Ω, u(Ω)) the Topological-Shape Sensitivity Method can be summarized in the following steps: 1. Compute the shape change derivative for the cost function F (Ω ǫ , u ǫ (Ω ǫ )), using the Lagrangian method.2. Identify the Eshelby momentum tensor Σ ǫ and write the sensitivity expression as an integral defined on the boundary ∂ω ǫ .3. Do the asymptotical analysis to study the behavior of u ǫ and q ǫ when ǫ → 0. 4. From the asymptotical analysis, chose f (ǫ). 5. Compute the D T using Eq. ( 17).This is a general and systematic way to compute the D T for an arbitrary cost function.The particular case presented in this chapter is applied to the introduction of a perturbation in the domain Ω in the shape of a crack, which is then presented as part of an algorithm for the restoration of degraded images.

Topological Derivative in image restoration
The concept of the TD, allows to quantify the sensitivity of a performance measure of cost function when the problem definition domain is perturbed.Therefore, if the cost function is associated to the noise present in the image, it will be possible to use its TD to develop appropriate image restoration algorithms.In this context, two state equations are studied: the first one based on a stationary problem and the second one on the evolutionary problem.The TD will allow to determine the location and orientation of the cracks that should be introduced in the domain to minimize the cost function.From the point of view of the state equation, the cracks will stop diffusion in the orthogonal direction, only allowing diffusion in the tangent direction.In other words, the TD provides a procedure to compute the diffusion anisotropic tensorial field that will be used to restore the image preserving the image edges.

Continuous approach -R D T -continuous
This approach is based in introducing cracks in an image under the effect of diffusion.These cracks are introduced at specific locations in the image using the information provided by the TD.For a small ǫ, let Ω ǫ = Ω \ γ ǫ be the domain perturbed by the insertion of a crack γ ǫ = x + ǫγ, where x ∈ Ω and γ(m) is a straight crack, being m the normal direction of γ (Fig. 2).Let u ǫ be the solution of the same variational problem in the perturbed domain Ω ǫ and J (u ǫ ) a cost function associated to it.Then, we obtain the following asymptotic topological expansion of J ǫ (u ǫ ) when ǫ → 0, i.e.,  ( 18) and q the solution of the adjoint equation For any point x, D T ( x, m) reaches its minimum when m is an eigenvector associated to the smallest eigenvalue κ min of M.Then, it is considered as the optimal direction of the crack γ ǫ (m) the eigenvector corresponding to the eigenvalue κ min .This value will be adopted as the TD associated to the creation of a crack at the point x.In Fig. 3 is presented an example for the Lena image (SNR = 26dB).
As mentioned, for any x , D T ( x, m) takes the minimal value when m is the eigenvector associated to the smallest eigenvalue κ min of M.Then, by considering the orientation of the crack γ ǫ (m) the eigenvector corresponding to the eigenvalue κ min .This minimal value will be adopted as the TD associated to the creation of a crack at the point x.The algorithm proposed here consists in computing the TD and introducing small cracks in the locations where the derivative is smaller than a given value D T Lim .Two algorithms are presented: isotropic and anisotropic.
To solve the numerical problem, the introduction of a small diffusion coefficient (or a conductivity tensor that acts on one direction) is interpreted as the presence of a crack.In the proposed algorithms, the tensor K(x) for the isotropic and anisotropic case are computed as a function of the TD, namely K = K(x, D T ): • Isotropic diffusion based on the TD (R D T -Continuous (Iso)):   for k ǫ ≪ 1, k 0 a positive real number, n and t the normal ad tangent directions to the crack, respectively.
With this diffusion tensor the restored image is obtained by solving the following variational problem: Determine As the solution u * of the variational problem given by Eq. ( 19) cannot be explicitly known in general, its necessary to compute an approximate solution.The Finite Element Method is used for this purpose Hughes (2000).Then, using the simplest finite element given by quadrilateral bilinear element (for the 2 dimensional case) or a trilinear parallelepiped (for the 3 dimensional case) whose nodal points coincide with the centers of the image elements,  approximate solutions u h of u, q h of q and u * h of u * can be easily obtained for any image u 0 ∈ L 2 (Ω).Using these solution, an approximation by finite elements M h of the tensor M is given by To find the restored image, three boundary value problems need to be computed.These correspond to the scalar fields of u h the isotropic problem, q h the adjoint problem and u * h the problem with the diffusivity tensor K(x) computed using R D T -Continuous (Iso or Aniso).Ensure: Restored image u * ∈ H 1 (Ω).compute u and q, solutions of the state and adjoint equation, respectively, compute the matrix 2 × 2 M and its minimal eigenvalue κ min for each point in Ω, find K using R D T -Continuous (Iso or Aniso), compute u * , a restored image, using K(x) previously obtained.
This algorithm uses one parameter (D T Lim ) to select the elements in the image that will have their coefficients with a modified diffusivity.Then, depending on the TD value of the image being processed, this coefficient will be modified in a different number of points (e.g., two similar images with different contrast will produce a different TD, as it depends on ∇u, see Fig. 5).It is easy to verify that the parameter D T Lim must be a value in the interval [D T MIN ,0 ) (being D T MIN the minimum value of the D T ).D T MIN and the distribution of values of the TD in this interval can vary considerably for different images.This will require readjusting this parameter for different images making its estimation more complex.
A different alternative, presented in Larrabide, Feijóo, Novotny & Taroco (2008), is to modify the diffusivity coefficient using a fixed point algorithm.This consists in sorting in growing order the values of the TD and selecting a percentage α of the most negative values and use this value as threshold for the insertion of cracks in the image.We then define the set M α as For s between one and the number of elements in the image.This alternative provides a better control on the algorithm (Algorithm 2).

Algorithm 2 Image restoration based on the continuous TD II
Require: Degraded image u 0 ∈ L 2 (Ω), parameters α e k 0 .Ensure: Restored image u * ∈ H 1 (Ω).compute u and q, solutions of the state and adjoint equation, respectively, compute the matrix 2 × 2 M and its minimal eigenvalue κ min for each point in Ω, find K using R D T -Continuous (Iso or Aniso) and Eq. ( 21), compute u * , the restored image, using K(x) previously obtained.

Continuous approach results
It can be observed in Fig. 3 the result obtained with Algorithm 2 for D T -Iso using the tensor K for the Lena image degraded wit white Gaussian noise (degraded image is presented in Fig. 3).For all the experiments presented k 0 = 2 was used.Visually, we observe that the noise is removed in the 3 cases.The results obtained in the three experiments, namely α = 0.10, 0.20 and 0.30, present considerable improvement in the SNR, going from 26.92 in the degraded image to 29.57, 30.36 and 30.88 in the processed images, respectively.By analyzing in detail these images, we observe that as other non-linear isotropic methods, it has difficulties to remove noise along edges.In  0.30) where used.Again, the SNR improves going from 26.92 in the degraded image to 28.47, 28.85 and 29.05 in the processed images, respectively.This time, and even if the SNR of the isotropic case are not reached, the noise along edges is more efficiently removed.Finally, in Fig. 6 is presented a detail of the TD, the vectors normal (in black) and tangent (in white) to the cracks and the restored image.

Discrete approach -R D T -discrete
The discrete approach relies on an auxiliary transient heat equation.In this case, cracks are introduced in the image to stop the diffusion in given points and directions.In this way, the information provided by the TD is used to determine the location of these cracks.In the discrete approach, the image is characterized by a matrix of pixels with an intensity associated to them.We consider a bi-dimensional image u given by a set of M × N pixels s.I n each pixel s, the image u intensity will be denoted as u s .Then, the image belongs to the space with Ω = ∪ s ω s , being ω s the domain of s.The set of neighbors n s of pixel s was defined as the four pixels 1 The cost function adopted by the discrete approach is which can be interpreted as a discrete approximation of the energy norm of the field u.I n this expression the term k s,p is the diffusion coefficient of pixel s with its neighbor p, n s = {n, s, e, w} are the neighbors of pixel s and ∆u s,p t is defined as In this case, u s t is explicitly computed as where the index t ≥ 1 represents the iteration number, being u s 0 the intensity at pixel s, k s = {k s,o , k s,l , k s,n , k s,s } characterizes the set of coefficients associated to pixel s, ∆t is the artificial step size in time.
As opposed to the continuous case, because u s t is an explicit function and given the discrete set k s , it is possible to compute the exact total variation of the cost function for each perturbation in k s,p .Also, we call k s ǫ the perturbed configuration of pixel s diffusivity coefficients.The value of the cost function when the perturbation is introduced is given by where D T (s, k s ǫ ) represents the total variation of the cost function due to a perturbation in the diffusivity coefficients of pixel s characterized by the set k s ǫ .Likewise in the continuous case, the introduction of a perturbation to pixel s where D T is negative, will produce a decrease in the value of the cost function J d .Using this information we can select the best candidate pixels for perturbations.We assume k s,p ∈{k ε , k 0 }, so the set of all possible configurations for k s is defined as C(s) :={k s =(k s,w , k s,e , k s,n , k s,s ); k s,p ∈{k ε , k 0 }, p = {w,e,n,s}}.
We see that 16 possible combinations for k s exist (values are k s,p = k ε or k s,p = k 0 , for each neighbor, then 2 4 = 16 cases are possible).The case k = k s,e = k s,n = k s,s = k ε is not taken into account as it does not change the value of the cost function.The 15 remaining combinations are 1 i.e. north, south, east e west of pixel s.  27) end for end for for each pixel s ∈M α do make D T (s)=min ǫ * {D T (s, k s ǫ ), k s ǫ ∈C σ } make k s = k s ǫ * the diffusivity coefficient associated to D T (s) end for compute u s t (k s ) using Eq. ( 24).
• no diffusion with one neighbor, • no diffusion with two neighbors sharing one vertex, • no diffusion with three neighbors, • diffuse on x direction, • diffuse in y direction, • diffuse in all directions.
The last case corresponds to isotropic diffusion and is defined as k s iso = {k 0 , k 0 , k 0 , k 0 }.To compute the value of D T for a determined pixel, its necessary to introduce a perturbation.This is done by changing, for one pixel s, the set k s for k s ǫ ∈C σ .Then, the cost function for u s t = u s t (k s ) and ǫ u s t = u s t (k s ǫ ) computed using Eq. ( 24) and ∆ ǫ u s,p t = u p t − ǫ u s t , respectively.
For Eqs. ( 25) and ( 26) the total variation of the cost function J d due to the perturbation k s ǫ is written as As in the continuous case, it will be considered a perturbation k s ǫ that minimizes the value of D T (s, k s ǫ ).Using this information is proposed the following discrete image restoration algorithm based on the TD 3).

Discrete approach results
The set M α is defined as in Eq. ( 21).As in the continuous case, the parameter α allows to control the values of the TD that will introduce changes in k s .In all the cases it was used △t = 1 4 e k 0 = 1, the maximum values that warrant the Courant-Friedrichs-Levy (stability) of the iterative solution algorithm.In Fig. 7 are presented some results obtained with this technique.The three images presented (corresponding to the results for values of α = 0.05, 0.15 e 0.25 respectively) present a considerable improvement of SNR (going from 26.91 in the degraded image to 28.49, 29.61 and 29.53 in the processed images, respectively).
In Figs. 8 are presented the results after different number of iterations.The edges introduced in the image are highlighted in white.We observe that after some iterations, the image has edges in almost all the edges.In this way, the variation of the cost function is practically null in two consecutive iterations, stopping the algorithm.The number of edges that are added to the image in each iteration is controlled by parameter α.In Fig. 9 is presented the detail of a region of the image before and after the processing.-Semi-quadratic minimization (Kornprobst et al. (1997)), -Continuous R D T (Belaid et al. (2007)).
In Fig. 10 are presented results for classical and TD image restoration methods.The processed image corresponds to the artificially degraded image of Lena (see Fig. 3) with uniform noise (r = 20), with an approximate SNR of 27 dB.All the classical methods depend on a parameter σ with equivalent meaning.This parameter is used to control the diffusion on the edges of the image to preserve relevant characteristics.For the classical methods, σ was adjusted using the technique proposed by Black et al. (1998), which allows to estimate the value of σ as a function of the gradient of the image processed.In the case of the evolutionary methods, the number of iterations was fixed to 10, 20 and 30.The results for 20 iterations where found to shield the lowest SNR and, thus, these are reported.For the semi-quadratic minimization, the same analysis was performed for parameter λ and the best SNR was obtained with λ ≈ 10 in this case.
Parameter selection in the case of TD methods is different.Both, the Discrete and Continuous R D T methods, depend on parameter α (a real value between 0 and 1), that determines the amount of cracks to be introduced in the image.This parameter determines the quantity (as a %) of the pixels that will be introduced in cracks to stop diffusion.As before, parameter α was analyzed to select the value that provides the best SNR (in the case of Discrete R D T α = 0.18 and for Continuous R D T α = 0.80.
As presented, the proposed restoration methods use information of the cost function sensitivity to a change in the topology.This information is used to find the optimal domain topology that, in the presence of diffusion, will eliminate noise preserving image features.In the continuous case, we observe that this technique eliminates most of the noise but has difficulties to remove noise in regions of elevated gradients.For the Discrete R D T , we observe that the noise is removed from the whole image, even from the edges.This algorithm is also capable of improving the sharpness of the edges, enhancing the image features.
In Table 1 are presented qualitative results for the Lena image.The different columns present the PSNR, SNR, μ(e) (mean error) and σ(e) (error standard deviation) between the processed image and the original one (i.e., without degradation).We observe similar results for the different methods, where the best performers were the evolutive method proposed by Black et al. with respect to PSNR (30.58dB) and SNR (30.37dB) and the one by Perona & Malik with respect to σ(e) (6.807).The proposed methods, Discrete R D T and Continuous R D T , provide results of similar quality to the existing methods.
The proposed methods provide an intuitive tool for image restoration based on the concept of the TD.These methods are intended to modify the topology of the image by inserting cracks in selected location that, in the presence of diffusion, will improve the quality of a degraded image.The diffusion will eliminate noise while preserving edges and details in the image.

Online material
The computational implementation of these methods in Matlab is available online at Matlab Central2 .A more complete description of the mathematical and numerical methods used in this work can be found in the work of Larrabide (2007).

Fig. 1 .
Fig. 1.Topological derivative -Change in topology 2.1 Topological shape sensitivity analysis Let a problem be defined in Ω where the quality/performance is characterized by a cost function J (Ω)=F (Ω, u(Ω)) where Ω ⊂ R n , n = 2, 3, is a regular domain of open boundary ∂Ω with exterior normal n.With the notation (Ω, u(Ω)) we empathize that F depends on Ω explicitly and implicitly through u(Ω), solution of the variational equation (state equation), that can be written in the abstract form as: determine u ∈U= U (Ω) such that

100
Image Restoration -Recent Advances and Applications www.intechopen.com

Fig. 2 .
Fig. 2. Topological Derivative concept for cracks.being M a symmetric tensor given by

104
Image Restoration -Recent Advances and Applications www.intechopen.com
(a) Lena image.(b) Lena image with low contrast.(c)Topological derivative for the previous images.

Fig. 5 .
Fig. 5. Lena image with different contrast.The second row corresponds to the TD for each image using the same color scale in both cases.

Algorithm 1
Image restoration based in the continuous topological derivative -R D T -Continuous Require: Degraded image u 0 ∈ L 2 (Ω), parameters D T Lim e k 0 .
Fig. 4 are presented the results corresponding to K computed using D T -Aniso and the same Lena image.The same values of α (namely α = 0.10, 0.20 and 108 Image Restoration -Recent Advances and Applications www.intechopen.com(a) Detail of degraded Lena image.(b) Detail of the TD.(c) Detail of the crack normal and tangent directions.(d) Detail of the restored image.

110
Image Restoration -Recent Advances and Applications www.intechopen.comAlgorithm 3 Image restoration based on a discrete version of the TD -R D T -Discrete Require: The 2D image u 0 ∈U, a diffusivity coefficient k 0 and a parameter α.Ensure: The restored image u s ∈U.make t=1; stop = FALSE; k s = k s iso , s = 1..M × N while stop = FALSE do for each pixel s do for each k s ǫ ∈C σ do compute D T (s, k s ǫ ) following Eq.(

Fig. 7 .
Fig. 7. Image restoration with the R D T -Discrete algorithm of the Lena image.Results correspond to k 0 = 1 and α = 0.05, 0.15 e 0.25 respectively.

Fig. 8 .
Fig. 8. Detail of the crack configuration introduced to the image during the processing (k 0 = 1 and α = 0.20).

Table 1 .
Comparison between the proposed and classical image restoration methods.