On the identification of piecewise constant coefficients in optical diffusion tomography by level set
aa r X i v : . [ m a t h . NA ] D ec On the identification of piecewise constant coefficients inoptical diffusion tomography by level set
J. P. Agnelli † A. De Cezaro ‡ A. Leit˜ao § M. Marques Alves ¶ December 23, 2020
Abstract
In this paper, we propose a level set regularization approach combined with a split strategyfor the simultaneous identification of piecewise constant diffusion and absorption coefficientsfrom a finite set of optical tomography data (Neumann-to-Dirichlet data). This problem isa high nonlinear inverse problem combining together the exponential and mildly ill-posedness ofdiffusion and absorption coefficients, respectively. We prove that the parameter-to-measurementmap satisfies sufficient conditions (continuity in the L topology) to guarantee regularizationproperties of the proposed level set approach. On the other hand, numerical tests consideringdifferent configurations bring new ideas on how to propose a convergent split strategy for thesimultaneous identification of the coefficients. The behavior and performance of the proposednumerical strategy is illustrated with some numerical examples.2000 Mathematics Subject Classification: 49N45, 65N21, 74J25.Key words: Optical Tomography, Parameter Identification, Level Set Regularization, Numer-ical Strategy . Optical tomography has demonstrated to be a powerful technique to obtain relevant physiologicalinformation of tissues in a non-invasive manner. The technique relies on the object under studybeing at least partially light-transmitting or translucent, so it works best on soft tissues such asbreast and brain tissue [24, 19]. By monitoring spatial-temporal variations in the light absorptionand scattering properties of tissue, regional variations in hemoglobin concentration or blood oxygen † FaMAF-CIEM, Universidad Nacional de C´ordoba, Medina Allende s/n 5000, C´ordoba, Argentina.( [email protected] ). ‡ Institute of Mathematics Statistics and Physics, Federal University of Rio Grande, Av. Italia km 8, 96201-900Rio Grande, Brazil ( [email protected] ). § Department of Mathematics, Federal University of St. Catarina, P.O. Box 476, 88040-900 Florian´opolis, Brazil( [email protected] ). ¶ Department of Mathematics, Federal University of St. Catarina, 88040-900 Florian´opolis, Brazil( [email protected] ). u , the equation to consider isthe following: −∇ · ( a ( x ) ∇ u ) + c ( x ) u = 0 in Ω (1) a ( x ) ∂u∂ν = g on Γ , (2)where Ω ⊂ R N , N ∈ { , , } , is open, bounded and connected with Lipschitz boundary denoted byΓ, the diffusion and absorption coefficients a ( x ) and c ( x ), respectively, are measurable real-valuedfunctions and ν is the outer-pointing normal vector. Moreover, g ∈ H − / (Γ) is the Neumannboundary data. Such boundary condition can be interpreted as the exitance on Γ.It is worth mentioning that equation (2) is a simplified way of modeling light fluence boundarycondition in diffuse optical tomography, since a more realistic description is to consider a Robinboundary condition [28, 27]. However, we believe that our simplified boundary model alreadycontains the essential aspects for the theoretical study that we present in this work. Further,this setting agrees with the uniqueness identification result derived by Harrach in [18]. For moredetails about boundary conditions in light propagation models we recommend to consult [28, 27]and references therein.Since the optical properties within tissue are determined by the values of the diffusion andabsorption coefficients, the problem of interest in DOT is the simultaneous identification of bothcoefficients from measurements of near-infrared diffusive light along the tissue boundary.In this contribution, we proposed a level set regularization approach [15, 10, 7, 8] combinedwith a split strategy for the simultaneous identification of piecewise constant diffusion a ( x ) andabsorption c ( x ) coefficients in (1)–(2), from a finite set of available measurements of the photondensity h := u | Γ , corresponding to inputs g ∈ H − / (Γ) in (1)–(2). Related works:
In [22], a Levenberg–Marquardt method for recovering internal boundaries ofpiecewise constant coefficients of an elliptic PDE as (1) was implemented. The proposed methodis based on the series expansion approximation of the smooth boundaries and on the finite elementmethod. However, in [22], there is not a theoretical result that guarantees regularizing properties ofthe iterated approximated solution. Indeed, as far as the authors are aware, there is not theoreticalregularization approaches in the literature for recovering the pair of coefficients ( a, c ) in (1) fromboundary data.In [31], the authors did a carefully designed experiment aimed to provide solid evidence that bothabsorption and scattering images of a heterogeneous scattering media can be reconstructed inde-pendently from diffuse optical tomography data. The authors also discuss the absorption scatteringcross-talk issue.Although it is well known that the identification of a and c simultaneously is not possiblein a general case [2], recently B. Harrach [18] obtained a uniqueness result for the simultaneousrecovery of a and c in (1)–(2) assuming that a ≥ a > c ∈ L ∞ + 1 is The subscript ′ + ′ denotes positive essential infima. a∂ ν u | ˜Γ and u | ˜Γ on anarbitrarily small open set ˜Γ of the boundary Γ. The difference between the work of B. Harrach [18]and our work is that here we are considering a more practical approach: we only have access to afinite number of Neumann-Dirichlet pairs.We also remark that the quantitative photoacoustic tomography problem (QPAT), in the diffu-sive approach, also aims to simultaneous recover ( a, c ) of an elliptic boundary valued problem. Seefor example [29] and references therein. However, in the QPAT situation the solution of the “firstinverse problem” generate internal data for the reconstruction. In this sense, the QPAT problem isvery different to the identification problem that we are facing here. Novelties:
The novelties of this contribution are divided as follow: • We prove the continuity of the parameter-to-measurement (forward) map F (defined in (3))in the [ L (Ω)] topology. It is done in Theorem 3 of Section 2, and it is possible thanks to ageneralization of Meyers’ Theorem [23] that prove the regularity of the solution of (1)–(2) in W ,p (Ω) for some p > H (Ω) functions concatenated with arestriction of the search space using nonlinear constraints. Such approach allows us to enforcethe desired additional properties on the pair of parameters ( a, c ) (namely: ( a, c ) is a pair ofpiecewise constant functions describing the high diffusion and absorption contrast betweenthe optical properties of the object) that are not smooth.Given the continuity of F in the [ L (Ω)] topology, it is now a standard result to prove thatthe level set approach is a regularization method as in the classical theory of regularization [7,8, 10, 11]. Therefore, we only point out the convergence and stability results without a formalproof in Subsection 3.2. • Another contribution of the proposed level set approach is related to the numerical implemen-tation presented in Section 5. It is worth to remind the reader that we aim to simultaneousreconstruct the pair ( a, c ) of piecewise constant functions from a finite set of optical measure-ments. With this aim, we first run several numerical experiments in order to recover theabsorption coefficient c , based on either total or partial knowledge of a . From this experi-ments we observed that the level set method for identifying c performs well, even if a goodapproximation of the exact value of a is not known. This is presented in Subsection 5.1. Afterthat, in Subsection 5.2 we run another set of experiments but now concerning the identifi-cation of the diffusion coefficient a , based on either total or partial knowledge of c . In thiscase, we observed that the level set method for identifying a performs well if a good approx-imation of the exact value of c is available, but may generate a sequence a k that does notapproximate the exact a if the initial guess for the coefficient c is far from its exact value.Such features of the identification problems suggested one of the main results related to thenumerical perspective presented in Subsection 5.3. Given a initial guess ( a , c ), we adoptedthe strategy to “freeze” the coefficient a k = a during the first iterations, and to iterate the3lgorithm only with respect to the coefficient c . We follow this strategy until the iteratedsequence c k stagnates. Then, we freeze the absorption coefficient c = c k and iterate the algo-rithm only with respect to a until the iterated sequence a k stagnates. Finally, we iterate bothcoefficients simultaneously. This numerical strategy has not only demonstrated that givesvery good results but also reduces significantly the computational effort.This article is organized as follows. In Section 2 we first introduce the parameter-to-measurement(forward) map F and after that, in Theorem 3, the continuity of this forward map is demonstrated.Then, in Section 3 we present the level set approach, we introduce the concept of generalizedminimizers for an appropriate energy functional and we establish the regularization properties. Inother words, we prove the well-posedness result and also convergence results for exact and noisydata. In Section 4, we introduce a smooth functional that is used in the numerical examples. Weprove that the minimizers of such functional converge to a minimizer of the early energy functionalin appropriated topologies. Section 5 is devoted to numerical experiments and a split strategy isdeveloped. We end this contribution in Section 6 with some conclusions and further developments.In the Appendix, we give a proof for a generalization of a Meyers’ type theorem (see Theorem 1)about the regularity of the solution of (1)–(2) that are used in the proof of the continuity of theforward map F . General Notation.
We denote by R N , N ≥
2, the N -dimensional Euclidean space endowedwith the usual scalar product x · y = P Ni =1 x i y i and norm | x | = √ x · x , where x = ( x i ) Ni =1 and y = ( y i ) Ni =1 . Given two normed vector spaces ( X , k·k X ) and ( Y , k·k Y ) we always consider the productspace X × Y endowed with the product topology generated by the norm k ( x, y ) k := k x k X + k y k Y (or the equivalent norms (cid:0) k x k X + k y k Y (cid:1) / or max {k x k X , k y k Y } ), where ( x, y ) ∈ X × Y . We alsouse the short notation X = X × X . We start this section assuming that coefficient a ( x ) is known for all x ∈ Γ. Then, for each input g ∈ H − / (Γ) in (1)–(2), we define the parameter-to-measurement (forward) map F := F g : D ( F ) ⊂ L (Ω) × L (Ω) → H / (Γ) (3)( a, c ) h := u | Γ , where u = u ( g ) is the unique solution of (1)–(2) given the boundary data g and the pair ( a, c ) inthe parameter space D ( F ) defined as: Definition 1.
Denote by D ( F ) the set of pairs of L (Ω) functions ( a, c ) on Ω satisfying the followingcondition: < a ≤ a ( x ) ≤ a, < c ≤ c ( x ) ≤ c ∀ x a.e. in Ω , (4) where a, a , c and c are known positive real numbers.
4e now make some comments about Definition 1 and the definition of the forward map F . First,it is easy to check that D ( F ) is a convex subset of [ L (Ω)] . Second, the forward map F is well-defined because for each ( a, c ) ∈ D ( F ) there exists a unique solution u ∈ H (Ω) of (1)–(2) (see [6]).Third, since D ( F ) depends on the scalars a, a , c and c , it turns out that F also depends on thelatter scalars. However, we are assuming that the scalars are known, fixed and independent of eachgiven Neumann data g in (2). Fourth, we are not assuming any smoothness condition on the pair( a, c ) ∈ D ( F ). In particular the latter fact allow us to consider solutions of (1)–(2) correspondingto piecewise constant coefficients.This section is devoted to prove the continuity of the forward map F in the [ L (Ω)] topology.In order to make such proof easier to understand we will consider the parameter-to-solution map G := G g : D ( F ) ⊂ [ L (Ω)] −→ H (Ω)( a, c ) G g ( a, c ) := u , (5)where u = u ( g ) ∈ H (Ω) is the unique solution of (1)–(2) for each input data g ∈ H − / (Γ) andparameters ( a, c ) ∈ D ( F ). Moreover, we will use the fact that any solution of (1)–(2) satisfies thefollowing weak formulation [6]: Z Ω a ∇ u · ∇ ϕ dx + Z Ω cuϕ dx = Z Γ gϕ dσ ∀ ϕ ∈ H (Ω) . (6) Remark 1.
Given the definition of the forward map F in (3) , and using the map defined in (5) ,we have that F can be written as F = γ ◦ G, (7) where γ : H (Ω) → H / (Γ) is the trace operator of order zero [6]. Since the operator γ is linearand continuous [6], the continuity of F follows from the continuity of G . In order to prove the continuity of the operator G defined in (5) in the desired topology, we willuse the following generalization of Meyers’ Theorem [16] on the regularity of the solution of (1)–(2).The proof of Theorem 1 is presented in the Appendix. Theorem 1 (Generalized Meyers’ Theorem) . Let Ω ⊂ R N , N ∈ { , , } , be a connected boundedopen set with a Lipschitz boundary Γ and let ( a, c ) ∈ D ( F ) . Then, there exists a real number p M > (depending only on Ω , a, a, c and c ) such that the following condition hold for every p ∈ (2 , p M ) : If g ∈ W − (1 /q ) ,q (Γ) ′ , where q := p/ ( p − , then the unique solution u of (1) – (2) belongs to W ,p (Ω) . Next, we present a lemma that is used in the main theorem of this section.
Lemma 2.
Let h be a measurable function such that | h ( x ) | ≤ M for all x a.e. in Ω , for someconstant M > . Then, h ∈ L s (Ω) for all ≤ s < ∞ and k h k L s (Ω) ≤ M ( s − /s k h k /sL (Ω) . Proof.
Note that k h k sL s (Ω) = Z Ω | h ( x ) || h ( x ) | s − dx ≤ M s − k h k L (Ω) , which readily implies the desired result. 5n the next theorem we prove the continuity of the parameter-to-solution map G in the [ L (Ω)] topology. Then, by Remark 1 we obtain the desired continuity of the parameter-to-measurement map F . Theorem 3.
Let p ∈ (2 , p M ) , where p M > is given by Theorem 1, and let q := p/ ( p − . Then,for any g ∈ W − (1 /q ) ,q (Γ) ′ the operator G defined in (5) is continuous in the [ L (Ω)] topology.As a consequence, for any g ∈ W − (1 /q ) ,q (Γ) ′ , the forward map F defined in (3) is also continuousin the [ L (Ω)] topology.Proof. Let g ∈ W − (1 /q ) ,q (Γ) ′ and consider the corresponding solutions u ′ = u ( a ′ , c ′ ) and u = u ( a, c )of (1)–(2) with parameters ( a ′ , c ′ ) , ( a, c ) ∈ D ( F ), respectively.Since ( a ′ , c ′ , u ′ ) and ( a, c, u ) satisfy the identity (6) for all ϕ ∈ H (Ω) we have Z Ω ( a ∇ u − a ′ ∇ u ′ ) · ∇ ϕ dx + Z Ω ( cu − c ′ u ′ ) ϕ dx = 0 . (8)Defining w := u − u ′ ∈ H (Ω) and using (8) with ϕ = w we obtain (after some algebraic manipula-tions) Z Ω (cid:0) a − a ′ (cid:1) ∇ u ′ · ∇ w dx + Z Ω a ∇ w · ∇ w dx + Z Ω (cid:0) c − c ′ (cid:1) u ′ w dx + Z Ω cww dx = 0 , which in turn is equivalent to Z Ω a ( x ) |∇ w | dx + Z Ω c ( x ) | w | dx = Z Ω (cid:0) a ′ − a (cid:1) ∇ u ′ · ∇ w dx + Z Ω (cid:0) c ′ − c (cid:1) u ′ w dx. (9)In view of Theorem 1 (for ( a ′ , c ′ )) we have u ′ ∈ W ,p (Ω). Thus, defining s := 2 p/ ( p − /s + 1 /p + 1 / { a, c }k w k H ≤ k a ′ − a k L s k∇ u ′ k L p k∇ w k L + k c ′ − c k L s k u ′ k L p k w k L ≤ (cid:0) k a ′ − a k L s k∇ u ′ k L p + k c ′ − c k L s k u ′ k L p (cid:1) k w k H ≤ { a − a, c − c } ) ( s − /s k u ′ k W ,p (cid:0) k a ′ − a k L + k c ′ − c k L (cid:1) /s k w k H . The latter inequality combined with the facts that G g ( a ′ , c ′ ) = u ′ , G g ( a, c ) = u (see (5)) and w = u − u ′ give k G g ( a, c ) − G g ( a ′ , c ′ ) k H ≤ f M k u ′ k W ,p (cid:0) k a − a ′ k L + k c − c ′ k L (cid:1) /s , (10)which proves the continuity of G g in the [ L (Ω)] topology, where f M := 2 (max { a − a, c − c } ) ( s − /s min { a, c } .The last statement of the theorem now follows easily from the first one and Remark 1.We now make a few comments about Theorem 3. First, according to Theorem 1, the real number p M > a , a , c and c . Second, since q <
2, it followsthat W − (1 /q ) ,q (Γ) ′ ⊂ H − / (Γ). As a consequence of the latter inclusion, we have that the conditionon g required in Theorem 3 is stronger than the usual inclusion g ∈ H − / (Γ). Third, condition (10)gives that both operators G and F are (locally) H¨older continuous in the [ L (Ω)] topology.6 The level set framework with a finite number of experi-ments
It is already known that in diffuse optical tomography the full Neumann-to-Dirichlet map (equiv-alently, the boundary data h corresponding to the boundary condition a ∂u∂ν on Γ) is required toobtain uniqueness of the parameters ( a, c ) in (1)–(2) [18]. However, in real applications, only afinite number of observations/measurements are available. Therefore, in this work we consider thatwe only have access to a quantity ℓ ∈ N of well-placed experiments. In other words, the inverseproblem we tackle consist in given a finite number of inputs g m = a ∂u m ∂ν | Γ and corresponding data h m = u m | Γ , reconstruct simultaneously the diffusion and absorption coefficients ( a, c ). As indicatedpreviously, the photon density u m satisfies ∇ · ( a ∇ u m ) + c u m = 0 in Ω ,a ∂u m ∂ν = g m on Γ , m = 1 , . . . , ℓ. This problem is known in the literature as the inverse problem for the Neumann-to-Dirichlet operatorwith a finite number of experiments. In this context, the identification problem can be written interms of the system of nonlinear equations F m ( a, c ) = h m , m = 1 , . . . , ℓ , (11)where F m := F g m is defined as in (3), for each m ∈ { , . . . , ℓ } .Moreover, given the nature of the measurements, we can not expect that exact data h m ∈ H / (Γ)are available. Instead, one disposes only an approximate measured data h δm ∈ L (Γ) satisfying (cid:13)(cid:13) h δm − h m (cid:13)(cid:13) L (Γ) ≤ δ , for m = 1 , . . . , ℓ (12)where δ > Remark 2.
From Theorem 3, we know that each forward map F m in (11) is continuous in the [ L (Ω)] topology. In contrast with the previous section, from now on we consider that the pair of parameters ( a, c ) arepiecewise constant function assuming two distinct values, i.e. a ( x ) ∈ { a , a } and c ( x ) ∈ { c , c } a.e. in Ω ⊂ R N , but we still consider ( a, c ) ∈ D ( F ). Hence, one can assume the existence of openand mensurable sets A ⊂⊂ Ω and C ⊂⊂ Ω with H ( ∂ A ) < ∞ and H ( ∂ C ) < ∞ , and suchthat a ( x ) = a if x ∈ A and a ( x ) = a if x ∈ A := Ω − A ; c ( x ) = c if x ∈ C and c ( x ) = c if x ∈ C := Ω − C . Consequently, the pair of parameters can be modeled as( a ( x ) , c ( x )) = ( a + ( a − a ) χ A ( x ) , c + ( c − c ) χ C ( x )) , where χ S is the indicator function of the set S . Here H ( S ) denotes the one-dimensional Hausdorff-measure of the set S . evel set framework: In order to model the space of admissible parameters, that is the pairof piecewise constant functions ( a, c ), we consider a standard level set (SLS) approach proposed in[15, 7, 8, 10]. In particular, our analysis of a level set approach for piecewise constant parametersfollows essentially from the techniques derived in [8]. We notice that many other level set approachesare known in the literature, see for instance [11, 12, 5, 22, 32, 30]. For the case where not only thediscontinuities but also the values of a and c are unknown then one can use the ideas of the level setapproach presented in [7]. Recently, in [9, 10, 11], piecewise constant level set approaches (PCLS)were derived for identification of piecewise constant parameters. The PCLS approach consists inintroducing constraints in the admissible class of level set functions in order to enforce these levelset functions to become piecewise constant. In this context, we do not need to introduce theHeaviside projector H (see below) to model the parameter space. However, the introduction ofconstraints imply different difficulties in the level set regularization analysis [9, 10]. Advantagesand disadvantages of SLS and PCLS approaches were discussed in [9, 10].According to the SLS representation strategy, level set functions φ a , φ c : Ω → R , in H (Ω), arechosen in such a way that its zero level set Γ φ a := { x ∈ Ω ; φ a ( x ) = 0 } and Γ φ c := { x ∈ Ω ; φ c ( x ) =0 } define connected curves within Ω and the discontinuities of the parameters ( a, c ) are located“along” Γ φ a and Γ φ c , respectively.Introducing the Heaviside projector H ( t ) := ( , if t > , if t ≤ , the diffusion and absorption parameters can be written as( a, c ) = (cid:0) a + ( a − a ) H ( φ a ) , c + ( c − c ) H ( φ c ) (cid:1) =: P ( φ a , φ c ) . (13)Notice that ( a ( x ) , c ( x )) = ( a i , c j ), x ∈ A i ∩ C j for i, j ∈ { , } , where the sets A i and C j aredefined by A = { x ∈ Ω : φ a ( x ) ≥ } , A = { x ∈ Ω : φ a ( x ) < } , C = { x ∈ Ω : φ c ( x ) ≥ } and C = { x ∈ Ω : φ c ( x ) < } . Thus, the operator P establishes a straightforward relation betweenthe level sets of φ a and φ c and the sets A i and C j that characterize the coefficients ( a, c ).As already observed in [8], the operator H maps H (Ω) into the space V , := { z ∈ L ∞ (Ω) | z = χ S , S ⊂ Ω measurable , H ( ∂S ) < ∞} . Therefore, the operator P in (13) maps H (Ω) × H (Ω) into the admissible class V defined by V := { ( z , z ) ∈ [ L ∞ (Ω)] | ( z , z ) = ( a + ( a − a ) χ A , c + ( c − c ) χ C ) , for some A , C ⊂ Ω } . Within this framework, the inverse problem in (11), with data given as in (12), can be writtenin the form of the operator equation F m ( P ( φ a , φ c )) = h δm m = 1 , . . . , ℓ. (14)Let us make the following general assumption: (A1) Equation (11) has a solution, i.e. there exists ( a ∗ , c ∗ ) ∈ L ∞ (Ω) × L ∞ (Ω) satisfying F ( a ∗ , c ∗ ) = h m , for m = 1 , . . . , ℓ . Moreover, there exists a pair of functions ( φ a ∗ , φ c ∗ ) ∈ [ H (Ω)] satisfying P ( φ a ∗ , φ c ∗ ) = ( a ∗ , c ∗ ), with |∇ φ a ∗ | 6 = 0 and |∇ φ c ∗ | 6 = 0 in a neighborhoodof { φ a ∗ = 0 } and { φ c ∗ = 0 } respectively and such that H ( φ a ∗ ) = z a = χ A ∈ L ∞ (Ω), H ( φ c ∗ ) = z c = χ C ∈ L ∞ (Ω). 8 .2 Level set regularization Since the unknown coefficients ( a, c ) are piecewise constant functions, a natural alternative to obtainstable solutions of the operator equation (11) is to use a least-square approach combined with atotal variation regularization. This corresponds to a Tikhonov-type regularization [7, 8, 10]. Withinthe level set framework presented above, the Tikhonov-type regularization approach for obtaininga regularized solution to the operator equation (14) is based on the minimization of the energyfunctional F α ( φ a , φ c ) := ℓ X m =1 k F m ( P ( φ a , φ c )) − h δm k L (Γ) + αR ( φ a , φ c ) , (15)where R ( φ a , φ c ) = (cid:16) β a | H ( φ a ) | BV (Ω) + β c | H ( φ c ) | BV (Ω) + k φ a − φ a k H (Ω) + k φ c − φ c k H (Ω) (cid:17) ,α > β j play the role of scaling factors.This approach is based on TV- H penalization. The H –terms act as a regularization for thelevel set functions on the space H (Ω) whereas the BV (Ω)-seminorm terms are well known forpenalizing the length of the Hausdorff measure of the boundary of the sets { x ∈ Ω : φ a ( x ) > } , { x ∈ Ω : φ c ( x ) > } (see [14]).In general, variational minimization techniques involve compact embedding arguments on theset of admissible minimizers and continuity of the operator in such set to guarantee the existence ofminimizers. The Tikhonov functional in (15) does not allow such characteristic, since the Heavisideoperator H and consequently the operator P are discontinuous. Therefore, given a minimizingsequence ( φ ak , φ ck ) for F α we cannot prove existence of a (weak-*) convergent subsequence. Conse-quently, we cannot guarantee the existence of a minimizer in [ H (Ω)] . In other words, the graphof F α is not closed in the desired topology.To overcome this difficulty in [15, 8] was introduced the concept of generalized minimizers werethe graph of F α becomes closed. It allow us to guarantee the existence of minimizers of the Tikhonovfunctional (15). For sake of completeness, we present the concept of generalized minimizers below. The concept of generalized minimizers:
For each ε >
0, we define the smooth approximationto H given by: H ε ( t ) := (cid:26) t/ε for t ∈ [ − ε, H ( t ) for t ∈ R \ [ − ε, P ε ( φ a , φ c ) := ( a H ε ( φ a ) + a (1 − H ε ( φ a )) , c H ε ( φ c ) + c (1 − H ε ( φ c ))) . (16) Definition 2.
Let the operators H , P , H ε and P ε be defined as above.(a) A vector ( z , z , φ a , φ c ) ∈ [ L ∞ (Ω)] × [ H (Ω)] is called admissible when there existsequences { φ ak } and { φ ck } of H (Ω) -functions satisfying lim k →∞ k φ ak − φ a k L (Ω) = 0 , lim k →∞ k φ ck − φ c k L (Ω) = 0 , nd there exists a sequence { ε k } ∈ R + converging to zero such that lim k →∞ k H ε k ( φ ak ) − z k L (Ω) = 0 and lim k →∞ k H ε k ( φ ck ) − z k L (Ω) = 0 . (b) A generalized minimizer of the Tikhonov functional F α in (15) is considered to be anyadmissible vector ( z , z , φ a , φ c ) minimizing G α ( z , z , φ a , φ c ) := ℓ X m =1 k F m ( Q ( z , z )) − h δm k L (Ω) + αR ( z , z , φ a , φ c ) (17) over the set of admissible vectors, where Q : [ L ∞ (Ω)] ∋ ( z , z ) ( a z + a (1 − z ) , c z + c (1 − z )) ∈ [ L ∞ (Ω)] , and the functional R is defined by R ( z , z , φ a , φ c ) := ρ ( z , z , φ a , φ c ) , with ρ ( z , z , φ a , φ c ) := inf n lim inf k →∞ (cid:16) β a | H ε k ( φ ak ) | BV (Ω) + β c | H ε k ( φ ck ) | BV (Ω) + k φ ak − φ a k H (Ω) + k φ ck − φ c k H (Ω) (cid:17)o . Here the infimum is taken over all sequences { ε k } and { φ ak , φ ck } characterizing ( z , z , φ a , φ c ) as anadmissible vector. In this subsection we present the regularization properties of the proposed level set approach tothe inverse problem of identifying ( a, c ) in the diffuse optical tomography model (1)–(2). Since theresults follow straightforward arguments presented in [7, 8, 10] we do not present their proofs here.
Theorem 4.
The following assertions hold true.i) The functional G α in (15) attains minimizers on the set of admissible vectors.ii) [Convergence for exact data] Assume that we have exact data, i.e. h δ = h . For every α > denote by ( z α , z α , φ aα , φ cα ) a minimizer of G α on the set of admissible vectors. Then, for everysequence of positive numbers { α k } converging to zero there exists a subsequence, denoted again by { α k } , such that ( z α k , z α k , φ aα k , φ cα k ) is strongly convergent in [ L (Ω)] × [ L (Ω)] . Moreover, the limitis a solution of (11) .iii) [Convergence for noisy data] Let α = α ( δ ) be a function satisfying lim δ → α ( δ ) = 0 and lim δ → δ α ( δ ) − = 0 . Moreover, let { δ k } be a sequence of positive numbers converging to zeroand { h δ k } ∈ L (Γ) be corresponding noisy data satisfying (12) . Then, there exists a subsequence,denoted again by { δ k } , and a sequence { α k := α ( δ k ) } such that ( z α k , z α k , φ aα k , φ cα k ) converges in [ L (Ω)] × [ L (Ω)] to a solution of (14) .Proof. The proof follows the arguments presented in [8], Theorem 6, Theorem 8 and Theorem 9respectively and therefore is omitted. 10
Numerical realization
In this section we introduce the functional G ε,α , which can be used for the purpose of numericalimplementations. This functional is defined in such a way that it’s minimizers are “close” to thegeneralized minimizers of F α in a sense that will be clear later (see Proposition 5). For each ε > G ε,α ( φ a , φ c ) := ℓ X m =1 k F m ( P ε ( φ a , φ c )) − h δm k L (Γ) + αR ε ( φ a , φ c ) , (18)where R ε ( φ a , φ c ) := (cid:16) β a | H ε ( φ a ) | BV (Ω) + β c | H ε ( φ c ) | BV (Ω) + k φ a − φ a k H (Ω) + k φ c − φ c k H (Ω) (cid:17) . (19)The next result guarantees that for ε → G ε,α attains a minimizer. Moreover,the minimizers of G ε,α approximate a generalized minimizer of F α . Proposition 5. i) Given α , β j , ε > and φ a , φ c in H (Ω) , then the functional G ε,α in (18) attains a minimizer on [ H (Ω)] .ii) Let α , β j be given. For each ε > denote by ( φ aε,α , φ cε,α ) a minimizer of G ε,α . There exists asequence of positive numbers { ε k } converging to zero such that ( H ε k ( φ aε k ,α ) , H ε k ( φ cε k ,α ) , φ aε k ,α , φ cε k ,α ) converges strongly in [ L (Ω)] × [ L (Ω)] and the limit is a generalized minimizer of F α in the setof admissible vectors.Proof. The proof follows from Lemma 10 and Theorem 11 presented in [8]. Therefore, we do notpresent the details it in this paper.Proposition 5 justifies the use of functional G ε,α in order to obtain numerical approximations tothe generalized minimizers of F α . It is worth noticing that, differently from F α , the minimizers of G ε,α can be actually computed. In the next subsection we derive the first order optimality conditionsfor the functional G ε,α , which will allow us to compute the desired minimizers. G ε,α For the numerical purposes we have in mind, it is necessary to derive the first order optimalityconditions for a minimizer of the functional G ε,α . To this end, we consider G ε,α in (18) and we lookfor the Gˆateaux directional derivatives with respect to φ a , φ c . In order to simplify the presentation,we will assume that the values a , a , c , c are known. Since H ′ ε ( ϕ ) is self-adjoint, the optimalityconditions for a minimizer of the functional G ε,α can be written in the form of the system of equations α (∆ − I )( φ a − φ a ) = L aε,α ( φ a , φ c ) , α (∆ − I )( φ c − φ c ) = L cε,α ( φ a , φ c ) , in Ω (20a) ∂∂ν ( φ a − φ a ) = 0 , ∂∂ν ( φ c − φ c ) = 0 , on Γ (20b) Notice that H ′ ε ( t ) = ( /ε t ∈ ( − ε, else . ν ( x ) is the external unit normal vector at x ∈ Γ and L aε,α ( φ a , φ c ) = ( a − a ) H ′ ε ( φ a ) " l X m =1 (cid:18) ∂F m ( P ε ( φ a , φ c )) ∂φ a (cid:19) ∗ ( F m ( P ε ( φ a , φ c )) − h δm ) − αβ a (cid:20) H ′ ε ( φ a ) ∇· (cid:18) ∇ H ε ( φ a ) |∇ H ε ( φ a ) | (cid:19)(cid:21) (21a) L cε,α ( φ a , φ c ) = ( c − c ) H ′ ε ( φ c ) " l X m =1 (cid:18) ∂F m ( P ε ( φ a , φ c )) ∂φ c (cid:19) ∗ ( F m ( P ε ( φ a , φ c )) − h δm ) − αβ c (cid:20) H ′ ε ( φ c ) ∇· (cid:18) ∇ H ε ( φ c ) |∇ H ε ( φ c ) | (cid:19)(cid:21) . (21b)Note that, in order to implement the numerical algorithm for solving the optimality conditions,we need to calculate the adjoint of the derivatives ∂F m ∂φ a and ∂F m ∂φ c . Remark 3.
Given the level set functions φ a , φ c ∈ H (Ω) and inputs g m ∈ H / (Γ) for m = 1 , . . . , ℓ ,denote the residual r m := F m ( P ε ( φ a , φ c )) − h δm ∈ L (Γ) . Then, (cid:18) ∂F m ( P ε ( φ a , φ c )) ∂φ a (cid:19) ∗ r m = ∇ u m · ∇ w m (22) and (cid:18) ∂F m ( P ε ( φ a , φ c )) ∂φ c (cid:19) ∗ r m = − u m w m (23) where u m and w m are the unique solutions of the following elliptic boundary problems −∇ ( a ∇ u m ) + c u m = 0 , in Ω (24) a ∂u m ∂ν = g m , on Γ , −∇ ( a ∇ w m ) + c w m = 0 , in Ω (25) a ∂w m ∂ν = r m , on Γ , for m = 1 , . . . , ℓ respectively. We have already introduced all the ingredients necessary to implement an algorithm based on thelevel set regularization approach to solve the identification problem in diffuse optical tomography(see Table 1). The iterative algorithm consists in minimizing, for k ≥
1, the functional G ( k ) ε,α ( φ a , φ c ) := ℓ X m =1 k F m ( P ε ( φ a , φ c )) − h δm k L (Γ) + αR ( k ) ε ( φ a , φ c ) , (26)where R ( k ) ε is the functional R ε defined in (19) with φ j replaced by φ jk − . The minimizer of eachfunctional can be computed solving the formal optimality conditions (20) with φ j replaced by φ jk − .12 . Evaluate the residual [ r k,m ] ℓm =1 := [ F m ( P ε ( φ ak , φ ck )) − h δm ] ℓm =1 = [ u k,m | Γ − h δm ] ℓm =1 ,where [ u k,m ] ℓm =1 ∈ [ H (Ω)] ℓ and each function solves (24) . Evaluate h(cid:16) ∂F m ( P ε ( φ ak ,φ ck )) ∂φ ak (cid:17) ∗ r k,m i ℓm =1 = [ ∇ w k,m · ∇ u k,m ] ℓm =1 ∈ [ L (Ω)] ℓ , and h(cid:16) ∂F m ( P ε ( φ ak ,φ ck )) ∂φ ck (cid:17) ∗ r k,m i ℓm =1 = − [ w k,m u k,m ] ℓm =1 ∈ [ L (Ω)] ℓ , where [ u k,m ] ℓm =1 are thefunctions computed in Step and [ w k,m ] ℓm =1 ∈ [ H (Ω)] ℓ solve (25) . Calculate L aε,α ( φ ak , φ ck ) and L cε,α ( φ ak , φ ck ) given by equations (21a) and (21b) . Evaluate the updates δφ ak , δφ ck ∈ H (Ω) by solving (∆ − I ) δφ jk = L jε,α ( φ ak , φ ck ) , in Ω ; ∂δφ jk ∂ν = 0 , on Γ . Update the level set functions φ ak +1 = φ ak + α δφ ak , φ ck +1 = φ ck + α δφ ck . Table 1: An explicit algorithm based on the iterative regularization method for solving the identi-fication problem in diffuse optical tomography.Each iteration of the proposed algorithm consists in the next five steps: • In the first step the residual vector [ r k,m ] ℓm =1 ∈ [ L (Γ)] ℓ , corresponding to the k th iteration( φ ak , φ ck ), is evaluated. This requires the solution of ℓ elliptic BVP’s given by (24). • The second step consists in computing the adjoint of the partial derivatives of F m appliedto the residuals. This is done by solving ℓ elliptic BVP given by (25) to get the solutions[ w k,m ] ℓm =1 ∈ [ H (Ω)] ℓ and then computing the products given by Remark 3. • In the third step, the terms L aε,α ( φ ak , φ ck ) and L cε,α ( φ ak , φ ck ) given by equations (21a) and (21b)are calculated. • The fourth step consists in computing the updates δφ ak , δφ ck ∈ H (Ω) for the level-set functions φ a and φ c . This corresponds to solving two non-coupled elliptic BVP’s, namely (20a) and (20b). • Finally, update the level set functions and go to step 1 until a stopping criteria is reached.A similar algorithm was successfully implemented in [15, 8] to solve the inverse potential problemunder the framework of level sets and multiple level sets respectively. Regarding our coefficientidentification problem in diffuse optical tomography the algorithm outlined above also seems to beeffective (see next Section), but in this case, it has the disadvantage that in each iteration step onehas to solve 2 ℓ + 2 elliptic BVP’s. Then, if the number ℓ of experiments is large the computationalcost will be high. 13 Numerical Experiments
In this section we implement a numerical algorithm based on the level set approach derived inthe previous sections for identifying the coefficient pair ( a, c ) in (1)–(2). First, the identificationof the absorption coefficient c , based on either total or partial knowledge of a , is considered inSection 5.1. Then, the separate identification of the diffusion coefficient a , based on either total orpartial knowledge of c , is considered in Section 5.2. Finally, the simultaneous identification of thepair ( a, c ) is investigated in Section 5.3.In all the numerical experiments of this Section we considered Ω = (0 , × (0 ,
1) and four ( ℓ = 4)different inputs g m ∈ L (Γ) were applied as Neumann boundary conditions in order to compute thecorresponding Dirichlet data h m . Each of these functions is supported at one of the four sides of Γ,for instance g ( x ) = (cid:26) , if x ∈ ( , ) × { } , else , and g , g and g are defined in a similar way. All boundary value problems were solved using aGalerkin Finite Element method in an uniform grid with 50 nodes at each boundary side. We useda custom implementation using MATLAB.(a) (b) (c) (d)Figure 1: (a) Exact coefficients for the first experiments in Sections 5.1 and 5.2. (b)
Number of iterationsneeded for the identification of the absorption coefficient c ∗ , starting from distinct initial guesses c ( a ∗ isgiven; see Section 5.1). (c) Number of iterations needed for the identification of the diffusion coefficient a ∗ ,starting from distinct initial guesses a ( c ∗ is given; see Section 5.2). (d) Exact coefficients for the secondexperiments in Sections 5.1 and 5.2.
In what follows we consider the identification of the absorption coefficient c , based on either totalor partial knowledge of a . The values assumed for these coefficients were (see Figure 1): a ∗ ( x ) = (cid:26) , inside blue inclusion1 , elsewhere , c ∗ ( x ) = (cid:26) , inside red inclusion1 , elsewhere . As initial guess for the level set method we have chosen distinct piecewise constant functions c ,whose supports are shown in Figure 1 (b). It is worth noting that each c corresponds to a level14et function φ c ∈ H (Ω). In all cases the initial level set function φ c was a paraboloid but withdifferent minima.The constant values assumed by the exact solution c ∗ are supposed to be known, as well asthe exact diffusion coefficient a ∗ . Moreover, exact data was considered for the reconstruction (i.e., δ = 0) and we tested the iterative level set regularization without the penalizing term | H ε ( φ j ) | BV (Ω) ,i.e., β j = 0 (see [15, Remark 5.1] ).In this and in all the following computed experiments of this Section, we considered the operator P ε defined in (16) with ε = 1 /
10. This election was motivated by the fact that as ε increases, thesupports of the functions appearing on the right-hand side of (20a) and (20b) become larger (dueto the term H ε ). Consequently, the updates δφ ak , δφ ck ∈ H (Ω) given by these equations have largevalues. If ε becomes too large, the level set method becomes unstable. Therefore, ε was chosen tomatch the mesh size considered to solve the boundary problems.The inverse problem we tackle here reduces to a shape identification problem for the absorptioncoefficient.Notice that, for each initial guess c in Figure 1 (b), a corresponding number of steps is plotted.It stands for the number of iterations needed to compute an approximation of c ∗ (starting from thecorresponding c ) with a precision of 10 − in the L -norm.This experiment allow us to determine the computational effort necessary for the reconstructionof c ∗ with respect to distinct choices of c . The identification problem for the absorptioncoefficient is known to be mildly ill-posed [13, 21]. This fact is in agreement with the valuesplotted in Figure 1 (b), in the sense that the number of iterations necessary to achieve a goodquality reconstruction do not strongly oscillate with the initial guess.We conduct yet another experiment for identifying only the absorption coefficient. This time,we assume the exact solution of problem (1)–(2) to be given by the coefficient pair ( a ∗ , c ∗ ) inFigure 1 (d). The setup of the inverse problem remains the same (domain, available data, parameterto output operator, etc.).On the first run of the algorithm, see Figure 2 (a)–(c), the diffusion coefficient a ∗ is assumed to beexactly known. In this situation, the level set method is able to identify the absorption coefficient(see Figure 2 (a) for the evolution of the iteration error), and the iteration stagnates after that.The corresponding differences between the exact solution c ∗ and the initial guess c and betweenthe exact solution c ∗ and the final iterate c are plotted in pictures (b) and (c) respectively.On the second run, see Figure 2 (d)–(f), we use the approximation a ( x ) ≡ a ∗ and iterate to recover c ∗ . In this case, the level set method is still able to identify the absorptioncoefficient, however with a poorer accuracy. Once again, the iteration stagnates after the numericalconvergence is reached (see Figure 2 (d) for the evolution of the error). The corresponding differencesfor the initial guess c and for the final iterate c are plotted in pictures (e) and (f) respectively.Notice that the number of iteration steps needed to recover c ∗ (approximately 2000 in both runs)is much larger than in the previous experiment. This can be explained by the complexity of thegeometry of the support of c ∗ in this experiment [15]. This complexity and non smooth geometryalso influence the quality of the reconstruction. In what follows we consider the identification of the diffusion coefficient a , based on either totalor partial knowledge of c . In the first set of experiments, we consider problem (1)–(2) in the unit15quare with four pairs of NtD experiments, and the same exact solution ( a ∗ , c ∗ ) as in Section 5.1(see Figure 1 (a)).As initial guess for the level set method we choose distinct piecewise constant functions a ,whose supports are shown in Figure 1 (c). Analogous as in Section 5.1, the constant values of theexact solution a ∗ are assumed to be known, as well as the exact absorption coefficient c ∗ . Moreover,exact data are used for the reconstruction (i.e., δ = 0) and the scaling factors β j = 0.This time, the inverse problem reduces to a shape identification problem for the diffusioncoefficient. Once again we plot, for each initial guess a , a corresponding number of steps (seeFigure 1 (c)). It stands for the number of iterations needed to compute an approximation of a ∗ (starting from the corresponding a ) with a precision of 10 − in the L -norm.This experiment allow us to determine the computational effort necessary for the reconstructionof a ∗ with respect to distinct choices of a . The identification problem for the diffusioncoefficient is known to be exponentially ill-posed [13, 20]. This fact is in agreement with thevalues plotted in Figure 1 (c), meaning that the number of iterations necessary to achieve a goodquality reconstruction does strongly oscillate with the initial guess.We conduct yet another experiment for identifying only the diffusion coefficient. This time,we assume the exact solution of problem (1)–(2) to be given by the coefficient pair ( a ∗ , c ∗ ) inFigure 1 (d). The setup of the inverse problem remains the same (domain, available data, parameter L e rr o r || c * − c k || L X Y Difference c * − c k −− Iteration k = 0 Y Difference c * − c k −− Iteration k = 2500 (a) (b) (c) L e rr o r || c * − c k || L X Y Difference c * − c k −− Iteration k = 0 Y Difference c * − c k −− Iteration k = 2500 (d) (e) (f)Figure 2: Section 5.1, 2nd example. (a) – (c) Identification of c ∗ from exact knowledge of a ∗ . (a) Evolutionof the L error. (b) Difference between the initial guess c and c ∗ . (c) Difference between c and c ∗ . (d) – (f ) Identification of c ∗ from partial knowledge of a ∗ . (d) Evolution of the L error. (e) Differencebetween the initial guess c and c ∗ . (f ) Difference between c and c ∗ .
16o output operator, etc.).On the first run of the algorithm, see Figure 3 (a)–(c), the absorption coefficient c ∗ is assumed tobe exactly known. In this situation, the level set method is able to identify the diffusion coefficient(see Figure 3 (a) for the evolution of the iteration error), and the iteration stagnates after that.The corresponding differences between the exact solution a ∗ and the initial guess a and betweenthe exact solution a ∗ and the final iterate a are plotted in pictures (b) and (c) respectively.On the second run, see Figure 3 (d)–(f), we considered the approximation (¸ x ) ≡ c ∗ and iterate to recover a ∗ . In this case, the level set method is no longer able to identifythe diffusion coefficient. The iteration once again stagnates, but this time at some configurationfar from the exact solution (see Figure 3 (d) for the evolution of the L error). The correspondingdifferences for the initial guess a and for the final iterate a are plotted in pictures (e) and (f)respectively. In this last set of experiments we consider the level set algorithm for the simultaneous identificationof the coefficient pair ( a, c ) in (1)–(2). Three examples are considered and the corresponding exactsolutions are shown in Figure 4. The setup of the inverse problem is the same as in Sections 5.1 L e rr o r || a * − a k || L X Y Difference a * − a k −− Iteration k = 0 Y Difference a * − a k −− Iteration k = 5000 (a) (b) (c) L e rr o r || a * − a k || L X Y Difference a * − a k −− Iteration k = 0 Y Difference a * − a k −− Iteration k = 5000 (d) (e) (f)Figure 3: Section 5.2, 2nd example. (a) – (c) Identification of a ∗ from exact knowledge of c ∗ . (a) Evolutionof the L error. (b) Difference between the initial guess a and a ∗ . (c) Difference between a and a ∗ . (d) – (f ) Identification of a ∗ from partial knowledge of c ∗ . (d) Evolution of the L error. (e) Differencebetween the initial guess a and a ∗ . (f ) Difference between a and a ∗ . a ∗ , c ∗ ) is the one shown in Figure 4 (a). In order to devisean efficient iteration strategy for the simultaneous identification of both coefficients, we must takesome facts into account: F1)
From the 2nd example in Section 5.1, we have learned that the method for identifying c ∗ performs well, even if a good approximation for a ∗ is not known (see Figure 2 (d)–(f)). F2)
On the other hand, from the 2nd example in Section 5.2, we have learned that the level setmethod for identifying a ∗ performs well if a good approximation for c ∗ is available, but may generatea sequence a k that does not approximate a ∗ if a good approximation to c ∗ is not known. F3)
In the first run of the level set algorithm for the simultaneous identification of ( a ∗ , c ∗ ) weupdated both coefficients ( a k , c k ) in every step and observed that the iteration error k c k − c ∗ k decreases from the very first iteration. However, the iteration error k a k − a ∗ k only starts improvingwhen k c k − c ∗ k is sufficiently small.Thus, in order to save computational effort, we adopted the strategy to “freeze” the coefficient a k ( x ) = a ( x ) ≡ c k . We follow thisstrategy until the sequence c k stagnates (this is an indication that the iteration error k c k − c ∗ k issmall). In Figure 5 (a) and (e) this stage corresponds to the first k = 250 iterative steps (noticethat k a k − a ∗ k remains constant for k = 0 , . . . , k , while the difference c k − c ∗ is plotted in (g)).After this first iteration stage, we freeze c k = c k and iterate only with respect to a k . Thischaracterizes the second stage of the method. A natural question at this point would be: Why notto iterate with respect to both ( a k , c k ) for k ≥ k ? We tried to proceed in this way, but what weobserved is that: as long as k a k − a ∗ k does not significantly improve, the iterates c k stagnate with k c k − c k k almost constant.This second stage of the iteration can be observed in Figure 5 (a) and (e). Notice that k a k − a ∗ k decreases significantly, while k c k − c ∗ k remains constant for k = k , . . . , k = 750 (the difference a k − a ∗ is plotted in Figure 5 (c)).After the conclusion of the second iteration stage, the pair ( a k , c k ) is already a good approxi-mation for ( a ∗ , c ∗ ) (see Figure 5 (c) and (g)). As a matter of fact, this approximation is so good(a) (b) (c)Figure 4: Section 5.3. (a) – (c) Exact solution for the first, second and third example. a k , c k ), the iteration errors k a k − a ∗ k and k c k − c ∗ k are monotone decreasing. However, having in mind the convergence of thediffusion coefficient a k to the correct solution a ∗ takes more iterations than the absorption coeffi-cient c k , in this third stage we decided that each iteration step consists in one iteration with respectto the absorption coefficient c k and two iterations with respect to the diffusion coefficient a k (seeFigure 5 (a) and (e) for k ≥ k ).The introduction of this 3-stage iteration is motivated by the above mentioned facts (F1) –(F3). The calculation of optimal transition indexes k , k between the three stages is a difficulttask. However, since the degree of ill-posedness of the separate inverse problems for a and c is verydistinct from each other it’s not hard to get approximate values for k and k that will lead to alarge gain in computational effort by using this 3-stage strategy.The following second and third examples in this section do belong together. The correspondingexact solutions are shown in Figure 4 (b) and (c) respectively. Our goal is to investigate how thedistance between the supports of the exact solution pair ( a ∗ , c ∗ ) may interfere with the quality of thereconstruction of each single coefficient. In the second example there is a positive distance betweenthe supports, while in the third example both supports overlap.Although the distance between supp( a ∗ ) and supp( c ∗ ) in example 2 is smaller than in example 1above, the 3-stage iteration behaves similarly in both examples. The 1st-stage is ended after k =250 iterations, when the error k c k − c ∗ k has decreased considerably (see Figure 6 (e)). The 2nd-stagecorresponds to k ≤ k ≤ k = 750; at this point the difference between the exact solution a ∗ and theiteration a has visibly decreased (see Figure 6 (c)). The 3rd-stage of the iteration corresponds to k ≥ k . In this final stage, each iteration step consists in two iterations with respect to the diffusion L e rr o r || a * − a k || L X Y Difference a * − a k −− Iteration k = 0 Y Difference a * − a k −− Iteration k = 750 Y Difference a * − a k −− Iteration k = 2500 (a) (b) (c) (d) L e rr o r || c * − c k || L X Y Difference c * − c k −− Iteration k = 0 Y Difference c * − c k −− Iteration k = 250 Y Difference c * − c k −− Iteration k = 2500 (e) (f) (g) (h)Figure 5: Section 5.3, 1st example. (a) – (d) Iterative reconstruction of a ∗ . (a) Evolution of the L error. (b) Difference a − a ∗ . (c) Difference a − a ∗ . (d) Difference a − a ∗ . (e) – (h) Iterative reconstructionof c ∗ . (e) Evolution of the L error. (f ) Difference c − c ∗ . (g) Difference c − c ∗ . (h) Difference c − c ∗ . a k and one iteration with respect to the absorption coefficient c k . The final results canbe observed in Figure 6 (d) and Figure 6 (h) respectively.The third and last example reveled itself as the most difficult identification problem among allthree considered in this section. The solution pair ( a ∗ , c ∗ ) is chosen such that the supports of a ∗ and c ∗ intersect (see Figure 4 (c)). We start the iteration once again keeping a k constant during thefirst stage. This part of the method is successful, since after k = 500 iterations c k delivers a goodapproximation for the exact solution c ∗ (Figure 7 (g)). After that, we start iterating with respectto a k . After k = 750 iterations we observe that the error k a k − a ∗ k has decreased considerably(Figure 7 (a)). Finally, we start with 3rd-stage of the algorithm and this is the point where thedifficulties arise. No matter how many iterations we compute with respect to c k , the approximationdoes not get better than the one plotted in Figure 7 (g), which is computed after k = 250 steps.After 1500 steps, no significant improvement can be observed in the reconstruction of the absorptioncoefficient (compare Figure 7 (g) and (h)). In this last example, the reconstruction of the diffusioncoefficient is very precise, but the approximation obtained for the absorption coefficient is not sogood.It is worth noticing that the poor reconstruction of c ∗ is not due to non-stable behavior of our3-stage algorithm. The particular exact solution ( a ∗ , c ∗ ) in this example (with intersecting supports)leads to a very hard identification problem already reported in [31, 22, 3].It is worth mentioning that all problems presented in this Section were solved using the standardlevel set method described in Section 3.2, i.e., updating both ( a k , c k ) in every iterative step (andneglecting the 3-stage strategy). The final results of these iterations were basically the same as theones presented here. However, the computational effort involved in the computation was by far L e rr o r || a * − a k || L X Y Difference a * − a k −− Iteration k = 0 Y Difference a * − a k −− Iteration k = 750 Y Difference a * − a k −− Iteration k = 1500 (a) (b) (c) (d) L e rr o r || c * − c k || L X Y Difference c * − c k −− Iteration k = 0 Y Difference c * − c k −− Iteration k = 250 Y Difference c * − c k −− Iteration k = 1500 (e) (f) (g) (h)Figure 6: Section 5.3, 2nd example. (a) – (d) Iterative reconstruction of a ∗ . (a) Evolution of the L error. (b) Difference a − a ∗ . (c) Difference a − a ∗ . (d) Difference a − a ∗ . (e) – (h) Iterative reconstructionof c ∗ . (e) Evolution of the L error. (f ) Difference c − c ∗ . (g) Difference c − c ∗ . (h) Difference c − c ∗ . L e rr o r || a * − a k || L X Y Difference a * − a k −− Iteration k = 0 Y Difference a * − a k −− Iteration k = 750 Y Difference a * − a k −− Iteration k =1500 (a) (b) (c) (d) L e rr o r || c * − c k || L X Y Difference c * − c k −− Iteration k = 0 Y Difference c * − c k −− Iteration k = 250 Y Difference c * − c k −− Iteration k = 1500 (e) (f) (g) (h)Figure 7: Section 5.3, 3rd example. (a) – (d) Iterative reconstruction of a ∗ . (a) Evolution of the L error. (b) Difference a − a ∗ . (c) Difference a − a ∗ . (d) Difference a − a ∗ . (e) – (h) Iterative reconstructionof c ∗ . (a) Evolution of the L error. (f ) Difference c − c ∗ . (g) Difference c − c ∗ . (h) Difference c − c ∗ . In this paper, we develop a level set regularization approach for simultaneous reconstruction ofthe piecewise constant coefficients ( a, c ) from a finite set of boundary measurements of opticaltomography in the diffusive regime. From the theoretical point of view, we prove that the forwardmap F is continuous in the [ L (Ω)] topology. Hence, following standard arguments presented bythe authors in previous papers (see [8]) we get that the proposed level set strategy is a regularizationmethod. The main result behind the continuity of F is a generalization of Meyers’ Theorem for ourparticular case.On the other hand, we propose a numerical algorithm to reconstruct simultaneously the diffusionand absorption coefficients. Both coefficients are computed by minimizing a regularized energyfunctional. Motivated by the fact that the reconstruction of the absorption coefficient c is a mildlyill-posed inverse problem whereas the reconstruction of the diffusion coefficient a is exponentiallyill-posed, we present a split strategy that consists in freezing a = a and first iterate with respectto c until the iteration stagnate. Then, keep c = c k and start to iterate with respect to a untilstagnation of the iteration. Finally, iterate both coefficient. This numerical strategy has not onlydemonstrated that gives very good results but also reduces significantly the computational effort.The situation of non-convergence of the level set algorithm, that is when coefficients ( a, c ) havea crossing section (as in Subsection 5.3) is not an easy problem and it has already been reportedin [31, 22, 3]. We conjecture that the level set algorithm will improve its performance if enough21airs of Neumann-to-Dirichlet data are available. Since the situation with many measurements isnumerically demanding, a strategy like the one proposed in [25] could be more appropriated. Welet this problem for future and careful investigation. Acknowledgments
The work of J.P.A. was partially supported by grants from CONICET and SECY-UNC.A.D. acknowledges support from CNPq - Science Without Border grant 200815/2012-1, ARD-FAPERGS grant 0839 12-3 and CNPq grant 472154/2013-3.The work of M.M.A. was partially supported by CNPq grants no. 406250/2013-8, 237068/2013-3and 306317/2014-1.The authors would like to thanks Prof. Dr. Uri M. Ascher for the all discussions and valuablesuggestions.
A Proof of Theorem 1
The main purpose of this appendix is to show that under mild assumptions on the boundary(Neumann) data g the solution u of (1)–(2) belongs to W ,p (Ω) for some p > u ∈ H (Ω)).As far as we know, this type of regularity, namely u ∈ W ,p (Ω) for p >
2, goes back to thepioneering work of Meyers [23], for elliptic BVPs with Dirichlet boundary conditions. Later on,Gallouet and Monier [16] generalized Meyers’ result to Neumann BVPs. However, for the best ofthe authors knowledge there is no proof of such a result for the problem (1)–(2).The following proof was suggested by one of the anonymous referees. The authors are gratefulto him for this suggestion.
Proof of Theorem 1
Let u ∈ H (Ω) be the unique solution of (1)–(2). It clearly satisfies theweak formulation Z Ω a ∇ u · ∇ ϕdx + Z Ω cuϕdx = Z Γ gϕdσ ∀ ϕ ∈ H (Ω) . (27)Define now ˜ u := u − | Ω | R Ω udx , which in particular satisfies the weak formulation (5) in [16], i.e., ˜ u ∈ H ∗ (Ω) , Z Ω a ∇ ˜ u · ∇ ϕdx = h f, ϕ i ( H ) ′ ,H ∀ ϕ ∈ H (Ω) , where f is defined as h f, ϕ i ( H ) ′ ,H := Z Γ gϕdσ − Z Ω cuϕdx ∀ ϕ ∈ H (Ω) . Note that f naturally satisfies h f, i ( H ) ′ ,H = 0. Hence, to finish the proof we only need to applythe regularity result given in [16, Theorem 2]. To this end, it remains to show that the distribution22 is in W ,q (Ω) ′ . As g belongs to W − /q,q (Γ) ′ the distribution ϕ R Γ gϕdσ is in W ,q (Ω) ′ thanks tothe trace theorem in Sobolev spaces. Next we shall prove that the distribution h : ϕ R Ω cuϕ dx also belongs to W ,q (Ω) ′ . Consider first the case N = 2. In this case, since q <
2, we have thecontinuous embedding [4, Corollary 9.14] W ,q (Ω) ֒ → L q ∗ (Ω), where q ∗ = 2 q/ (2 − q ) >
2. Letting s := q ∗ / ( q ∗ −
1) be the conjugate of q ∗ we have s < q ∗ >
2) and, as a consequence,the continuous embedding [4] H (Ω) ֒ → L (Ω) ֒ → L s (Ω). Using the latter inclusions, the fact that u ∈ H (Ω), the second inequality in (4) and the H¨older’s inequality, we obtain (for all ϕ ∈ W ,q (Ω)): |h h, ϕ i| = (cid:12)(cid:12)(cid:12)(cid:12)Z Ω c ( x ) uϕ dx (cid:12)(cid:12)(cid:12)(cid:12) ≤ c k u k L s k ϕ k L q ∗ ≤ c k u k L s k ϕ k W ,q , (28)which proves that h ∈ W ,q (Ω) ′ . Consider now the case N ∈ { , } and let q ∗ := qN/ ( N − q ) > s := q ∗ / ( q ∗ −
1) its conjugate. In this case, we have also the continuous embeddings [4,Corollary 9.14] W ,q (Ω) ֒ → L q ∗ (Ω) and, since 1 ≤ s ≤ ∗ := 2 N/ ( N − H (Ω) ֒ → L s (Ω). Usingthe same reasoning as in the case N = 2 we find that (28) also holds when N ∈ { , } whichconcludes the proof of the desired regularity to the distribution h . Altogether, we obtain that f iswell-defined and belongs to W ,q (Ω) ′ . References [1] S. R. Arridge. Optical tomography in medical imaging.
Inverse Problems , 15(2):R41–R93,1999.[2] S. R. Arridge and W.R.B. Lionheart. Nonuniqueness in diffusion-based optical tomography.
Opt. Lett. , 23:882–4, 1998.[3] S. R. Arridge and M. Schweiger. A general framework for iterative reconstruction algorithmsin optical tomography, using a finite element method. In
Computational radiology and imaging(Minneapolis, MN, 1997) , volume 110 of
IMA Vol. Math. Appl. , pages 45–70. Springer, NewYork, 1999.[4] Haim Brezis.
Functional analysis, Sobolev spaces and partial differential equations . Universi-text. Springer, New York, 2011.[5] M. Burger and S. Osher. A survey on level set methods for inverse problems and optimaldesign.
European J. Appl. Math. , 16(2):263–301, 2005.[6] R. Dautray and J.-L. Lions.
Mathematical analysis and numerical methods for science andtechnology. Vol. 2 . Springer-Verlag, Berlin, 1988.[7] A. De Cezaro, A. Leit˜ao, and X.-C. Tai. On level-set type methods for recovering piecewiseconstant solutions of ill-posed problems. In X.-C. Tai, K. Mørken, K. Lysaker, and K.-A.Lie, editors,
Scale Space and Variational Methods in Computer Vision , volume 5667 of
LectureNotes in Comput. Sci. , pages 50–62. Springer, Berlin, 2009.238] A. De Cezaro, A. Leit˜ao, and X.-C. Tai. On multiple level-set regularization methods forinverse problems.
Inverse Problems , 25:035004, 2009.[9] A. De Cezaro, A. Leit˜ao, and X.-C. Tai. On piecewise constant level-set (pcls) methods for theidentication of discontinuous parameters in ill-posed problems.
Inverse Problems , 29:015003,2013.[10] A. De Cezaro and A. Leit˜ao. Level-set of L type for recovering shape and contrast in inverseproblems. Inverse Problems in Science and Enginnering , 20(4):517–587, 2012.[11] A. De Cezaro and A. Leit˜ao. Corrigendum: Level-set of L type for recovering shape andcontrast in inverse problems. Inverse Problems in Science and Enginnering , 21:1–2, 2013.[12] O. Dorn and D. Lesselier. Level set methods for inverse scattering—some recent developments.
Inverse Problems , 25(12):125001, 11, 2009.[13] H. W. Engl, M. Hanke, and A. Neubauer.
Regularization of inverse problems , volume 375 of
Mathematics and its Applications . Kluwer Academic Publishers Group, Dordrecht, 1996.[14] L.C. Evans and R.F. Gariepy.
Measure theory and fine properties of functions . Studies inAdvanced Mathematics. CRC Press, Boca Raton, FL, 1992.[15] F. Fr¨uhauf, O. Scherzer, and A. Leit˜ao. Analysis of regularization methods for the solutionof ill-posed problems involving discontinuous operators.
SIAM J. Numer. Anal. , 43:767–786,2005.[16] T. Gallouet and A. Monier. On the regularity of solutions to elliptic equations.
Rend. Mat.Appl. (7) , 19(4):471–488 (2000), 1999.[17] A. P. Gibson, J. C. Hebden, and S. R. Arridge. Recent advances in diffuse optical imaging.
Physics in Medicine and Biology , 50(4):R1–R43, 2005.[18] B. Harrach. On uniqueness in diffuse optical tomography.
Inverse Problems , 25(5):055010, 14,2009.[19] Jeremy C Hebden, Simon R Arridge, and David T Delpy. Optical imaging in medicine: I.experimental techniques.
Physics in Medicine and Biology , 42(5):825, 1997.[20] V. Isakov.
Inverse problems for partial differential equations , volume 127 of
Applied Mathe-matical Sciences . Springer, New York, second edition, 2006.[21] B. Kaltenbacher, A. Neubauer, and O. Scherzer.
Iterative regularization methods for nonlinearill-posed problems , volume 6 of
Radon Series on Computational and Applied Mathematics .Walter de Gruyter GmbH & Co. KG, Berlin, 2008.[22] V. Kolehmainen, S. R. Arridge, W. R. B. Lionheart, M. Vauhkonen, and J. P. Kaipio. Recoveryof region boundaries of piecewise constant coefficients of an elliptic PDE from boundary data.
Inverse Problems , 15(5):1375–1391, 1999. 2423] N. G. Meyers. An L p -estimate for the gradient of solutions of second order elliptic divergenceequations. Ann. Scuola Norm. Sup. Pisa (3) , 17:189–206, 1963.[24] Vasilis Ntziachristos, XuHui Ma, A. G. Yodh, and Britton Chance. Multichannel photoncounting instrument for spatially resolved near infrared spectroscopy.
Review of ScientificInstruments , 70(1):193–201, 1999.[25] F. Roosta-Khorasani, Kees van den Doel, and Uri M. Ascher. Stochastic algorithms for inverseproblems involving pdes and many measurements.
SIAM J. Scient. Comput. , 36:s3–s22, 2014.[26] F. Santosa. A level-set approach for inverse problems involving obstacles.
ESAIM ContrˆoleOptim. Calc. Var. , 1:17–33, 1995/96.[27] M. Schweiger and S. R. Arridge. Application of temporal filters to time resolved data in opticaltomography.
Phys Med Biol. , 44(7):1699–717, 1999.[28] M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy. The finite element model forthe propagation of light in scattering media: boundary and source conditions.
Med. Phys. ,22(11):1779–1792, 1995.[29] T. Tarvainen, B.T. Cox, J.P. Kaipo, and S.R. Arridge. Reconstructing absorption and scat-tering distributions in quantitative photoacoustic tomography.
Inverse Problems , 28:084009,2012.[30] K. van den Doel and U. M. Ascher. On level set regularization for highly ill-posed distributedparameter estimation problems.
J. Comput. Phys. , 216(2):707–723, 2006.[31] Yong Xu, Xuejun Gu, Taufiquar Khan, and Huabei Jiang. Absorption and scattering images ofheterogeneous scattering media can be simultaneously reconstructed by use of dc data.
Appl.Opt. , 41(25):5427–5437, 2002.[32] A. D. Zacharopoulos, S.R. Arridge, O. Dorn, V. Kolehmainen, and J. Sikora. Three-dimensionalreconstruction of shape and piecewise constant region values for optical tomography usingspherical harmonic parametrization and a boundary element method.