- About this Journal
- Abstracting and Indexing
- Aims and Scope
- Annual Issues
- Article Processing Charges
- Articles in Press
- Author Guidelines
- Bibliographic Information
- Citations to this Journal
- Contact Information
- Editorial Board
- Editorial Workflow
- Free eTOC Alerts
- Publication Ethics
- Reviewers Acknowledgment
- Submit a Manuscript
- Subscription Information
- Table of Contents
Abstract and Applied Analysis
Volume 2013 (2013), Article ID 816134, 20 pages
Hybrid Topological Derivative-Gradient Based Methods for Nondestructive Testing
1Departamento de Matemática Aplicada, Universidad Complutense de Madrid, 28040 Madrid, Spain
2Departamento de Fundamentos Matematicos, Universidad Politécnica de Madrid, 28040 Madrid, Spain
Received 26 December 2012; Accepted 9 June 2013
Academic Editor: Ağacik Zafer
Copyright © 2013 A. Carpio and M.-L. Rapún. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This paper is devoted to the reconstruction of objects buried in a medium and their material properties by hybrid topological derivative-gradient based methods. After illustrating the techniques in time-harmonic acoustic problems with different boundary conditions and in electrical impedance tomography problems with continuous Neumann conditions, we extend the hybrid method for a realistic model in tomography where the boundary conditions are given at a discrete set of electrodes.
In many practical aspects of life we face the need of reconstructing the geometry of objects buried in a medium and identifying their material properties. In medicine, precise detection systems are required to locate tumors, clots, or damaged tissues without harming the patient. Magnetic resonance imaging together with different kinds of tomography, X-rays, or ecography have become powerful diagnosis tools. Geophysicists have long relied on measuring scattered waves to obtain information on the structure of the ground: from locating water, oil, and gas reservoirs or mineral mines to identifying the nature of earth quake sources or fault structure. Assessing the structural integrity of bridges, planes, power plants, and buildings after natural disasters also requires sharp techniques for nondestructive testing to identify cracks or damaged areas.
A classical strategy to obtain information about the inner structure of a medium is to illuminate it with some kind of radiation; see Figure 1. The total wave is usually measured at a set of receptors distributed over a curve or surface . From these measurements, we aim to reconstruct the inner geometry of the medium and their material coefficients. Objects are distinguished from the background because they differ in their material parameters. The nature of the incident radiation is chosen to enhance that contrast. Acoustic waves are used when the parameters of interest are the elastic constants. Electromagnetic scattering tracks the electric conductivity and permittivity. Thermal waves allow to distinguish objects by their thermal diffusivity and conductivity.
There are two important mathematical problems involved in this process: the forward problem and the inverse problem. In the forward problem, the material properties of the background medium and the shape, size, location, and material parameters of all the inclusions are known. The goal is to compute the wave field at the receptors. It is usually a well-posed boundary value problem with a unique solution that depends continuously on the initial data. The equations governing the wave field will be Navier, Maxwell, or heat equations depending on the nature of the chosen incident radiation. In the inverse problem, measurements of the wave field at the receptors are known. The objective is to identify the geometry and nature of the objects located inside the background medium. Those objects are characterized with the fact that solving the forward problem using the known incident wave and placing those objects inside the medium, the resulting wave field evaluated at the receptors agrees with the measured data. This inverse scattering problem is nonlinear, and severely ill posed. It may not have a solution, the solution may not be unique, or it may not depend continuously on the data.
The inverse problem can be reformulated in a more regular way as a constrained optimization problem. In practice, one does not need to know the exact geometry and material parameters of the objects for which the wave field evaluated at the receptors agrees exactly with the measured data. It is enough to find approximations making the error small enough. This motivates a weaker variational reformulation of the inverse problem: Find domains , and parameters minimizing where is the solution of the forward problem with objects immersed in the exterior media and parameters ( stands for “interior”). The design variables are and . The forward problem is the constraint.
Developing mathematical tools to solve this type of inverse problems is an active field, in need of even better techniques. Many beautiful mathematical theories have been proposed to produce reasonable reconstructions in a large variety of contexts from different standpoints, but the demand of more refined methods remains. Some of the first theories relied on linear approximations, such as the Kirchhoff or physical optics approximation  or the Born approximation . Linear sampling methods became a step forward [3, 4]. These methods formulate an equation for a parametrization of the contour of the unknown domains in terms of its representative coefficients. Backpropagation principles  and modified gradient methods  have also been employed. More recently, different optimization strategies have been proposed to deform initial guesses of the contours of the objects in such a way that adequate cost functionals decrease. The techniques differ by the way the objects are parameterized and deformed. Classical shape deformation along vector fields was tried first [7–10]. This strategy has a drawback: it does not allow for topological changes in the objects. The number of boundaries (connected components, holes, etc.) has to be known from the beginning, and convergence of the method is unlikely without a priori information. Representing the objects as level sets of level set functions deformed following Hamilton-Jacobi equations, with velocities chosen to ensure the decrease of the cost functional, solves that drawback [11, 12]. Topological changes are allowed during the deformation process. However, level set methods need an external guess to start, and convergence may be slow unless the guess is good enough. Topological derivative techniques have arisen as a promising alternative. They have the advantage of providing good first guesses for the number, size, and location of the objects, serving also as a basis for iterative schemes that improve those initial approximations. Their potential to provide good first approximations was explored in  for exterior acoustic problems involving rigid bodies, in [14, 15] for three-dimensional transmission problems, and in  for waveguides with Dirichlet boundary conditions. Reference  discussed extensions to scattering by elastic waves. Reference  and references therein suggested possible applications in electromagnetism and tomography. A hybrid level set-topological derivative scheme was proposed in . Iterative schemes entirely based on topological derivatives were introduced in  and later developed in [21–23]. Modified cost functionals allow to prove theoretical convergence results and study the effect of noise [24, 25].
A vast amount of work in the field refers to time-harmonic incident waves, for which the boundary value problem governing the wave field becomes stationary. Less variety of methods are available when the incident waves are not time harmonic, for instance, for pulse-like excitations. For incident waves governed by time-dependent wave equations, see the classical work  and references therein. Methods for detection of inclusions near the surface by photothermal techniques are developed in .
The previuosly mentioned references focus mostly on reconstructing the distribution of the scatterers. In practice, one also needs to predict their material parameters to know their nature. In medical applications, for instance, this is essential to be able to distinguish benign from malignant tumors. There is a vast amount of the literature about reconstructing the spatial variations of parameters dealing in particular with electrical impedance tomography; see the reviews [28, 29] and references therein. The complex case and the problem with discrete electrodes are much less studied. Schemes following discontinuities or objects with sharp interfaces usually consider a piecewise-constant conductivity [30–35]. Some of them seek for the objects without being able to precise the value of the conductivity [30, 36]. Others may require a large number of iterations [31–33, 35]. Level set techniques have been recently extended to track interfaces between regions of different conductivity [31, 35]. These methods are initialized resorting to topological derivatives in [32, 33]. For smooth conductivities, [37–39] describe promising methods.
In elastic scattering, spatial dependence is unavoidable if the background matrix is heterogeneous or as a way to detect damage; see  for tests in functionally graded reference materials. Hybrid topological derivative-gradient based methods to reconstruct both the geometry of objects and their material parameters were introduced in [21, 22] for acoustic problems and in  for the electrical impedance tomography problem with continuous Neumann conditions.
In this paper, we explore the potential of some specific optimization techniques to solve imaging problems. We focus on hybrid descent techniques that combine gradient methods to approximate the parameters and topological derivative based schemes to calculate the domains. We review some recent results and present new developments for tomography. The paper is organized as follows. Section 2 exemplifies the hybrid schemes for time-harmonic incident waves. We focus on acoustic waves, though some types of polarized electromagnetic waves lead to inverse problems with exactly the same mathematical structure. Electrical impedance tomography is considered in Section 3, both the continuous and the discrete versions of the electrode model. We reformulate the problem with discrete electrodes as a constrained optimization problem. Expressions for the topological derivative with respect to the inclusions and the correctors of their parameters are given in terms of adjoint fields satisfying adequate variational equations.
2. Methods for Inverse Acoustic Problems
In this section, we deal with inverse acoustic scattering problems. We assume that the incident wave is time harmonic; that is, where is the frequency. In this case, after a sufficiently long time, the total wave field is also time harmonic, , and the total amplitude is governed by a boundary value problem for the Helmholtz equation.
Depending on the nature of the objects, different boundary conditions model the interaction between the scatterers, the medium, and the incident wave. For penetrable objects transmission conditions are imposed at the interfaces between the objects and the surrounding medium, whereas Dirichlet and Neumann conditions model sound-soft and sound-hard objects, respectively.
We assume that the incident wave is a planar wave propagating in the direction , with wave number . Other types of incident waves can be considered, for instance, in , models a point source excitation at . Here, is the Hankel function of the first kind and order zero (see ).
Let us begin with describing the geometrical configuration of our inverse acoustic scattering problems (see Figure 1). We have a finite number of bounded and simply connected open objects having no pairwise contact. Their boundaries are assumed to be smooth, for instance . For simplicity, we will use the notation . These objects are immersed in the exterior media . We focus on two-dimensional problems to illustrate the methods, though the theory extends to any dimension.
Let us start by considering penetrable nonhomogeneous obstacles. To simplify, we assume that the medium is homogeneous. Then, the total wave field where and denote the scattered wave outside the obstacles and the transmitted wave inside , respectively, is the solution of Here, is the wave number of the incident wave and is related to the material properties of the inclusions . We assume that is differentiable. When is constant inside each , we are dealing with homogeneous materials.
The condition at infinity is the standard Sommerfeld radiation condition with , which has to be satisfied uniformly in all directions. It implies that only outgoing waves are allowed. denotes radial derivatives. stands for outward normal derivatives, with pointing outside . and denote the limits of from the exterior and interior of , respectively. The continuity condition for the normal derivative can be extended to when is replaced by in and by in , as explored in .
The optimization problem we want to solve is then as follows: Find and minimizing where is the solution of the forward problem (4) with objects and parameters . This cost functional is positive but vanishes for the true scatterers and the true parameters, which are characterized as a global minimum of the functional.
A common strategy to approximate solutions of optimization problems is to resort to descent techniques. The idea is to find initial guesses for the parameters and the objects and then generate sequences of approximations , along which the functional decreases. In the next subsections, we explain how to generate those sequences in several stages. First, we assume that the domains are known and present a gradient strategy to guess . Next, we will assume that the parameters are known and discuss how to approximate by topological derivative techniques. Then, we will combine both methods performing a double iteration with respect to the domains and the parameters to approximate both the objects and their material constants. Once the method is presented for penetrable objects, we will discuss how it is modified to handle rigid or sound-soft obstacles.
2.1. Gradient Techniques to Identify Coefficients
Let us assume that the objects are known. The optimization problem becomes as follows: Find minimizing where is the solution of the forward problem (4) with objects and parameters .
Fixed an initial guess , we generate a sequence providing a decreasing sequence by a gradient method . At each stage, we fix the current guess and study the variations of the functional along directions : , where is a real parameter. The derivative with respect to is given by being the solution of the forward problem (4) with coefficient inside the objects . Computing is avoided introducing an adjoint field ; see .
Theorem 1. The derivative with respect to of the function defined by (6) evaluated at zero is given by where is the solution of the adjoint problem and the solution of the forward problem (4) setting in both cases .
Choosing in such a way that the derivative , we ensure that for small enough. For instance, we may set The updated guess is then This choice is reasonable if we are looking for space-dependent parameters. When is composed of several disjoint domains , and is known to be almost piecewise constant, we may set Technical details and proofs can be found in , where more general Helmholtz type constraints are considered: The coefficients , , , are bounded from below by positive constants and depend smoothly on , becoming constant at infinity. The constant at the Sommerfeld radiation condition is defined as . It is the wave number of the incident wave. The interior parameter can also be identified following the same technique.
2.2. Topological Derivative Techniques to Reconstruct Domains
We have seen how to generate descent directions for the parameters fixing the objects . Let us assume now that we know the parameters . The optimization problem becomes as follows: Find minimizing where is the solution of the forward problem (4) with objects and parameters . To generate descent directions for the domains, we must rely on some kind of derivative with respect to the domains. We choose a specific type of derivative for shape functionals: the topological derivative. It has the advantage of providing a first guess for the objects too.
The topological derivative of a shape functional at the point is defined as  where is a ball centered at with small radius ; see Figure 2. The scalar function is chosen in such a way that the limit (15) is finite and does not vanish. For instance, for the transmission or Neumann conditions in two dimensions, we take . For Dirichlet conditions in 2D, we may choose .
The topological derivative is a scalar function defined for every which measures the sensitivity of the functional to removing from a small object located at . In regions where takes large negative values, we have a large probability to find an object. This remark is based on the following expansion: Whenever , for small. When is our functional (14), this suggests a strategy to find first guesses for the objects and then improve them iteratively. Our algorithm to recover objects when the parameters are known is as follows:(i) Compute the topological derivative when and choose an initial guess for the objects as where is a positive constant. (ii) Check that . Otherwise increase the constant . (iii) For (a) Compute the topological derivative when and select for decreasing thresholds . (b) Check that . Otherwise increase the constant . (c) Check the stopping criteria: the iteration stops if either meas is small enough, or is small enough, or the difference between consecutive values of the functional is small, or the discrepancy principle is satisfied, where describes measurement errors and (in our numerical experiments ). Here, is the solution of the forward problem when the objects are .
Computing topological derivatives at grids of points at each stage using (15) is computationally very expensive. Explicit formulas that require only the knowledge of forward and adjoint fields are used instead. For our particular Helmholtz type constraint (4), the explicit expression is [14, 21].
Notice that when and are the true objects and parameters for which the data are measured, the adjoint field vanishes since the forward field agrees with the measured data at the receptors. Therefore, the topological derivative vanishes at the exact solution, which is a global minimum of the cost functional, as explained before.
In real reconstruction tests, one usually has access to measurements for several incoming directions , . The shape cost functional (14) is then replaced with: Here, is the solution to the forward problem (4) with incident field and is the measured total field on . The topological derivatives for this functional are obtained adding up the topological derivatives for each single incident wave.
We test the algorithm for a geometrical configuration with two objects. The interior parameter in them is whereas . We have computed the topological derivative for data given at 14 sampling points uniformly distributed on the circle centered at with radio 3, for 10 incident planar waves with uniformly distributed angles in . The geometrical setting for our experiments is represented in Figure 3. Our synthetic data were generated by solving the corresponding forward problems by boundary element techniques using a finer mesh and adding a one percent relative random noise at each observation point.
The values of the topological derivative when at the sampling region are represented in Figure 4(a). The largest negative values are attained inside each object. The first computation of the topological derivative provides the correct number of objects and their location, although they look like balls. Our initial guess is represented in Figure 4(b). We compute now the topological derivative when , obtaining the color map in Figure 4(b). New regions close to the current objects appear where the topological derivative attains large negative values. Those regions should belong to the scatterers. Our objects are bigger and they are not balls. Figures 4(c)–4(j) represent the subsequent iterations. The values of the cost functional versus the number of iterations are given in Figure 4(k). After few iterations, our reconstructions are satisfactory.
The quality of the reconstructions as well as the number of iterations depends on the values of the constants appearing in (18). Guidelines for their selection are given in . When are too big, the number of iterations required is high. A sequence of smaller provides faster reconstructions, but we are prone to include spurious regions in that cannot be removed in the next iterations. We illustrate this fact in Figure 5. Figure 5(a) is exactly the same as Figure 4(a) and represents the values of the topological derivative when . Choosing a smaller constant , our initial guess , represented in Figure 5(b), has two objects that are bigger than the ones in Figure 4(b). The object on the right includes too many spurious points. Computing the topological derivative when , we detect a small region at the top of the left object that is included in the next iteration (see Figure 5(b)). No improvement for the rightmost object is pointed out. The reconstruction is given in Figure 5(c). In the next iteration (see Figure 5(d)), only the object on the left is slightly modified and the algorithm stopped because the difference between the consecutive iterations in the value of the functional was negligible, as observed in Figure 5(e). The quality of the reconstruction of the object on the left is reasonable, but, comparing with the previous example, the reconstruction of the object on the right is not acceptable.
The monotone algorithm cannot remove spurious regions in the rightmost object. The topological derivative is only defined at points not belonging to the current approximation. However, observing the values of the topological derivative in Figures 5(b)–5(d), we realize that there are regions close to object where it takes large positive values. Computing the topological derivative when (see Figure 5(a)), we did not see that effect. This indicates that spurious regions are included and that we should use a nonmonotone method if we want to correct them.
The classical notion of topological derivative we used previously is defined only in . Therefore, any time we update our approximations , we have a criterion to add points. However, if at some stage we include in our approximation points that do not belong to the objects, we cannot remove them. This artifact is eliminated observing that for our particular type of functional we can extend the definition of topological derivative to points inside the objects. The topological derivative of a functional was extended to the points inside in  as From this definition, the expansion below follows: Since is a positive function, if is large and negative, the functional decreases when we remove the ball from . This suggests a strategy to remove points from . Our algorithm can generate a nonmonotone sequence of objects by substituting the definition of in the previous method (see (18)) with with decreasing sequences , .
Notice that the expressions (19) and (24) are identical except for a minus sign. Moreover, since the solutions to the problems (4) and (9) are continuous functions in , the expression defines a continuous function in . This motivated in  to modify the extended definition of the topological derivative inside of a functional as which is exactly the same as (21) multiplied by . Using this notion, and whenever is large and positive, the functional decreases when removing the ball from . This definition also allows to define a non-monotone scheme by substituting the definition of in (23) by for decreasing thresholds , .
With this definition, the topological derivative of the functional for the transmission problem is as follows.
There is a missprint in . Definition (21) was given, but in fact (26) was used. The topological derivative inside was computed as in (29) and in the numerical examples the implemented formula is (29). The modified algorithm described by (28) was used to update the domains.
Let us revisit now the example of Figure 5. The topological derivative for provided the initial guess represented in Figure 6(b) by a dashed line. The values of the topological derivative inside in Figure 6(b) are calculated using the definition (21), that is, formula (24). The object on the right includes spurious regions (the points with large negative values inside ) that will be removed in the next iterations. Some points at the top of the object on the left should be included in that object. If we use the alternative definition (26), the topological derivative is a continuous function (see Figure 6(c)). A region inside the rightmost object where large positive values are attained should be removed. A region outside the leftmost object where it takes large negative values has to be included. The second definition provides a color map that is easier to interpret. Therefore, from now on we choose the definition (26) and update the domains using (28). The final reconstruction at the 11th iteration is represented in Figure 6(d) by a dashed red line. For comparison, we have included in Figure 6(e) the final reconstruction obtained for the same problem in Figure 4 with the monotone method (for a well-calibrated sequence of constants ) after nine iterations. Our non-monotone method is not very sensitive to the initial guess.
A wide gallery of reconstructions using the monotone and the non-monotone methods can be found in [21, 43], illustrating the ability of the method to reconstruct objects with different parameters , (simulating high and low frequencies), problems with different values of inside each object, scatterers of different sizes, different levels of noise in the data, reconstructions when data are given on a limited region of a circle, source point excitations instead of planar incident waves, and so forth.
2.3. Hybrid Approach to Reconstruct Objects and Identify Their Material Parameters
Combining the two strategies to approximate objects and parameters discussed in the previous two sections, we obtain a reconstruction scheme [21, 22]. (i) Initialization (a) Choose an initial guess for the parameters , for example, a perturbation of the background parameter . (b) Choose an initial guess for the objects , using the topological derivative: set and compute the topological derivative with object : where and are forward and adjoint fields with object and parameter . Then, where is a certain large positive threshold. (ii) Iteration(a) Update using the gradient technique: where , solve forward and adjoint problems with objects and parameter . (b) Check that . Otherwise, reduce .(c) Update computing the topological derivative with objects and parameter . We may resort to a monotone or nonmonotone strategy with decreasing thresholds , . The second one is more convenient if we wish to be able to remove spurious points at any stage, for instance, when reconstructing objects with holes. (d) Check that . Otherwise, increase the thresholds , . (e) Stop when meas is small enough, or is small enough, or the difference between consecutive values of the functional is small, or the discrepancy principle is satisfied.
The iterative procedure described perviously generates sequences of domains and parameters along which the functional decreases. The magnitude of the topological derivative is also observed to decrease. Theoretically, convergence to the true solution is not guaranteed, since our cost functional might have local minima different from the global minimum we are looking for. This situation can be improved or worsened by the amount of incident waves we use, the number of receptors, and how we distribute them. In the tests we have performed, the thresholds were calibrated so that the scheme worked quite well in a few iterations. See [21, 22] for details about calibration, numerical schemes and extensions to constraints of the form (13).
The forward and adjoint problems have explicit solutions when is constant and we are computing the first guess on the whole space. They become The solution to (35) is the incident wave and the solution to (36) can be computed using the fundamental solution of the Helmholtz equation. The forward and adjoint problems (4) and (9) are solved by either boundary element or finite element methods depending on whether is piecewise constant or vary smoothly with space [11, 13, 21, 44–46].
In practice, updating the domains is more expensive than updating the parameters. Any time we change the guess of the objects, we must mesh again the new domains to be able to compute forward and adjoint fields. Therefore, we usually correct the parameters several times for a fixed approximate domain before proceeding to update , to reduce the computational cost.
We have tested the hybrid topological derivative-gradient based method considering the configuration with two objects of the previous examples, keeping the same observation points and incident waves (see Figure 3). We also take the same values of the parameters: is assumed to be known and is assumed to be unknown. We computed the topological derivative when choosing as initial guess for . Notice that since the forward and adjoint fields when do not depend on , the topological derivative is the same as in Figure 4(a) except for the scale (we have now in (19) the multiplying factor ( instead of ). Therefore, our initial guess , represented in Figure 7(a), is the same as in Figure 4(b). We fix now and iterate five times to update the value of without modifying the scatterers. Then, we fix the value of and compute the topological derivative to update the domains and so on. The reconstruction at the ninth iteration with respect to the domain is shown in Figure 7(b). The values of through the iterative procedure are represented in Figure 7(c). Two identical values of the parameter in the plot mark each iteration to improve the domains ( is not updated). Our final approximation is 1.23 while the true value is 1. Some reconstructions with two or three objects with different constant parameters can be found in [21, 43].
In the next example we test the hybrid algorithm for the reconstruction of an object with space-dependent coefficient . The scatterer is the ball centered at and radius 0.35 and is the function represented in Figure 8(c). It grows radially from 1 to 2 inside . The observation points and the incident waves are the same as in all the previous examples. To start the algorithm, we took as initial guess the constant function . The topological derivative for provided the initial guess represented by a red dashed line in Figure 8(a). We applied the hybrid method alternating one iteration with respect to the domain with five steps of the gradient method. The final reconstructions of and are shown in Figures 8(b) and 8(d), respectively. The parameter is overestimated with an average error about 0.2, but we clearly recognize a function that grows radially. The location, shape, and size of the object are reconstructed with accuracy. The interested reader can find in  some numerical experiments including the reconstruction of two homogeneous objects immersed in an heterogeneous material and of two heterogeneous scatterers buried in an homogeneous media.
2.4. Sound-Hard and Sound-Soft Objects
The methods presented in the previous sections apply when the objects are penetrable by the incident radiation. Part of the incident wave is transmitted inside the object and transmission boundary conditions are imposed at the interface between the objects and the surrounding medium.
When the scatterers are sound-hard (rigid), no transmitted wave is generated. The incident radiation is scattered and a Neumann boundary condition is imposed at the surface of the objects. The constraint for the cost functional (14) becomes Topological derivative methods to find first guesses of the number, shape, and location of rigid objects were developed in [13, 44]. However, to the best of the authors’ knowledge, iterative methods entirely based on the computation of topological derivatives for the reconstruction of sound-hard objects had not been previously used.
For sound-soft scatterers, a Dirichlet boundary condition is imposed at the interface between the objects and the background medium: Topological derivative methods were first employed to reconstruct sound-soft objects in [16, 20, 47].
For rigid and sound-soft objects, we cannot predict their material parameters, since they do not appear in the constraints. The cost functionals only depend on the geometry of the inclusions . The topological derivative based iterative methods sketched in the previous sections can be extended to these two types of objects, replacing the explicit expressions of the topological derivatives with other adequate formulas, which involve modified forward and adjoint fields. The forward and adjoint problems are formally similar to the transmission ones, suppressing the equation for the transmitted wave inside and imposing either Neumann or Dirichlet boundary conditions on .
The new expressions for the topological derivatives are: (i) For the Neumann problem as follows where the forward field is a solution of (37) and the adjoint field is a solution of This formula was obtained in , correcting the previous paper  where the multiplying factor 2 in (39) was missed. (ii) For the Dirichlet problem (see ) where the forward field is a solution of (38) and the adjoint field is a solution of
We test the performance of our iterative topological derivative method for rigid and sound-soft objects in Figures 9 and 10, respectively. The geometrical setting, sketched in Figure 3, is the same as in the previous examples for the transmission problem with constant coefficients. We also keep the same observation points and incident waves. The (exterior) wave number is too.
In Figure 9, we illustrate the behavior of the method when the objects are sound-hard. Computing the topological derivative when (Figure 9(a)), we detect the correct position and approximate size of the two objects. The initial guess is represented in Figure 9(b). After nine steps of the iterative procedure, the reconstruction is rather satisfactory as shown in Figure 9(c).
When applying the method for sound-soft objects (see Figure 10(a)) the topological derivative when attains the largest negative values in a region inside the object on the right. This means that in this case the topological derivative ignores the object on the left and the first approximation consists of one object, the ball , represented in Figure 10(b). In the same plot we also represent the values of the topological derivative when . The object on the right is bigger. The object on the left is ignored again. In the next iteration (see Figures 10(c)-10(d)), we find the leftmost object. After thirteen steps both objects are reconstructed with accuracy, see Figure 10(e). A number of theoretical results [47–49] suggest that good reconstructions of sound-soft objects can be achieved with just one or very few incident waves. Numerically, this was tested in .
Some inverse scattering problems involving polarized electromagnetic waves (TM or TE) can be reformulated as constrained optimization problems similar to the ones studied in Section 2. The constraints are Helmholtz type equations, eventually with complex wave numbers . The methods described in the previous section apply, with a slight variation of the formulas for topological derivatives about handling constraints involving complex wave numbers .
3. Methods for Impedance Tomography
In this section, we deal with another imaging technique that uses time-harmonic electromagnetic waves: impedance imaging.
A derivation of the impedance imaging problem from Maxwell’s equations is given in . Inside the body , the electric potential satisfies where is the admittivity, , is the electric conductivity, is the electric permittivity, and is the angular frequency of the applied current. When , we are left with an inverse conductivity problem.
Currents are applied through electrodes located on the boundary of the body, where the resulting voltage is measured. In real experiments, only the values of the current and the voltages at a “discrete” set of electrodes are known. Realistic models supplement (43) with a “discrete” Neumann boundary condition: where satisfies is the current sent to the th electrode and denotes the part of that corresponds to the th electrode. We assume that each is an open subset of with positive surface measure and that for . and have smooth boundaries. is continuous up to the boundary, is continuous, and and define functions.
At the electrodes, we measure the voltages . Dirichlet conditions have been shown not to be realistic  and are commonly replaced with the constraints where is the effective contact impedance or surface impedance and As shown in , problem (43)-(44) admits the variational formulation: Find such that where The subscript refers to vectors of vanishing mean; that is, . Existence and uniqueness are guaranteed by the Lax-Milgram lemma for bounded coefficients and smooth domains. A discretization scheme by FEM has been proposed in ; see also references therein.
A simpler model imposes a “continuous” Neumann boundary condition: and measures the voltage on the boundary
Here, denotes the outward unit normal and stands for the outward pointing normal component of the current density generated by the electrodes on the surface. This Neumann problem (43)–(51) admits solutions when the current satisfies the compatibility condition , which turns to be the law of conservation of charge. To select a unique potential , we impose . The measured voltage also satisfies the condition .
The goal is to reconstruct the structure of the admittivity inside from measurements at the boundary. Two slightly different settings may be considered. If contains a number of inclusions , the admittivity is a piecewise smooth function in with discontinuities at the boundaries of the inclusions. The admittivity outside the objects is known. The admittivity inside is unknown and must be reconstructed, together with the geometry of , also unknown. Otherwise, one may consider admittivities to vary smoothly over , being only known at its boundary. For brevity, we assume here that the background contains a number of inclusions. The second setting was studied in . First, we illustrate the idea applying hybrid topological derivative-gradient based techniques to the simplified model (43)–(51) following . Reasonable reconstructions of the geometry of the objects are obtained in a few iterations, though the values of the admittivity inside inclusions might be better estimated replacing the gradient approach with other techniques [37–39]. Then, we extend the theory to the realistic tomography problem with discrete electrodes, much less studied.
3.1. Tomography with the Simplified Model
The “continuous” impedance tomography problem can be reformulated as an optimization problem. We assume that contains a number of inclusions and the admittivity is a piecewise function in with discontinuities at the boundaries of the inclusions. For simplicity, we take to be a bounded set of with smooth boundary. We set with simply connected open bounded sets with smooth boundaries without pairwise contact; see Figure 11. The admittivity in the background medium is known and is unknown inside . The functional to be optimized becomes where is a solution of The unit normal points outside but inside . and denote the limit values of on from outside and inside , respectively; see Figure 11.
As observed in Section 2, to apply hybrid topological derivative-gradient based methods to this functional, we need two ingredients: simple expressions for its topological derivatives and formulas for variations with respect to the admittivity. Formulas for topological derivatives of this functional are given in [23, 32, 33].
When we fix and perturb guesses of the admittivity in a certain direction, the derivatives of the resulting functional are given below. See  for the proof. Related differential calculus and formulas are presented in [28, 37, 52].
Choosing we ensure if is small enough. We have defined only in , but formula (59) holds for all . When is piecewise constant, we expect to be almost piecewise constant. In this case, another choice for the corrector is a piecewise constant approximation:
The reconstruction algorithm for the impedance tomography problem is similar to the algorithm described in Section 2.3 with , playing the role of , , and using (55)-(56) for the topological derivatives, (59) or (60) for the parameter correctors. We test this method in the next experiment where is the unit ball. Data are given at 30 uniformly distributed electrodes on (represented by crosses in Figure 12) for the ten current patterns We consider complex admittivities of the form with . Inside the two objects represented by solid blue lines in Figures 12(a)-12(b), the unknown admittivity is and is assumed to be known. We started the hybrid algorithm with . The first guess for the domains was obtained observing the topological derivative with and . It has only one object, as shown in Figure 12(a). We applied the algorithm alternating one iteration with respect to the domain with five iterations with respect to the parameter. Although we initialized the method with an incorrect number of objects, after 17 iterations with respect to the domain we obtained a reasonable reconstruction of the two defects; see Figure 12(b). When the algorithm stops, the true admittivity is approximated by . The values of the conductivity (the real part of ) and of the permittivity (the imaginary part of ) throughout the iterative procedure are represented in Figures 12(c) and 12(d), respectively. We refer to  for an example of the reconstruction of an object with unknown spatial dependent admittivity. The interested reader can also find in  a gallery of reconstructions mainly focused on the conductivity problem (with ).
3.2. Tomography with Discrete Electrodes
As before, we assume that contains inclusions of a different material. For simplicity, we work in the two dimensional setting, but the method remains valid in 3D with straightforward modifications. The admittivity of the background medium is known. To determine its spatial variation inside the inclusions we must optimize (62) with respect to the interior domains and their admittivity. As shown before, even if the admittivity is unknown, we may locate the regions where it undergoes noticeable changes by computing the topological derivative of the shape functional in the whole domain. Once a first guess for is found, we may approximate the admittivity optimizing with respect to for fixed by a gradient technique. The guesses for can be updated computing new topological derivatives. Below, we give expressions for the correctors to be used in the gradient method and for the topological derivatives.
Let us fix and seek to minimize when is a solution of the forward problem with discrete electrodes (63). When the nature of the inclusions we look for is known, an initial approximate value is usually available. Otherwise, we select a first approximation by perturbing . We may improve this first approximation iteratively by a gradient method. Given a guess , we generate a descent direction to improve this guess studying the function of a real variable . As in the acoustic setting and in the “continuous” model, and are selected in such a way that . An explicit expression for this derivative provides the following formula for the descent directions .
Theorem 7. The derivative with respect to of the function , with defined in (62) is where is a solution of the forward problem (63) with parameter in and is a solution of the adjoint problem (65) given below where with in . Choosing and small enough, one ensures .
Proof. The derivative with respect to is
where and is a solution of (63) with ; that is,
and is defined in (50). Notice that when and , the sesquilinear forms (67) and (72) are the same.
We eliminate by introducing the adjoint problem. Let us set for . Then, choosing , it follows that and therefore for all . Selecting as a solution of (65), we obtain the desired formula
Notice that if , then and . Thus, the integrals appearing in the definition (68) of are in fact duality products.
We have defined only in , but formula (69) holds for all . When is piecewise constant, we replace defined in (69) with piecewise constant approximations inside each object , An initial approximation of the objects buried in the medium as well as their successive corrections is computed using the topological derivative. For the tomography problem with discrete electrodes, the explicit expression for the topological derivative in terms of forward and adjoint fields is given in the theorem below. Combining the expressions for the approximation of the admittivity with the topological derivative computation to approximate the domains , we can define a hybrid reconstruction scheme for the discrete electrode model as in the previous sections. The implementation of the numerical codes is out of the scope of this paper and will be the subject of our future work.
Proof. We adapt the strategy used in [13, 21] for inverse scattering problems. The topological derivative of a shape functional can be computed as a limit of shape derivatives . Consider the family of deformations depending on the real parameter ,
where the vector field is defined as
and extended to in such a way that it vanishes away from a narrow neighborhood of . Then,
where is the derivative of the function that appears in the definition of the topological derivative (15). For the electrical impedance tomography problem in 2D, we can take as we will shortly see.
The proof is organized in three steps. We assume first that and are constant and compute the topological derivative when . Afterwards, we adapt the proof to smoothly varying coefficients. Finally we consider the case .
Step 1 (computation of the topological derivative when and and are constant parameters). Let us compute first Since the field decays fast enough away from , it follows that and , where . The cost functional in the transformed domains is then where solves with defined in (67) and defined in (50). Notice that when , , and the solution of (84) is the solution to (63) with .
The derivative of (83) with respect to involves the derivative of with respect to . To avoid its computation we introduce the Lagrangian functional Then, for all , If satisfies (65) with object and parameter , then To obtain the last equality in (87), we have integrated by parts and used that and solve (63) and (65) with , respectively, and that vanishes on .
By (81), the topological derivative at can be computed now as The asymptotic behavior of and at is obtained expressing them as corrections of and , the solutions to (63) and (65) with , and expanding the remainders in powers of . Reasoning as in  we find that uniformly when . For there is a slight technical difference because is defined by the variational problem (65). However, this is due to the source term at , that is included in . Once that source term is eliminated, the remainder satisfies equations similar to the equations for and we find uniformly when .
Replacing the limits (90) and (91) in (88) and choosing , formula (78) with follows.
Step 2 (variable parameters). The previous derivation holds when the parameters are constant. When they vary smoothly, we expand them about and repeat the previous computations. Notice that due to the choice of , only the values of and in a narrow band about are relevant for the computation of the shape derivative. For the topological derivative, we set and an analogous expression for . Then, we check that the limits (90)-(91) also hold with and . For continuous and , we use a regularizing approximation, see  for details.
Step 3 (topological derivative with inclusions ). We have to redo the previous steps taking as objects for , and when depending on whether we are computing the derivative inside or outside ; see  for details.
The authors are partially supported by the Spanish Government Research Project TRA2010-18054 and the Spanish Ministerio de Economia y Competitividad Grants no. FIS2011-28838-C02-02 and no. FIS2010-22438-E.
- A. Liseno and R. Pierre, “Imaging of voids by means of a physical-optics-based shape-reconstrution algorithm,” Journal of the Optical Society of America A, vol. 21, no. 6, pp. 968–974, 2004.
- A. J. Devaney, “Geophysical diffraction tomography,” IEEE Transactions on Geoscience and Remote Sensing, vol. 22, no. 1, pp. 3–13, 1984.
- D. Colton, K. Giebermann, and P. Monk, “A regularized sampling method for solving three-dimensional inverse scattering problems,” SIAM Journal on Scientific Computing, vol. 21, no. 6, pp. 2316–2330, 2000.
- D. Colton and A. Kirsch, “A simple method for solving inverse scattering problems in the resonance region,” Inverse Problems, vol. 12, no. 4, pp. 383–393, 1996.
- F. Natterer and F. Wübbeling, “A propagation-backpropagation method for ultrasound tomography,” Inverse Problems, vol. 11, no. 6, pp. 1225–1232, 1995.
- R. E. Kleinman and P. M. van den Berg, “A modified gradient method for two-dimensional problems in tomography,” Journal of Computational and Applied Mathematics, vol. 42, no. 1, pp. 17–35, 1992.
- F. Hettlich, “Fréchet derivatives in inverse obstacle scattering,” Inverse Problems, vol. 11, no. 2, pp. 371–382, 1995.
- A. Kirsch, “The domain derivative and two applications in inverse scattering theory,” Inverse Problems, vol. 9, no. 1, pp. 81–96, 1993.
- M. Masmoudi, Outils pour la conception optimale des formes [Thése d'Etat en Sciences Mathématiques], Université de Nice, 1987.
- R. Potthast, “Fréchet differentiability of the solution to the acoustic Neumann scattering problem with respect to the domain,” Journal of Inverse and Ill-Posed Problems, vol. 4, no. 1, pp. 67–84, 1996.
- A. Litman, D. Lesselier, and F. Santosa, “Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set,” Inverse Problems, vol. 14, no. 3, pp. 685–706, 1998.
- F. Santosa, “A level-set approach for inverse problems involving obstacles,” Optimisation et Calcul des Variations, vol. 1, pp. 17–33, 1996.
- G. R. Feijoo, “A new method in inverse scattering based on the topological derivative,” Inverse Problems, vol. 20, no. 6, pp. 1819–1840, 2004.
- M. Bonnet and B. B. Guzina, “Sounding of finite solid bodies by way of topological derivative,” International Journal for Numerical Methods in Engineering, vol. 61, no. 13, pp. 2344–2373, 2004.
- B. B. Guzina and M. Bonnet, “Small-inclusion asymptotic of misfit functionals for inverse problems in acoustics,” Inverse Problems, vol. 22, no. 5, pp. 1761–1785, 2006.
- B. Samet, S. Amstutz, and M. Masmoudi, “The topological asymptotic for the Helmholtz equation,” SIAM Journal on Control and Optimization, vol. 42, no. 5, pp. 1523–1544, 2003.
- M. Bonnet and A. Constantinescu, “Inverse problems in elasticity,” Inverse Problems, vol. 21, no. 2, pp. R1–R50, 2005.
- O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” Inverse Problems, vol. 22, no. 4, pp. R67–R131, 2006.
- M. Burger, B. Hackl, and W. Ring, “Incorporating topological derivatives into level set methods,” Journal of Computational Physics, vol. 194, no. 1, pp. 344–362, 2004.
- S. Garreau, P. Guillaume, and M. Masmoudi, “The topological asymptotic for PDE systems: the elasticity case,” SIAM Journal on Control and Optimization, vol. 39, no. 6, pp. 1756–1778, 2001.
- A. Carpio and M.-L. Rapún, “Solving inhomogeneous inverse problems by topological derivative methods,” Inverse Problems, vol. 24, no. 4, Article ID 045014, 2008.
- A. Carpio and M.-L. Rapún, “An iterative method for parameter identification and shape reconstruction,” Inverse Problems in Science and Engineering, vol. 18, no. 1, pp. 35–50, 2010.
- A. Carpio and M.-L. Rapún, “Hybrid topological derivative and gradient-based methods for electrical impedance tomography,” Inverse Problems, vol. 28, no. 9, Article ID 095010, 2012.
- H. Ammari, J. Garnier, V. Jugnon, and H. Kang, “Direct reconstruction methods in ultrasound imaging of small anomalies,” in Mathematical Modelling in Biomedical Imaging II, vol. 2035 of Lecture Notes in Mathematics, Springer, Berlin, Germany, 2012.
- H. Ammari, J. Garnier, V. Jugnon, and H. Kang, “Stability and resolution analysis for a topological derivative based imaging functional,” SIAM Journal on Control and Optimization, vol. 50, no. 1, pp. 48–76, 2012.
- L. Borcea, G. Papanicolaou, C. Tsogka, and J. Berryman, “Imaging and time reversal in random media,” Inverse Problems, vol. 18, no. 5, pp. 1247–1279, 2002.
- A. Carpio and M.-L. Rapún, “Domain reconstruction using photothermal techniques,” Journal of Computational Physics, vol. 227, no. 17, pp. 8083–8106, 2008.
- L. Borcea, “Electrical impedance tomography,” Inverse Problems, vol. 18, no. 6, pp. R99–R136, 2002.
- M. Cheney, D. Isaacson, and J. C. Newell, “Electrical impedance tomography,” SIAM Review, vol. 41, no. 1, pp. 85–101, 1999.
- M. Brühl and M. Hanke, “Numerical implementation of two noniterative methods for locating inclusions by impedance tomography,” Inverse Problems, vol. 16, no. 4, pp. 1029–1042, 2000.
- E. T. Chung, T. F. Chan, and X.-C. Tai, “Electrical impedance tomography using level set representation and total variational regularization,” Journal of Computational Physics, vol. 205, no. 1, pp. 357–372, 2005.
- M. Hintermüller and A. Laurain, “Electrical impedance tomography: from topology to shape,” Control and Cybernetics, vol. 37, no. 4, pp. 913–933, 2008.
- M. Hintermüller, A. Laurain, and A. A. Novotny, “Second-order topological expansion for electrical impedance tomography,” Advances in Computational Mathematics, vol. 36, no. 2, pp. 235–265, 2012.
- B. Hofmann, “Approximation of the inverse electrical impedance tomography problem by an inverse transmission problem,” Inverse Problems, vol. 14, no. 5, pp. 1171–1187, 1998.
- K. Ito, K. Kunisch, and Z. Li, “Level-set function approach to an inverse interface problem,” Inverse Problems, vol. 17, no. 5, pp. 1225–1242, 2001.
- M. Hanke and M. Brühl, “Recent progress in electrical impedance tomography,” Inverse Problems, vol. 19, no. 6, pp. S65–S90, 2003.
- L. Borcea, “A nonlinear multigrid for imaging electrical conductivity and permittivity at low frequency,” Inverse Problems, vol. 17, no. 2, pp. 329–359, 2001.
- K. A. Dines and R. J. Lytle, “Analysis of electrical conductivity imaging,” Geophysics, vol. 46, no. 7, pp. 1025–1036, 1981.
- T. J. Yorkey, J. G. Webster, and W. J. Tompkins, “Comparing reconstruction algorithms for electrical impedance tomography,” IEEE Transactions on Biomedical Engineering, vol. 34, no. 11, pp. 843–852, 1987.
- B. B. Guzina and I. Chikichev, “From imaging to material identification: a generalized concept of topological sensitivity,” Journal of the Mechanics and Physics of Solids, vol. 55, no. 2, pp. 245–279, 2007.
- M. Abramovitz and I. A. Stegun, Handbook of Mathematical Functions, Dover, New York, NY, USA, 1972.
- J. Sokołowski and A. Żochowski, “On the topological derivative in shape optimization,” SIAM Journal on Control and Optimization, vol. 37, no. 4, pp. 1251–1272, 1999.
- A. Carpio and M. L. Rapún, “Topological derivatives for shape reconstruction,” in Inverse Problems and Imaging, vol. 1943 of Lecture Notes in Mathematics, pp. 85–133, Springer, Berlin, Germany, 2008.
- G. R. Feijóo, A. A. Oberai, and P. M. Pinsky, “An application of shape optimization in the solution of inverse acoustic scattering problems,” Inverse Problems, vol. 20, no. 1, pp. 199–228, 2004.
- M.-L. Rapún and F.-J. Sayas, “Indirect methods with Brakhage-Werner potentials for Helmholtz transmission problems,” in Numerical Mathematics and Advanced Applications: Proceedings of ENUMATH 2005, the 6th European Conference on Numerical Mathematics and Advanced Applications Santiago de Compostela, Spain, July 2005, pp. 1146–1154, Springer, Berlin, Germany, 2006.
- M.-L. Rapún and F.-J. Sayas, “A mixed-FEM and BEM coupling for the approximation of the scattering of thermal waves in locally non-homogeneous media,” ESAIM: Mathematical Modelling and Numerical Analysis, vol. 40, no. 5, pp. 871–896, 2006.
- A. Carpio, B. T. Johansson, and M.-L. Rapún, “Determining planar multiple sound-soft obstacles from scattered acoustic fields,” Journal of Mathematical Imaging and Vision, vol. 36, no. 2, pp. 185–199, 2010.
- D. Colton and B. D. Sleeman, “Uniqueness theorems for the inverse problem of acoustic scattering,” IMA Journal of Applied Mathematics, vol. 31, no. 3, pp. 253–259, 1983.
- T. Johansson and B. D. Sleeman, “Reconstruction of an acoustically sound-soft obstacle from one incident field and the far-field pattern,” IMA Journal of Applied Mathematics, vol. 72, no. 1, pp. 96–112, 2007.
- M. Cheney, D. Isaacson, and E. Somersalo, “Existence and uniqueness for electrode models for electric current computed tomography,” SIAM Journal on Applied Mathematics, vol. 52, no. 4, pp. 1023–1040, 1992.
- A. Lechleiter and A. Rieder, “Newton regularizations for impedance tomography: a numerical study,” Inverse Problems, vol. 22, no. 6, pp. 1967–1987, 2006.
- D. C. Dobson, “Convergence of a reconstruction method for the inverse conductivity problem,” SIAM Journal on Applied Mathematics, vol. 52, no. 2, pp. 442–458, 1992.