Linearized Reconstruction for Diffuse Optical Spectroscopic Imaging
LLinearized Reconstruction for Diffuse Optical SpectroscopicImaging
Habib Ammari ∗ Bangti Jin † Wenlong Zhang ‡ Abstract
In this paper, we present a novel reconstruction method for diffuse optical spectroscopic imagingwith a commonly used tissue model of optical absorption and scattering. It is based on lineariza-tion and group sparsity, which allows recovering the diffusion coefficient and absorption coefficientsimultaneously, provided that their spectral profiles are incoherent and a sufficient number of wave-lengths are judiciously taken for the measurements. We also discuss the reconstruction for imperfectlyknown boundary and show that with the multi-wavelength data, the method can reduce the influenceof modelling errors and still recover the absorption coefficient. Extensive numerical experiments arepresented to support our analysis.
Mathematics Subject Classification (MSC2000).
Keywords. diffuse optical spectroscopic imaging, reconstruction algorithm, group sparsity, imper-fectly known boundary.
Diffuse optical spectroscopy (DOS) is a noninvasive and quantitative medical imaging modality, for re-constructing absorption and scattering properties from optical measurements at multiple wavelengthsexcited by the near-infrared (NIR) light. NIR light is particularly attractive for oncological applicationsbecause of its deep tissue penetrance and high sensitivity to haemoglobin concentration and oxygenationstate [14]. DOS imaging effectively exploits the wavelength dependences of tissue optical properties (e.g.absorption, scattering, anisotropy, reduced scattering and refractive index), and such dependences havebeen measured and tabulated for various tissues [10, 21]. It was reported that the spectral dependenceof tissue scattering contains much useful information for functional imaging [13, 28, 29]. For example,dual-wavelength spectroscopy has been widely used to determine the absorption coefficient and hencethe concentrations of reduced hemoglobin and oxygenated hemoglobin in tissue [28]. Thus, DOS imag-ing holds significant potentials for many biomedical applications, e.g., breast oncology functional brainimaging, stroke monitoring, neonatal hymodynamics and imaging of breast tumours [12, 13, 14, 24].In many biomedical applications, it is realistic to assume that the absorption coefficient is linked tothe concentrations and the spectra of chromophores through a linear map (cf. (2.2) below), and thespectra of chromophores are known from experiments. Then the goal in DOS is to recover the individualconcentrations from the measurements taken at multiple wavelengths. This task represents one of themost fundamental problems arising in accurate functional and molecular imaging. Theoretically, verylittle is known about uniqueness and stability of DOS (see [5, 6, 7] for related work in electrical impedancetomography (EIT) and [1, 3, 4, 9] for quantitative photoacoustic imaging). ∗ Department of Mathematics, ETH Z¨urich, R¨amistrasse 101, CH-8092 Z¨urich, Switzerland( [email protected] ). † Department of Computer Science, University College London, London, WC1E 6BT, UK ( [email protected] ). ‡ Department of Mathematics, Southern University of Science and Technology, 1088 Shenzhen, Guangdong, P.R. China( [email protected] ). a r X i v : . [ m a t h . NA ] A ug ote that despite the linear dependence of absorption on the concentrations, the dependence of imag-ing data on the absorption coefficient is highly nonlinear. Further, it may suffer from severe ill-posednessand possible nonuniqueness; the latter is inherent to diffuse optical tomography of reconstructing simul-taneously the absorption and diffusion coefficients [8]. Thus, the imaging problem is numerically verychallenging. Nonetheless, there have been many important efforts in developing effective reconstructionalgorithms using multi-wavelength data to obtain images and estimates of spatially varying concentra-tions of chromophores inside an optically scattering medium such as biological tissues. These includestraightforward least-squares minimization [13], models-based minimization [23] and Bayesian approach[25]. Generally, there are two different ways to use the spectral measurements. One is to recover theoptical parameters at each different wavelength separately and then fit the spectral parameter model tothese optical parameters [9, 16, 26]; and the other is to express the optical parameters as a function ofthe spectral parameter model and then estimate the spectral parameter directly [23, 25, 31]. We referthe interested readers to the survey [15] and the references therein for detailed discussions.In this paper, we shall develop a simple and efficient linearized reconstruction method for DOS torecover the absorption and diffusion coefficients. We employ the diffusion approximation to the radiativetransfer equation for light transport, which has been used widely in biomedical optical imaging [8]. Ourmain contributions are as follows. First, we show that within the linearized regime, incoherent spectraldependence allows recovering the concentrations and diffusion coefficient simultaneously, provided that asufficient number of measurements are judiciously taken. However, generally there is no explicit criterionon the number of wavelength, except in a few special cases (see Remark 1). Second, we demonstratethat with multi-wavelength measurements, the chromophore concentrations can be still be reasonablyrecovered even if the domain boundary is only imperfectly known. Thus, DOS can partially alleviate thedeleterious effect of modeling errors, in a manner similar to multifrequency EIT [2]. Third and last, theseanalytical findings are verified by extensive numerical experiments, where the reconstruction is performedvia a group sparse type recovery technique developed in [2].The rest of the paper is organized as follows. In Section 2, we derive the linearized model, and discussthe conditions for simultaneous recovery (and also the one group sparse reconstruction technique). Thenin Section 3, we demonstrate the potential of the multi-wavelength data for handling modeling errors,especially imperfectly known boundaries. We show that the chromophore concentrations can still bereasonably recovered from multi-wavelength data, but the diffusion coefficient is lost due to the corruptionof domain deformation. Extensive numerical experiments are carried out in Section 4 to support thetheoretical analysis. Finally, some concluding remarks are provided in Section 5. In this section, we mathematically formulate the linearized multi-wavelength method in DOS.
First, we introduce the diffuse optical tomography (DOT) model. Let Ω ⊂ R d ( d = 2 ,
3) be an openbounded domain with a smooth boundary ∂ Ω. The photon diffusion equation for the photon fluence rate u in the frequency domain takes the following form: −∇ · ( D ( x, λ ) ∇ u ) + µ a ( x, λ ) u = 0 in Ω ,D ( x, λ ) ∂u∂ν + αu = S on ∂ Ω , (2.1)where ν is the unit outward normal vector to the boundary ∂ Ω, the nonnegative functions D and µ a denote the photon diffusion coefficient and absorption coefficient, respectively. In practice, the sourcefunction S ( x ) is often taken to be a smooth approximation of the Dirac function δ y ( x ) located at y ∈ ∂ Ω[27]. The parameter α in the boundary condition is formulated as α = − R R ) in DOT model, where R is a directionally varying refraction parameter [8]. Throughout, the parameter α is assumed to be2ndependent of the wavelength λ . The weak formulation of problem (2.1) is to find u ∈ H (Ω) := { v ∈ L (Ω) : ∇ v ∈ L (Ω) } such that (cid:90) Ω D ( x, λ ) ∇ u · ∇ v + µ a ( x, λ ) uv d x + (cid:90) ∂ Ω αuv d s = (cid:90) ∂ Ω Sv d s, ∀ v ∈ H (Ω) . In practice, the coefficients µ a and D are actually depending on the light wavelength λ . The opticalproperties of the tissue can be expressed using their spectral representations. Commonly used spectralmodels for optical properties can be written as [9, 23, 26, 31]: µ a ( x, λ ) = K (cid:88) k =1 µ k ( x ) s k ( λ ) , (2.2) µ s ( x, λ ) = µ s, ref ( λ/λ ref ) − b . (2.3)That is, the absorption coefficient µ a is expressed as a weighted sum of chromophore concentrations µ k and the corresponding absorption spectra s k of K known chromophore. The scattering coefficient µ s is given, according to Mie scattering theory, as being proportional to µ s, ref and (scattering) power − b of a relative wavelength λ/λ ref [11, 18, 29]. The coefficient µ s, ref is known as the reduced scatteringcoefficient at a reference wavelength λ ref , and it can be spatially dependent. In many types of tissues(e.g., muscle and skin tissues), the wavelength dependence of µ s has been measured and can often beaccurately approximated by µ s ( x, λ ) = a ( x ) λ − b , where the exponent b is recovered from experiments[10, 21].The optical diffusion coefficient D ( x, λ ) is given by D ( x, λ ) = [3( µ a + µ s )] − . The condition µ s (cid:29) µ a is usually considered valid in order to ensure the accuracy of the diffusion approximation to the radiativetransfer equation [15, 16, 8]. Hence, we assume below that the diffusion coefficient D ( x, λ ) has the form: D ( x, λ ) = d ( x ) s ( λ ) , (2.4)where the wavelength dependence s ( λ ) is known from experiments.In DOT experiments, the tissue under consideration is illuminated with M sources, and measurementsare taken at detectors. In this work, we assume for the sake of simplicity that the positions x n ofthe sources and detectors are the same, and are distributed over the boundary ∂ Ω. The spectroscopicinverse problem is to recover the spatially-dependent coefficient d ( x ) and the concentrations µ k ( x ) ofthe chromophores given the measured data u (corresponding to the known sources S ) on the detectorsdistributed on the boundary ∂ Ω measured at several wavelengths λ i . It is well-known that the inverseproblem of recovering both coefficients d ( x ) and µ k ( x ) is quite ill-posed [16], since two different pairsof scattering and absorption coefficients can lead to identical measured data. The multi-wavelengthmethod is a promising approach to resolve this challenging nonuniqueness issue. It is also reasonable toassume that the wavelength dependence s ( λ ) and the absorption spectra s k ( λ ) are linearly independentso as to distinguish the diffusion coefficient and the chromophore concentrations by effectively using theinformation contained in the multi-wavelength data. Now we derive the linearized DOS model, which plays a crucial role in the reconstruction technique. Wediscuss the cases of known and unknown diffusion coefficient separately.
First, we derive the linearized model with both diffusion coefficient D ( x, λ ) and concentrations µ k ( x ) beingunknown. For simplicity, we assume that the coefficient d ( x ) is a small perturbation of the background,which is taken to be 1, i.e., d ( x ) = 1 + δd ( x ) , δd ( x ) has a compact support in the domain Ω and is small (in suitable L p (Ω) norms).For the inversion, smooth approximations S n of the Dirac masses at { δ x n } Nn =1 are applied and thecorresponding fluence rates u n are measured on the detectors located at all x n over the boundary ∂ Ωto gain sufficient information about the diffusion coefficient D ( x, λ ) and absorption coefficient µ a ( x, λ ).That is, let { u n ≡ u n ( x, λ ) } Nn =1 ⊂ H (Ω) be the corresponding solutions to (2.1), i.e., (cid:90) Ω D ( x, λ ) ∇ u n · ∇ v + µ a ( x, λ ) u n v d x + (cid:90) ∂ Ω αu n v d s = (cid:90) ∂ Ω S n v d s, ∀ v ∈ H (Ω) . (2.5)Next we derive the linearized multi-wavelength model for the DOT problem based on an integralrepresentation. Let v m ≡ v m ( λ ) ∈ H (Ω) be the background solution corresponding to D ( x, λ ) ≡ s ( λ )and µ a ≡ S m , i.e., v m fulfils (cid:90) Ω s ( λ ) ∇ v m · ∇ v d x + (cid:90) ∂ Ω αv m v d s = (cid:90) ∂ Ω S m v d s, ∀ v ∈ H (Ω) . (2.6)Note that unless s ( λ ) is independent of the wavelength λ , the dependence of the background solution v m on the wavelength λ cannot be factorized out. Taking v = v m in (2.5) and v = u n in (2.6) andsubtracting the two identities yield s ( λ ) (cid:90) Ω δd ( x ) ∇ u n · ∇ v m d x + K (cid:88) k =1 s k ( λ ) (cid:90) Ω µ k ( x ) u n v m d x = (cid:90) ∂ Ω ( S n v m − S m u n )d s. Since δd and µ k s are assumed to be small, we can derive the approximations ∇ u n ( x, λ ) ≈ ∇ v n ( x, λ )and u n ( x, λ ) ≈ v n ( x, λ ) in the domain Ω (valid in the linear regime), and hence arrive at the followinglinearized model s ( λ ) (cid:90) Ω δd ( x ) ∇ v n · ∇ v m d x + K (cid:88) k =1 s k ( λ ) (cid:90) Ω µ k ( x ) v n v m d x = (cid:90) ∂ Ω ( S n v m − S m u n )d s. (2.7)Note that since (cid:82) ∂ Ω S m u n d s is the measured data on the detector located at x m and (cid:82) ∂ Ω S n v m d s canbe computed given the background spectra s ( λ ), the right-hand side of the model (2.7) is completelyknown and can be readily computed.The DOT imaging problem for the linearized model is to recover δd and the chromophore concen-trations { µ k } Kk =1 from { u n ( x, λ ) } Nn =1 on the boundary ∂ Ω at several wavelengths { λ q } Qq =1 . For thereconstruction, we divide the domain Ω into a shape regular quasi-uniform mesh of elements { Ω l } Ll =1 such that Ω = ∪ Ll =1 Ω l , and consider a piecewise constant approximation of the coefficient δd ( x ) and theconcentrations { µ k } Kk =1 of the chromophores as follows δd ( x ) ≈ L (cid:88) l =1 ( δd ) l χ Ω l ( x ) ,µ k ( x ) ≈ L (cid:88) l =1 ( µ k ) l χ Ω l ( x ) , k = 1 , . . . , K, where χ Ω l is the characteristic function of the l th element Ω l , and ( µ k ) l denotes the value of the concen-tration µ k of the k th chromophore in the l th element Ω l , so is ( δd ) l . Upon substituting the approximationinto (2.7), we have a finite-dimensional linear inverse problem s ( λ ) L (cid:88) l =1 ( δd ) l (cid:90) Ω l ∇ v n · ∇ v m d x + K (cid:88) k =1 s k ( λ ) L (cid:88) l =1 ( µ k ) l (cid:90) Ω l v n v m d x = (cid:90) ∂ Ω ( S n v m − S m u n )d s. M ( λ ), M ( λ ) and the data vector X . We use a singleindex j = 1 , . . . , J with J = N for the index pair ( m, n ) with j = N ( m −
1) + n , and introduce thesensitivity matrix M ( λ ) = [ M jl ] ∈ R J × L and M ( λ ) = [ M jl ] ∈ R J × L with its entries M jl given by M jl ( λ ) = (cid:90) Ω l ∇ v n · ∇ v m d x ( j ↔ ( m, n )) ,M jl ( λ ) = (cid:90) Ω l v n v m d x ( j ↔ ( m, n )) , which is independent of the wavelength λ . Likewise, we introduce a data vector X ( λ ) ∈ R J with its j thentry X j ( λ ) given by X j ( λ ) = (cid:90) ∂ Ω ( S n v m − S m u n )d s ( j ↔ ( m, n )) . By writing the vectors A = ( δd ) l ∈ R L and A k = ( µ k ) l ∈ R L , k = 0 , . . . , K , we obtain the followinglinear system (parameterized by the light wavelength λ ) M ( λ ) s ( λ ) A + M ( λ ) K (cid:88) k =0 s k ( λ ) A k = X ( λ ) . (2.8) If the diffusion coefficient D ( x, λ ) is known, then the goal is to recover the concentrations { µ k } Kk =1 ofthe chromophores. As before, we assume the unknowns µ k are small (in suitable L p (Ω) norms). We canrepeat the above procedure, except that the background solution v m ∈ H (Ω) is now defined by (cid:90) Ω D ( x, λ ) ∇ v m · ∇ v d x + (cid:90) ∂ Ω αv m v d s = (cid:90) ∂ Ω S m v d s, ∀ v ∈ H (Ω) . (2.9)Taking v = v m in (2.5) and v = u n in (2.9) and subtracting the two identities gives K (cid:88) k =1 s k ( λ ) (cid:90) Ω µ k ( x ) u n v m d x = (cid:90) ∂ Ω ( S n v m − S m u n )d s. Using the approximation u n ( x, λ ) ≈ v n ( x, λ ) in the domain Ω (which is valid in the linear regime), wearrive at the following linearized model K (cid:88) k =1 s k ( λ ) (cid:90) Ω µ k ( x ) v n v m d x = (cid:90) ∂ Ω ( S n v m − S m u n )d s. As before, we divide the domain Ω into a shape regular quasi-uniform mesh of elements { Ω l } Ll =1 suchthat Ω = ∪ Ll =1 Ω l , and consider piecewise constant approximations of the chromophore concentrations µ k : µ k ( x ) ≈ L (cid:88) l =1 ( µ k ) l χ Ω l ( x ) , k = 1 , . . . , K. Then we obtain the following finite-dimensional linear inverse problem K (cid:88) k =1 s k ( λ ) L (cid:88) l =1 ( µ k ) l (cid:90) Ω l v n v m d x = (cid:90) ∂ Ω ( S n v m − S m u n )d s. Using the sensitivity matrix M and the data vector X in (2.8), we get the following parameterizedlinear system M ( λ ) K (cid:88) k =1 s k ( λ ) A k = X ( λ ) . (2.10)5 .3 The linearized DOT with multi-wavelength data In the two linearized DOT inverse problems in Section 2.2, the vectors A (if d ( x ) is unknown) and { A k } Kk =1 are the quantities of interest and are to be estimated from the wavelength dependent data X ( λ ),given the spectra s ( λ ) and s k ( λ ). These quantities directly contain the information of the locationsand supports of δd ( x ) and all the chromophores µ k . Now we describe a procedure for recovering thecoefficient d ( x ) (if unknown) and the concentrations µ k ( x ) of the chromophores simultaneously. Thediffusion wavelength dependence s ( λ ) is always known (e.g., s ( λ ) = cλ b , where the parameter b isknown from experiments [10]).We first formulate the inversion method for the case an unknown diffusion coefficient D ( x ). Weconsider the case when all the absorption coefficient spectra { s k ( λ ) } Kk =1 are known. Suppose that wehave collated the measured data at Q distinct wavelengths { λ q } Qq =1 . We write S = ( s k ( λ q )) ∈ R K × Q , S = ( s ( λ )) ∈ R × Q . We also introduce the measured matrix X = [ X t ( λ ) . . . X t ( λ Q )] t ∈ R J × Q, and the vector of unknowns A = [ A t . . . A tK , A t ] t ∈ R L × ( K +1) , . Here, the superscript t denotes thematrix/vector transpose. Then we define a total sensitivity matrix M constructed by Q × ( K + 1) blocks.The ( i, j )th block of M is M ( λ i ) s j ( λ i ) for 1 ≤ i ≤ Q , 1 ≤ j ≤ K and the ( i, K + 1)th block of M is M ( λ i ) s ( λ i ) for 1 ≤ i ≤ Q . Hence, (2.8) yields the following linear system: M A = X. (2.11)Similarly, when the diffusion coefficient D ( x, λ ) is known, we define another total sensitivity matrix M constructed by Q × K blocks. The ( i, j )th block of M is M ( λ i ) s j ( λ i ) for 1 ≤ i ≤ Q , 1 ≤ j ≤ K andthe vector of unknowns A = [ A t . . . A tK ] t ∈ R L × K, . We get from (2.10) the following linear system: M A = X. (2.12) Remark 1.
Under the condition that the wavelength dependence of M ( λ ) and M ( λ ) can be factorizedout (e.g., s is independent of λ ) (and thus can be absorbed into the spectra s k ( λ ) ), the linear systemscan be decoupled to gain further insight. To see this, we consider the case (2.11) . Since all the spectraare assumed to be linearly independent, when a sufficient number of wavelengths { λ q } Qq =1 are judiciouslytaken in the experiment, the corresponding spectral matrix ˜ S t = [ S t S t ] is incoherent in the sense that Q ≥ K + 1 and rank( S ) = K + 1 and ˜ S is also well-conditioned. Then the matrix ˜ S has a right inverse ˜ S − . By letting ˜ Y = X ˜ S − , we obtain [ M A , M A ] = ˜ Y .
These are K + 1 decoupled linear systems. By letting ˜ Y = [ ˜ Y . . . ˜ Y K ] ∈ R J × ( K +1) , we have K + 1 independent (finite-dimensional) linear inverse problems M A = ˜ Y ,M A k = ˜ Y k , k = 1 , . . . , K, (2.13) where A represents the diffusion coefficient δd ( x ) and A k (for ≤ k ≤ K ) represents the k th chromophore µ k ( x ) . Note that each linear system determines one and only one unknown concentration A k . Similarly,for the case of a known diffusion coefficient, the matrix S has a right inverse S − , under the givenincoherence assumption. By letting Y = XS − , we obtain M A = Y. These are K decoupled linear systems. By letting Y = [ Y . . . Y K ] ∈ R J × ( K +1) , we have K independent(finite-dimensional) linear inverse problems M A k = Y k , k = 1 , . . . , K, (2.14) where A k (for ≤ k ≤ K ) represents the k th chromophore concentrations µ k ( x ) . Below, we describe a group sparse reconstruction method developed in [2] to solve the ill-conditionedlinear systems (2.11) and (2.12). 6 .4 Group sparse reconstruction algorithm
Upon linearization and decoupling steps (see Remark 1), one arrives at decoupled linear systems of theform: Dx = b, (2.15)where D = M or M ∈ R J × L is the sensitivity matrix, x = A k ∈ R L (0 ≤ k ≤ K ) is the unknownvector, and b = Y k ∈ R J (0 ≤ k ≤ K ) is a known measured data. These linear systems are oftenunder-determined, and severely ill-conditioned, due to the inherent ill-posed nature of the DOT inverseproblem. We adapt a numerical sparse method developed in [2] to solve (2.15).The algorithm takes the following two aspects into consideration:(1) Under the assumption that the unknowns δd and µ k are small, we may assume that x is sparse.This suggests to solve the minimization problemmin x ∈ Λ (cid:107) x (cid:107) subject to (cid:107) Dx − b (cid:107) ≤ (cid:15), where (cid:107) · (cid:107) denotes the (cid:96) norm of a vector. Here, Λ represents an admissible constraint on theunknowns x , since they are bounded from below and above, and (cid:15) > b .(2) In DOT applications, it is also reasonable to assume that each concentration of chromophore µ k isclustered, and this refers to the concept of group sparsity. The grouping effect is useful to removethe undesirable spikes typically observed for the (cid:96) penalty alone.Now we describe the algorithm, i.e., group iterative soft thresholding, listed in Algorithm 1, adaptedfrom iterative soft thresholding for (cid:96) optimization [17]. Here, N is the maximum number of iterations, w lk are nonnegative weights controlling the strength of interaction, and N l denotes the neighborhood ofthe l th element. We take w lk = β , for some β > β = 0 . N l consists of all elements inthe triangulation sharing one edge with the l th element. Since the solution x is expected to be sparse, anatural choice of the initial guess x is the zero vector. The regularization parameter γ plays a crucial rolein the performance of the reconstruction quality: the larger the value γ is, the sparser the reconstructedsolution becomes. There are several possible strategies to determine its value, e.g., discrepancy principleand balancing principle, or a trial-and-error manner [20].Below we briefly comment on the main steps of the algorithm and refer to [2] for details.Step 3 g j is a gradient descent update of x j , and s j > s j = 1 / (cid:107) D (cid:107) .Step 4 This step takes into account the neighboring influence.Step 5 ¯ d j indicates a grouping effect: the larger ¯ d jl is, the more likely the l th element belongs to the group.Step 6 This step rescales γ to be inversely proportional to ¯ d jl .Step 7 This step performs the projected thresholding with a spatially variable ¯ α j . P Λ denotes the pointwiseprojection onto the constraint set Λ and S λ for λ > S λ ( t ) = max( | t | − λ,
0) sign( t ).In our numerical experiments, we apply Algorithm 1 to the coupled linear systems (2.11) and (2.12)directly. This can be achieved by a simple change to Algorithm 1. Specifically, at Step 3 of the algorithm,instead we compute the gradient of the least-squares functional (cid:107) M A − Y (cid:107) by g j = A j − s j M t ( M x j − Y ) . The remaining steps of the algorithm are applied to each component A i independently. Note that onecan easily incorporate separately a regularization parameter γ on each component, which is useful sincethe diffusion and scattering coefficients are likely to have different magnitudes.7 lgorithm 1 Group iterative soft thresholding. Input D , b , W , N , γ , N and x ; for j = 1 , . . . , N do Compute the proxy g j by g j = x j − s j D t ( Dx j − b ); Compute the generalized proxy d j by d jl = | g jl | + (cid:88) k ∈N l w lk | g jk | ; Compute the normalized proxy ¯ d j by ¯ d j = max( d j ) − d j ; Adapt the regularization parameter ¯ α j by¯ α jl = γ/ ¯ d jl , l = 1 , . . . , L ; Update x j +1 by the group thresholding x j +1 = P Λ ( S s j ¯ α j ( g j )); Check the stopping criterion. end for Now, in order to show the potentials of DOS for handling modelling errors, we consider the case wherethe boundary of the domain of interest is not perfectly known. This is one type of the modelling errorsthat occurs whenever the positions of the point sources and detectors or the domain of interest are notperfectly modelled.We denote the true but unknown physical domain by (cid:101)
Ω, and the computational domain by Ω, whichapproximates (cid:101)
Ω. Next, we introduce a forward map F : (cid:101) Ω → Ω, (cid:101) x → x , which is assumed to be a smoothorientation-preserving map with a sufficiently smooth inverse map F − : Ω → (cid:101) Ω. We denote the Jacobianof the map F by J F , and the Jacobian of F with respect to the surface integral by J SF .Suppose now that the function (cid:101) u n ( (cid:101) x, λ ) satisfies problem (2.1) in the true domain (cid:101) Ω with the diffusioncoefficient (cid:101) D ( (cid:101) x, λ ), absorption coefficient (cid:101) µ a ( (cid:101) x, λ ) and source (cid:101) S ( (cid:101) x ), namely −∇ (cid:101) x · ( (cid:101) D ( (cid:101) x, λ ) ∇ (cid:101) x (cid:101) u n ( (cid:101) x, λ )) + (cid:101) µ a ( (cid:101) x, λ ) (cid:101) u n = 0 in (cid:101) Ω , (cid:101) D ( (cid:101) x, λ ) ∂ (cid:101) u n ∂ (cid:101) ν ( (cid:101) x ) + α (cid:101) u n ( (cid:101) x ) = (cid:101) S n ( (cid:101) x ) on ∂ (cid:101) Ω . (3.1)Here, (cid:101) S n is a smooth approximation of the Dirac mass at the true position (cid:101) x n . The wavelength-dependentabsorption coefficient (cid:101) µ ( (cid:101) x, λ ) also has a separable form related to the true concentrations (cid:101) µ k ( (cid:101) x ) of thechromophres: (cid:101) µ a ( (cid:101) x, λ ) = K (cid:88) k =1 (cid:101) µ k ( (cid:101) x ) s k ( λ ) , (3.2)where (cid:101) µ k are assumed to be small. Furthermore, the diffusion coefficient D takes the linear form: (cid:101) D ( (cid:101) x, λ ) = s ( λ )(1 + (cid:102) δd ( (cid:101) x )) . (3.3)8he weak formulation of problem (3.1) is given by: find (cid:101) u n ( · , λ ) ∈ H ( (cid:101) Ω) such that (cid:90) (cid:101) Ω (cid:101) D ( (cid:101) x, λ ) ∇ (cid:101) x (cid:101) u n · ∇ (cid:101) x (cid:101) v + (cid:101) µ a ( (cid:101) x, λ ) (cid:101) u n (cid:101) v d (cid:101) x + (cid:90) ∂ (cid:101) Ω α (cid:101) u n (cid:101) v d (cid:101) s = (cid:90) (cid:101) Ω (cid:101) S n (cid:101) v d (cid:101) s, ∀ (cid:101) v ∈ H ( (cid:101) Ω) . (3.4)In the experimental settings, (cid:101) u n is assumed to be measured on the boundary ∂ (cid:101) Ω. However, becauseof the incorrect knowledge of ∂ (cid:101) Ω, the measured quantity is in fact u n := (cid:101) u n ◦ F − restricted to thecomputational boundary ∂ Ω. Below we consider only the case that the domain Ω is a small variation ofthe true physical one (cid:101)
Ω, so that the linearized regime is valid. Specifically, the map F : (cid:101) Ω → Ω is givenby F ( (cid:101) x ) = (cid:101) x + (cid:15) (cid:101) φ ( (cid:101) x ), where (cid:15) is a small scalar and the smooth function (cid:101) φ ( (cid:101) x ) characterizes the domaindeformation. Further, let F − ( x ) = x + (cid:15)φ ( x ) be the inverse map, which is also smooth.In order to analyze the influence of the domain deformation on the linearized DOT problem, weintroduce the solution v m ∈ H (Ω) corresponding to ¯ D ( λ, x ) ≡ s ( λ ) and µ a ≡ S m being asmooth approximation of δ x m , i.e., v m fulfils (cid:90) Ω s ( λ ) ∇ v m · ∇ v d x + (cid:90) ∂ Ω αv m v d s = (cid:90) ∂ Ω S m v d s, v ∈ H (Ω) . (3.5)Now we can state the corresponding linearized DOT problem with an unknown boundary. The resultindicates that even for an isotropic diffusion coefficient (cid:101) D in the true domain (cid:101) Ω, in the computationaldomain Ω the equivalent diffusion coefficient D is generally anisotropic, and there is an additional per-turbation factor on the boundary ∂ Ω. Proposition 3.1.
Let µ a = (cid:101) µ a ◦ F − . The linearized inverse problem on the domain Ω is given by s ( λ ) (cid:18)(cid:90) Ω ( δd ( x ) + (cid:15) Ψ) ∇ v n · ∇ v m d x (cid:19) + (cid:90) ∂ Ω α(cid:15)ψv n v m d s + K (cid:88) k =1 s k ( λ ) (cid:90) Ω µ k ( x ) v n v m d x = (cid:90) ∂ Ω ( S n v m − S m u n )d s, (3.6) for some smooth functions Ψ : Ω → R d × d and ψ : ∂ Ω → R , which are independent of the wavelength λ .Proof. First, we derive the governing equation for the variable u n = (cid:101) u n ◦ F − in the domain Ω from(3.4). Let v = (cid:101) v ◦ F − ∈ H (Ω). By the chain rule, we have ∇ (cid:101) x (cid:101) u n ◦ F − = ( J tF ◦ F − ) ∇ x u n , where thesuperscript t denotes the matrix transpose. Thus, we deduce that (cid:90) (cid:101) Ω (cid:101) D ( (cid:101) x, λ ) ∇ (cid:101) x (cid:101) u n ( (cid:101) x ) ·∇ (cid:101) x (cid:101) v ( (cid:101) x )d (cid:101) x = (cid:90) Ω ( (cid:101) D ◦ F − )( x )( J tF ◦ F − )( x ) ∇ u n ( x ) · ( J tF ◦ F − )( x ) ∇ v ( x ) | det J F ( x ) | − d x = (cid:90) Ω ( J F ◦ F − )( x )( (cid:101) D ◦ F − )( x )( J tF ◦ F − )( x ) ∇ u n ( x ) · ∇ v ( x ) | det J F ( x ) | − d x = (cid:90) Ω D ( x, λ ) ∇ u n ( x ) · ∇ v ( x ) d x, where the transformed diffusion coefficient D ( x, λ ) is given by [30, 22, 2] D ( x, λ ) = (cid:32) J F ( · ) (cid:101) D ( · , λ ) J tF ( · ) | det J F ( · ) | ◦ F − (cid:33) ( x ) . Similarly, we obtain (cid:90) (cid:101) Ω (cid:101) µ a ( (cid:101) x, λ ) (cid:101) u n ( (cid:101) x ) (cid:101) v ( (cid:101) x )d (cid:101) x = (cid:90) Ω µ a ( x, λ ) u n ( x ) v ( x ) | det J F ( x ) | − d x, (cid:90) ∂ Ω α (cid:101) u n (cid:101) v d (cid:101) s = (cid:90) ∂ Ω αu n v | det J SF ( x ) | − d s. (cid:101) D ≡ δd is compactlysupported in the domain. From (3.4), it follows that u n satisfies (cid:90) Ω D ( x, λ ) ∇ u n · ∇ v + µ a ( x, λ ) u n v | det J F ( x ) | − d x + (cid:90) ∂ Ω αu n v | det J SF ( x ) | − d s = (cid:90) ∂ Ω S n v d s, ∀ v ∈ H (Ω) . (3.7)Then by choosing v = v m in (3.7) and v = u n in (3.5), we arrive at (cid:90) Ω ( D ( x, λ ) − s ( λ )) ∇ v n · ∇ v m d x + (cid:90) Ω µ a ( x ) v n v m | det J F ( x ) | − d x + (cid:90) ∂ Ω α ( | det J SF ( x ) | − − u n v d s = (cid:90) ∂ Ω ( S n v m − S m u n )d s. (3.8)Note that J F = I + (cid:15)J (cid:101) φ , and J F − = I + (cid:15)J φ = I − (cid:15)J (cid:101) φ ◦ F − + O ( (cid:15) ), since (cid:15) is small. It is known that | det J F | = 1 − (cid:15) div (cid:101) φ + O ( (cid:15) ) [19, equation (2.10)], we similarly derive | det J SF | = 1 + (cid:15)ψ + O ( (cid:15) ). Then D ( x, λ ) is given by D ( x, λ ) = (cid:101) D ( · , λ )(1 + (cid:15) div (cid:101) φ ( · )) − ( I + (cid:15) ( J (cid:101) φ ( · ) + J t (cid:101) φ ( · ))) ◦ F − ( x ) + O ( (cid:15) )= (cid:101) D ( · , λ )((1 − (cid:15) div (cid:101) φ ( · )) I + (cid:15) ( J (cid:101) φ ( · ) + J t (cid:101) φ ( · ))) ◦ F − ( x ) + O ( (cid:15) )= (cid:101) D ( · , λ )(1 + Ψ (cid:15) ) ◦ F − ( x ) + O ( (cid:15) ) , where Ψ = ( J (cid:101) φ + J t (cid:101) φ − div (cid:101) φI ) is smooth and independent of λ . This and the linear form of (cid:101) µ a ( (cid:101) x, λ ) in(3.3) yield D ( x, λ ) ≈ s ( λ ) I + (cid:15)s ( λ )Ψ( x ) and µ a ( x ) | det J F ( x ) | − = K (cid:88) k =1 µ k ( x ) s k ( λ ) + o ( (cid:15) ) , where we have used the assumption that the µ k are small. Upon substituting the above expressions into(3.8) and the approximations ∇ u n ≈ ∇ v n , u n ≈ v n in the domain, we obtain (3.6).By Proposition 3.1, in the presence of an imperfectly known boundary with its magnitude (cid:15) be-ing comparable with the concentrations { µ k } Kk =1 and the perturbation δd , the perturbed sensitivitysystem contains significantly modeling errors resulting from the domain deformation. Consequently, adirect inversion of the linearized model (3.6) is unsuitable. This issue can be resolved using the multi-wavelength approach as follows. Since (3.6) is completely analogous to (2.7), with the only differencelying in the additional terms in s ( λ ) (corresponding to the diffusion coefficient) and the edge perturba-tion (cid:82) ∂ Ω α(cid:15)ψv n v m d s . However, the edge perturbation on the boundary ∂ Ω can be treated as unknownscorresponding to an additional spectral profile s ∗ ( λ ) ≡
1. Thus, one may apply the multi-wavelengthapproach to recover the quantities of interest.Specifically, assume that the spectral profiles s , s , ..., s K , and s ∗ are incoherent. Then the method inSection 2.2 may be applied straightforwardly, since the right-hand side is known. However, the diffusionperturbation δd will never be properly reconstructed, due to the pollution of the error term (cid:15) Ψ (resultingfrom the domain perturbation). The concentrations of chromophores µ k corresponding to the wavelengthspectrum s k , k = 1 , . . . , K may be reconstructed, since they are affected by the deformation only throughthe transformation µ k = (cid:101) µ k ◦ F − . That is, the location and shape can be slightly deformed, provided thatthe deformation magnitude (cid:15) is small. Only the information of the diffusion coefficient is affected, andcannot be reconstructed. In summary, multi-wavelength DOT is very effective to eliminate the modellingerrors caused by the boundary uncertainty, at least in the linearized regime.We have discussed that the influence of an uncertain boundary in the case where both diffusion andabsorption coefficients are unknown. We can also analyze for the case with a known diffusion coefficient10igure 4.1: The true boundary shape, the positions of sources and detectors.similarly. Specifically, one may assume the deformed diffusion coefficient on the domain ¯ D ( x, λ ) = (cid:101) D ( (cid:101) x, λ ) ◦ F − and repeat the procedure of Proposition 3.1. We just give the conclusion: when the diffusioncoefficient is known, the domain deformation contributes to a perturbation inside the spectrum s , andthe boundary deformation pollutes the known diffusion term, and the concentration of chromophores µ k corresponding to the wavelength spectrum s k , k = 1 , . . . , K could be reconstructed. Now we show some numerical results to illustrate our analytical findings. The general setting for thenumerical experiments is as follows. The domain Ω is taken to be the unit circle Ω = { ( x , x ) : x + x < } . There are 16 point sources uniformly distributed along the boundary ∂ Ω; see Figure 4.1 for a schematicillustration of the domain Ω, the point sources and the detectors.Furthermore, we assume that the spectral profile s ( λ ) for the diffusion coefficient is s ( λ ) = 0 . λ b ,where the parameter b is known from experiments [10]. In all the examples below, we take b = 1 .
5. Wewill also see that µ s (cid:29) µ a is fulfilled in all the numerical examples. We take a directionally varyingrefraction parameter R = 0 .
2, so that α = − R R ) = 1 /
3. We use a piecewise linear finite elementmethod on a shape regular quasi-uniform triangulation of the domain Ω. The unknowns are representedon a coarser finite element mesh using a piecewise constant finite element basis. We measure the data u n ( x m , λ )(:= (cid:82) ∂ Ω S m u n d s ) on the detectors located at x m . The noisy data u δn ( x m , λ ) is generated byadding Gaussian noise to the exact data u † n ( x m , λ ) corresponding to the true diffusion coefficient D ( x, λ )and absorption coefficient µ a ( x, λ ) by u δn ( x m , λ ) = u † n ( x m , λ ) + η max l | u † n ( x m , λ ) − v n ( x m , λ ) | ξ n,m , where η is the noise level, and ξ n,m follows the standard normal distribution.We present the numerical results for the cases with known boundary and with imperfectly knownboundary separately, and for each case, we also present the examples with the diffusion coefficient beingknown and unknown. In Algorithm 1, we take a constant step size to solve (2.15). First, we show numerical results for the case with a perfectly known boundary shape. We shall testthe robustness of the algorithm against the noise, and show that the multi-wavelength approach couldreduce the deleterious effects of the noise in the measured data. The regularization parameter γ was11a) true µ k (b) recovered µ (c) recovered µ (d) recovered µ (e) recovered µ Figure 4.2: Numerical results for Example 4.1: (a) exact µ and µ of two chromophores; (b)–(c) recoveredresults with η = 1% noise level; (d)–(e) recovered results with η = 10% noise level.determined by a trial-and-error manner, and it was fixed at γ = 5 × − for the diffusion coefficient and γ = 1 × − for the absorption coefficient in all the numerical examples with perfectly known boundary.This algorithm is always initialized with a zero vector. Example 4.1.
Consider a known diffusion coefficient D ( λ, x ) = s ( λ ) = 0 . λ . , and two chromophoresinside the domain: the wavelength dependence of the chromophore on the top is s ( λ ) = 0 . λ , and theone on the bottom is s ( λ ) = 0 . λ − ; See Figure 4.2 for an illustration. We take measurements at Q = 3 wavelengths with λ = 1 , λ = 1 . and λ = 2 . The numerical results for Example 4.1 are presented in Figure 4.2. It is observed that the recoveryis very localized within a clean background even with 10% noise in the data, and the supports of therecovered concentrations of the chromophores agree closely with the true ones and the magnitudes are well-retrieved. Remarkably, the increase of the noise level from 1% to 10% does not influence much the shapeof the recovered concentrations. Therefore, if the given spectral profiles s k ( λ ) are sufficiently incoherent,the corresponding unknowns can be fairly recovered. This example also show that the proposed multi-wavelength approach is very robust to data noise, due to strong prior imposed by Algorithm 1.The next example shows the approach for reconstructing three chromophores inside the domain. Example 4.2.
Consider the case with a known diffusion coefficient D ( λ, x ) = s ( λ ) = 0 . λ . , and 3chromophores inside the domain:(i) The two chromophores on the top share the wavelength dependence s ( λ ) = 0 . λ , and the one onthe bottom has a second spectral profile s ( λ ) = 0 . λ − ;(ii) The wavelength dependence of the chromophore on the top right is s ( λ ) = 0 . λ − , of the topleft one is s ( λ ) = 0 . λ and of the bottom is s ( λ ) = 0 . λ − .We take measurements at Q = 3 wavelengths with λ = 1 , λ = 1 . and λ = 2 and the noise level is setto be η = 1% . The reconstruction results for Example 4.2 are shown in Figure 4.3. Figure 4.3 indicates that theunknowns corresponding to two or three spectral profiles can be fairly recovered in terms of both thesupports and magnitudes. In case (i), the two chromophores on the top share the wavelength dependence,12nd they are recovered simultaneously; whereas in case (ii), the chromophores have three incoherentwavelength dependences, and they can be recovered separately.The next example aims at recovering both diffusion and absorption coefficients, which is known to bevery challenging in the absence of multi-wavelength data.(a) true µ k (b) recovered µ (c) recovered µ (d) recovered µ (e) recovered µ (f) recovered µ Figure 4.3: Numerical results for Example 4.2: (a) exact µ and µ of two chromophores; (b)-(c) recoveredresults for case (i) (noise level η = 1%); (d)-(f) recovered results for case (ii) (noise level η = 1%). Example 4.3.
Consider the case of an unknown diffusion coefficient given by D ( λ, x ) = s ( λ )(1 +0 . δd ( x )) = 0 . λ . (1+0 . δd ( x )) . Similar to Example 4.1, consider two chromophores inside the domain:the wavelength dependence of the chromophore on the top is s ( λ ) = 0 . λ , and that of the bottom is s ( λ ) = 0 . λ − . The measurements are taken at Q = 3 wavelengths with λ = 1 , λ = 1 . and λ = 2 , and the noise level η is fixed at η = 1% . The numerical results for Example 4.3 are shown in Figure 4.4. Simultaneously reconstructing thediffusion and absorption coefficients is more sensitive to data noise, when compared with the case ofa known diffusion coefficient. Recall that the problem of recovering both coefficients is quite ill-posed[16]: two different pairs of scattering and absorption coefficients can give rise to identical measured data.However, it is observed from Example 4.3 that the multi-wavelength approach allows overcoming thisnonuniqueness issue, provided that the spectra are indeed incoherent.The next example shows that multi-wavelength data can mitigate the effects of the noise.
Example 4.4.
Consider the setting of Example 4.3, but with a noise level η = 30% . We study twodifferent numbers of wavelengths.(i) The measurements are taken at Q = 3 wavelengths with λ i = 1 + ( i − / , i = 1 , . . . , ;(ii) The measurements are taken at Q = 30 wavelengths with λ i = 1 + ( i − / , i = 1 , . . . , . Numerical results for Example 4.4 are shown in Figure 4.5. When using only data for 3 wavelengths,the recovered images are blurred by 30% noise. However, using data for 30 wavelengths, both the diffusioncoefficient δd and two chromophore concentrations µ k are much better resolved than using data with 3wavelengths. Hence, more wavelength observations can greatly mitigate the effects of data noise; whichconcurs with the observations from the experimental study [29].13a) true µ k (b) recovered µ (c) recovered µ (d) true δd ( x ) (e) recovered δd ( x )Figure 4.4: Numerical results for Example 4.3: (a) exact µ and µ of two chromophores; (b)-(c) recoveredresults for two chromophores; (d)-(e) true and recovered δd ( x ).(a) recovered δd ( x ) (b) recovered µ (c) recovered µ (d) recovered δd ( x ) (e) recovered µ (f) recovered µ Figure 4.5: Numerical results for Example 4.4: (a)-(c) the recovered results for case (i) (with 3 wave-lengths); (d)-(f): results for case (ii) (with 30 wavelengths).
Now we illustrate the approach in the case of an imperfectly known boundary. The (unknown) true domain (cid:101)
Ω is an ellipse centered at the origin with semi-axes a and b , E a,b = { ( x , x ) : x /a + x /b < } , andthe computational domain Ω is the unit disk. In this part, the regularization parameter γ was determinedby a trial-and-error manner, and it was fixed at γ = 5 × − for the diffusion coefficient, γ = 1 × − for the absorption coefficient and γ = 1 × − for the edge perturbation in all the numerical exampleswith imperfectly known boundary. This algorithm is always initialized with a zero vector. Example 4.5.
Consider the case of a known diffusion coefficient D ( λ, x ) = s ( λ ) = 0 . λ . , and two ifferent shape deformations: ( i ) (cid:101) Ω is an ellipse with a = 1 . and b = 0 . and ( ii ) (cid:101) Ω is an ellipse with a = 1 . and b = 0 . . Consider two chromophores inside (cid:101) Ω : the wavelength dependence of the chromophoreon the top is s ( λ ) = 0 . λ − , and that of the bottom is s ( λ ) = 0 . λ − . The measurements aretaken at Q = 3 wavelengths with λ = 1 , λ = 1 . and λ = 2 , and the noise level is fixed at η = 1% . The numerical results for Example 4.5 are shown in Figure 4.6. This example illustrates the influenceof the deformation scale on the reconstruction. The numerical results show clearly the potential of themulti-wavelength approach: Even using the wrong domain for the inversion step, we can still recover theconcentrations of the chromophores (or more precisely the deformed concentrations µ k = (cid:101) µ k ◦ F − ). Thenumerical results also show even we assume known diffusion coefficient in case (i), we should also use allthe spectra s , s k , and s ∗ to recover the right concentrations of the chromophores (Figures 4.6 (d) and(e)), or the results will be ruined by the shape deformation (Figures 4.6 (b) and (c)).(a) true µ k (b) recovered µ (c) recovered µ (d) recovered µ (e) recovered µ (f) true µ k (g) recovered µ (h) recovered µ Figure 4.6: Numerical results for Example 4.5: (a),(f) exact µ and µ of two chromophores in (cid:101) Ω; (b)-(c)recovered results for two chromophores in Ω for case (i) only using the spectra s k of the chromophores;(d)-(e) recovered results for two chromophores in Ω for case (i) using the spectra s k of the chromophores,the spectrum s of the diffusion coefficient and the spectrum s ∗ ( λ ) ≡ s , s k , and s ∗ . Example 4.6.
Consider the case of an unknown diffusion coefficient D ( λ, x ) = s ( λ )(1 + 0 . δd ( x )) =0 . λ . (1 + 0 . δd ( x )) , and the unknown true domain (cid:101) Ω is an ellipse with a = 1 . and b = 0 . . Considertwo chromophores inside the domain (cid:101) Ω : the wavelength dependence of the chromophore on the top is s ( λ ) = 0 . λ − , and that of the bottom is s ( λ ) = 0 . λ − . The measurements are taken at Q = 3 wavelengths with λ = 1 , λ = 1 . and λ = 2 , and the noise level is fixed at η = 1% . The numerical results for Example 4.6 are shown in Figure 4.7. It is observed that the two chro-mophores are recovered well in spite of the imperfectly known boundary, while the diffusion coefficient15s totally distorted by domain deformation, and thus it cannot be accurately recovered. The empiricalobservations on Examples 4.5 and 4.6 concur with the theoretical predictions in Section 3: Proposition3.1 implies that the domain deformation will be added to the recovered results corresponding to thespectral profile s ( λ ) and the diffusion coefficient cannot be recovered due to the domain deformation,but the deformed chromophore concentrations µ k = (cid:101) µ k ◦ F − can still be recovered.(a) true µ k (b) recovered µ (c) recovered µ (d) true δd (e) recovered δd Figure 4.7: Numerical results for Example 4.6: (a),(d): exact µ and µ of two chromophores in (cid:101) Ω;(b)-(c): recovered results for two chromophores in Ω; (e) the recovered result corresponding to s ( λ ). In this work, we have introduced a novel reconstruction technique for diffuse optical imaging with multi-wavelength data. The approach is based on a linearized model and a group sparsity approach. We haveshown that within the linear regime, our reconstruction technique allows recovering the concentrationof individual chromophore and the diffusion coefficients, provided that their spectral profiles are knownand incoherent. Furthermore, we have demonstrated that the multi-wavelength data can significantlyreduce modelling errors associated with an imperfectly known boundary. In fact, it allows recoveringwell (deformed) concentrations of the chromophores. These findings are fully supported by extensivenumerical experiments.
References [1] G. S. Alberti and H. Ammari. Disjoint sparsity for signal separation and applications to hybridinverse problems in medical imaging.
Appl. Comput. Harmon. Anal. , 42(2):319–349, 2017.[2] G. S. Alberti, H. Ammari, B. Jin, J.-K. Seo, and W. Zhang. The linearized inverse problem inmultifrequency electrical impedance tomography.
SIAM J. Imaging Sci. , 9(4):1525–1551, 2016.[3] H. Ammari, E. Bossy, V. Jugnon, and H. Kang. Mathematical modeling in photoacoustic imagingof small absorbers.
SIAM Rev. , 52(4):677–695, 2010.[4] H. Ammari, E. Bossy, V. Jugnon, and H. Kang. Reconstruction of the optical absorption coefficientof a small absorber from the absorbed energy density.
SIAM J. Appl. Math. , 71(3):676–693, 2011.165] H. Ammari, J. Garnier, L. Giovangigli, W. Jing, and J.-K. Seo. Spectroscopic imaging of a dilutecell suspension.
J. Math. Pures Appl. , 105(5):603–661, 2016.[6] H. Ammari, J. Garnier, H. Kang, L. Nguyen, and L. Seppecher.
Multi-Wave Medical Imaging ,volume 2 of
Modelling and Simulation in Medical Imaging . World Scientific, London, 2017.[7] H. Ammari and F. Triki. Identification of an inclusion in multifrequency electric impedance tomog-raphy.
Comm. Partial Differential Equations , 42(1):159–177, 2017.[8] S. R. Arridge. Optical tomography in medical imaging.
Inverse Problems , 15(2):R41–R93, 1999.[9] G. Bal and K. Ren. On multi-spectral quantitative photoacoustic tomography in diffusive regime.
Inverse Problems , 28(2):025010, 13, 2012.[10] A. N. Bashkatov, E. A. Genina, and V. V. Tuchin. Optical properties of skin, subcutaneous andmuscle tissues: a review.
J. Innovat. Opt. Health Sci. , 04(01):9–38, 2011.[11] F. Bevilacqua, A. J. Berger, A. E. Cerussi, D. Jakubowski, and B. J. Tromberg. Broadband absorp-tion spectroscopy in turbid media by combined frequency-domain and steady-state methods.
Appl.Optics , 39(34):6498–6507, 2000.[12] G. Boverman, E. L. Miller, A. Li, Q. Zhang, T. Chaves, D. H. Brooks, and D. A. Boas. Quantita-tive spectroscopic diffuse optical tomography of the breast guided by imperfect a priori structuralinformation.
Phys. Med. Biol. , 50(17):3941–3956, 2005.[13] A. E. Cerussi, D. B. Jakubowski, N. Shah, F. Bevilacqua, R. M. Lanning, A. J. Berger, D. Hsiang,J. A. Butler, R. F. Holcombe, and B. J. Tromberg. Spectroscopy enhances the information contentof optical mammography.
J. Biomed. Optics , 7(1):60–71, 2002.[14] A. E. Cerussi, V. W. Tanamai, D. Hsiang, J. Butler, R. S. Mehta, and B. J. Tromberg. Diffuseoptical spectroscopic imaging correlates with final pathological response in breast cancer neoadjuvantchemotherapy.
Phil. Trans. Royal Sci. , 369(1955):4512–4530, 2011.[15] B. Cox, J. Laufer, S. R. Arridge, and P. C. Beard. Quantitative spectroscopic photoacoustic imaging:a review.
J. Biomed. Opt. , 17(6):061202, 2012.[16] B. T. Cox, S. R. Arridge, and P. C. Beard. Estimating chromophore distributions from multiwave-length photoacoustic images.
J. Opt. Soc. Am. A , 26(2):443–455, 2009.[17] I. Daubechies, M. Defrise, and C. De Mol. An iterative thresholding algorithm for linear inverseproblems with a sparsity constraint.
Comm. Pure Appl. Math. , 57(11):1413–1457, 2004.[18] R. Graaff, J. G. Aarnoudse, J. R. Zijp, P. M. A. Sloot, F. F. M. de Mul, J. Greve, and M. H. Koelink.Reduced light-scattering properties for mixtures of spherical particles: a simple approximation de-rived from Mie calculations.
Appl. Optics , 31(10):1370–1376, 1992.[19] F. Hettlich. Fr´echet derivatives in inverse obstacle scattering.
Inverse Problems , 11(2):371–382, 1995.[20] K. Ito and B. Jin.
Inverse Problems: Tikhonov Theory and Algorithms . World Scientific, Hackensack,NJ, 2015.[21] S. L. Jacques. Optical properties of biological tisses: a review. 58(11):R37–R61, 2013.[22] V. Kolehmainen, M. Lassas, and P. Ola. The inverse conductivity problem with an imperfectlyknown boundary.
SIAM J. Appl. Math. , 66(2):365–383, 2005.[23] J. Laufer, B. Cox, E. Zhang, and P. Beard. Quantitative determination of chromophore concentra-tions from 2D photoacoustic images using a nonlinear model-based inversion scheme.
Appl. Optics ,49(8):1219–1233, 2010. 1724] T. O. McBride, B. W. Pogue, E. D. Gerety, S. B. Poplack, U. L. ¨Osterberg, and K. D. Paulsen.Spectroscopic diffuse optical tomography for the quantitative assessment of hemoglobin concentrationand oxygen saturation in breast tissue.
Appl. Optics , 38(25):5480–5490, 1999.[25] A. Pulkkinen, B. T. Cox, S. R. Arridge, J. P. Kaipio, and T. Tarvainen. A Bayesian approach tospectral quantitative photoacoustic tomography.
Inverse Problems , 30(6):065012, 18, 2014.[26] D. Razansky, M. Distel, C. Vinegoni, R. Ma, N. Perrimon, R. W. K¨oster, and V. Ntziachristos. Mul-tispectral opto-acoustic tomography of deep-seated fluorescent proteins in vivo.
Nature Photonics ,3(7):412–417, 2009.[27] J. Ripoll and V. Ntziachristos. Quantitative point source photoacoustic inversion formulas for scat-tering and absorbing media.
Phys. Rev. E , 71:031912, 9 pp., 2005.[28] E. M. Sevick, B. Change, J. Leigh, S. Nioka, and M. Maris. Quantitation of time-resolved andfrequency-resolved optical spectra for the determination of tissue oxygenation.
Anal. Biochem. ,195(2):330–351, 1991.[29] S. Srinivasan, B. W. Pogue, S. Jiang, H. Dehghani, and K. D. Paulsen. Spectrally constrained chro-mophore and scattering near-infrared tomography provides quantitative and robust reconstruction.
Appl. Optics , 44(10):1858–1869, 2005.[30] J. Sylvester. An anisotropic inverse boundary value problem.
Comm. Pure Appl. Math. , 43(2):201–232, 1990.[31] Z. Yuan and H. Jiang. Simultaneous recovery of tissue physiological and acoustic properties and thecriteria for wavelength selection in multispectral photoacoustic tomography.