IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control ,

vol. 61, no. 10,

October

2014

1627

Bayesian Approximation Error Approach in Full-Wave Ultrasound Tomography Janne Koponen, Tomi Huttunen, Tanja Tarvainen, and Jari P. Kaipio Abstract—In ultrasound tomography, the spatial distribution of the speed of sound in a region of interest is reconstructed based on transient measurements made around the object. The computation of the forward problem (the full-wave solution) within the object is a computationally intensive task and can often be prohibitive for practical purposes. The purpose of this paper is to investigate the feasibility of using approximate forward solvers and the partial recovery from the related errors by employing the Bayesian approximation error approach. In addition to discretization error, we also investigate whether the approach can be used to reduce the reconstruction errors that are due to the adoption of approximate absorbing boundary models. We carry out two numerical studies in which the objective is to reduce the computational times to around 3% of the time that would be required by a numerically accurate forward solver. The results show that the Bayesian approximation error approach improves the reconstructions.

I. Introduction

I

n ultrasound tomography, spatial distribution of the speed of sound in an object is estimated based on ultrasound measurements made around the object. The reconstruction of the speed of sound is a nonlinear inverse problem. There are several different techniques to solve the acoustic inverse problem, for example, the algebraic reconstruction technique (ART) [1]–[4], the Born approximation [5], [6], and the full-wave approaches [7]–[11]. The ART approximates wave distribution by a ray-based model. The Born approximation makes linear approximation assuming the scattered field to be much weaker than the incident field. Although the ART and the Born approximation are computationally lighter than the full-wave methods, these assumptions make them less accurate in some cases; for example, the Born approximation does not model multiple scattering. The rapid increase of computing power during recent years has made full-wave acoustic inversion possible even on an ordinary PC [10]. In reconstruction problems, and inverse problems in general, the accuracy of the forward solver is usually a central requirement and is often a critical topic. Accurate

Manuscript received January 28, 2014; accepted June 28, 2014. This work has been supported by the Academy of Finland (projects 136220, 140984, 272803, and 250215 Finnish Centre of Excellence in Inverse Problems Research) and by the strategic funding of the University of Eastern Finland. The authors are with the Department of Applied Physics, University of Eastern Finland, Kuopio, Finland (e-mail: [email protected]) T. Tarvainen is also with Department of Computer Science, University College London, London, UK. J. Kaipio is also with the Department of Mathematics, University of Auckland, Auckland, New Zealand. DOI http://dx.doi.org/10.1109/TUFFC.2014.006319 0885–3010

solvers are, however, often slow and their use can result in prohibitively long overall computational times. The nonlinear estimation problem can usually be formulated as a minimization problem with a high dimensional unknown. The adoption of the adjoint method [8], [11]–[14] can be used with acoustic inverse problems to speed up the computation of the update direction during the iteration but the related computational complexity inherits mainly that of the forward solver. Recently, the reduction of estimation errors that are due to employing (highly) approximate forward solvers has been investigated with a large number of different types of inverse problems. In this paper, we consider the reconstruction problem in the Bayesian inversion framework [15]–[17]. This framework makes it possible, in principle, to provide feasible (realistic) estimates for the related estimation errors. With feasibility, we require that the actual values are supported by the posterior model of the unknowns. The particular approach that we adopt here is the Bayesian approximation error (BAE) approach, in which we first model all the underlying low-level uncertainties [16], [18], [19]. Then, the approximate models are specified and the modeling errors are embedded in an additive error term. Subsequent approximate pre-marginalization over the approximation error is carried out and we are left with a possibly heavily approximative, but computationally efficient, model that ideally should be a feasible model. The Bayesian approximation error approach was originally introduced to cope with numerical model reduction in distributed parameter estimation problems [16], [19]. For example, in electrical impedance (resistance) tomography (EIT) and deconvolution problems, it was shown that significant model reduction is possible without essentially sacrificing the feasibility of estimates. The approach has also been applied to handle other kinds of approximation and modeling errors as well as other inverse problems. For example, model reduction, domain truncation, and unknown anisotropy structures and optode coupling coefficients in diffuse optical tomography were treated in [20]–[24] and modeling of errors in quantitavive photoacoustic tomography was considered in [25]. Missing boundary data in the case of image processing and geophysical ERT/EIT were considered in [26] and [27], respectively. Furthermore, in [28] and [29], the problem of recovery from simultaneous domain truncation and model reduction was found to be possible, and in [30] and [31], the recovery from the errors related to inaccurately known domain geometry was shown to be feasible. The approximation error approach was extended to nonstationary inverse problems in [32], in which linear

© 2014 IEEE

1628

IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control ,

vol. 61, no. 10,

October

2014

from sparse spatial discretization of the forward model and numerically too long time stepping. Furthermore, we also truncate the computational domain to very near the region of interest, and investigate how well one can recover from the errors that are related to employing approximate absorbing boundary conditions. The structure of the paper is as follows. In Section II, we define the wave propagation problem and the (noiseless) observation model that we consider in this paper. In Section III, we review the Bayesian approximation error approach. In Section IV, we consider the numerical methods that we use for the computation of the maximum a posteriori estimate for the distribution of the speed of sound and, in particular, the implementation of the solver for the adjoint problem. The numerical studies are given in Section V. Fig. 1. Schematic figure of measurement setup in UT. The transmitter s sends waveform f (t)δ(p − ps), which propagates through the object. The propagated signal is measured by sensors located on the dashed circle.

nonstationary (heat transfer) problems were considered, and in [33] and [34], in which nonlinear problems and state space identification problems were considered, respectively. A modification in which the approximation error statistics can be updated with accumulating information was proposed in [35] and an application to hydrogeophysical monitoring in [36]. In [37], aquifer parameters were estimated based on earthquake-related data, and it was found that the poroelastic model for an aquifer could be replaced with an elastic wave propagation model. Using this approach, it may also be possible to replace an accurate physical model with a less accurate but computationally more efficient model. For example, in [38] the radiative transfer model (Boltzmann transfer equation) was replaced with the diffusion approximation and in [39] corrections to linear methods in diffuse optical tomography were studied. Generally, in diffuse optical tomography and quantitative photoacoustic tomography, the physical model consists of two random fields which both must be estimated simultaneously. In [40] and [41], it was shown that one can pre-marginalize over one of the fields and estimate only the other field. In [18], the BAE was shown to be feasible also for the local X-ray tomography problem. The most comprehensive treatments of the Bayesian approximation error approach are given in [18] and [40]. In addition, there also exist other inversion methods that are based on the less-expensive surrogate models, such as the Neighborhood algorithm [42] or neural network inversion [43]; the first is a derivative-free method in same class as simulated annealing, and the latter uses a trained neural network to obtain posterior distribution statistics. Both of these methods have been successfully utilized in geophysical inverse problems. In this paper, we use transient wave data to estimate the speed of sound in an object. We consider the recovery

II. The Forward Problem Let Ms ultrasound transmitters and Mr receivers be located around the region of interest in Ω ⊂ R 2. A schematic figure of the measurement setup is given in Fig. 1. A transient waveform f (t) : [0, T ] → R is sent one at time by a transmitter s ∈ Js = {1, 2, …, Ms} at ps ∈ R 2, and the pulse is recorded by all receivers which are located at pr ∈ R 2. The receiver r ∈ Jr = {1, 2, …, Mr} measures pressure u(pr, t) over the time interval t ∈ [0, T ]. The pressure field u(p, t) : R 2 × R → R is modeled by acoustic wave equation

−∇ 2u(p, t) +

1 ∂2 u(p, t) = f (t)δ(p − p s ), c 2(p) ∂t 2

where p ∈ R 2, with initial conditions

 ∂ u(p, 0) = 0  ∂t   u(p, 0) = 0,

and where δ is the Dirac delta distribution and c is the speed of sound. The computational domain is truncated by an artificial exterior boundary and absorbing boundary conditions (ABC) are used on this boundary. The ABC simulates the unbounded domain by attempts to absorb the outgoing waves. There are several models for ABCs, such as the perfectly matched layer [44], [45], one-way wave equations [46], [47], and absorbing layers near boundaries [48]. In this paper, we employ the absorbing layers near boundaries and a first-order absorbing boundary condition which satisfies

1 ∂   ∂  − u  ∂n c ∂t 

p ∈∂Ω

= 0,

where n is the outward normal on the boundary.

koponen et al.: bayesian approximation error approach in full-wave ultrasound tomography

Let us( p, t) be the pressure field that is generated when the waveform is sent from transmitter s. For each s, we have

−∇ 2u s(p, t) +

1 ∂2 u s(p, t) = f (t)δ(p − p s ). (1) c (p) ∂t 2 2

The measured data are collected at times tk = (k − 1)Δt, where k = 1, …, K and Δt = T/(K − 1). We do not model the hydrophones explicitly and assume that the noiseless measurements (model predictions) can be written in the form R M r  ys(t k ) = Mu s(p, t k ), where the matrix M samples or interpolates the wave field at points pr. We collect all model predictions in a single vector y ∈ R M sM rK . Following [10], we formulate the estimation problem in terms of the (squared) slowness x(p) = c−2(p), which gives the simplest form to compute steepest descent with the adjoint method.

III. Bayesian Approximation Error Approach In the Bayesian framework for inverse problems, all unknowns are modeled explicitly as random variables, and operators that are not known exactly as stochastic operators. The modeling tasks include the construction of the likelihood model (conditional distribution of the measurements given the unknown) π (y | x), which includes the forward model and a model for all errors, and the prior model π (x), which is a probabilistic model for the unknowns. All inference is based on the posterior model π (x | y) ∝ π (y | x)π (x); see, for example, [15]–[17]. In this paper, we will use computationally (numerically) highly approximative models for the forward map. We denote accurate representations of variables, mappings, and operators with a bar, whereas the representations that are used in the inversion do not have the bar. Therefore, we write the measurement model in the form y = Aλ(x ) + e

= Aλ0(x ) + (Aλ(x ) − Aλ0(x )) + e (2) = Aλ0(x ) + ε + e,

where we define the approximation error ε = ψ(x , λ) = Aλ(x ) − Aλ0(x ), e denotes measurement noise, and λ denotes other (uninteresting) unknowns that we do not intend to estimate simultaneously with the primary unknowns x. For more details, see [18] and [40]. Thus, the approximation error is the difference of predictions of the measurements when using the accurate model Aλ(x ) and the approximate model Aλ0(x ). The accurate model, in this paper, refers to a numerical model, preferably of as high quality as is possible, which can be utilized in computing sample-based estimates of modeling error; see Section IVA. More generally, the preferred accurate model is the

1629

best available physical and numerical model which is computationally feasible for the problem at hand. In this paper, Aλ(x ) is the numerically accurate wave propagation model and λ refers to accurate absorbing boundary conditions. In practice, we can also have such a large computational domain that the errors that are due to using approximate absorbing boundary models are smaller than the measurement errors, and they can therefore be neglected. The model Aλ0(x ) is a numerically inaccurate model in a (possibly too) small computational domain and we can expect significant errors from the approximate boundary conditions λ0. When computing the estimates, we approximate the accurate representation of the primary unknown x by x = Px, where P is typically a projection operator. When the models A and A are fixed, we have formally π(ε | x ) = δ (ε − ψ(x )) where δ (·) is the Dirac distribution and the random variable λ is embedded in ε. We assume that the measurement errors are mutually independent with x. Then, the Bayes’ formula gives

π(y | x ) =

∫ πe(y − Aλ (x) − ε)π ε|x(ε |x)dε; 0

for derivation, see [18] and [40]. In the Bayesian approximation error approach, we approximate both πe and πε|x with normal distributions

e ∼ N (e ∗, Γe), ε | x ∼ N (ε∗,x, Γ ε|x ),

where the notation z ~ N (z ∗, Γ z ) means that the random variable (vector) z is normally distributed with mean z ∗ and covariance (matrix) Γz, and

1 ε∗,x = ε∗ + Γ εx Γ − x (x − x ∗) 1 Γ ε|x = Γ ε − Γ εx Γ − x Γ x ε.

Using normally distributed modeling error is a highly approximate model. The same approximation has been, however, utilized successfully in numerous studies applying the BAE, for example, [19], [22], [26], and [27]. We define the normal random variable ν with ν | x = e + ε | x as follows: ν | x ∼ N (ν ∗|x , Γ ν |x )

1 ν ∗|x = e ∗ + ε∗ + Γ εx Γ − x (x − x ∗) 1 Γ ν |x = Γe + Γ ε − Γ εx Γ − x Γ xε .

Thus, we obtain the approximate likelihood model

y | x ∼ N (y − Aλ0(x ) − ν ∗|x, Γ ν |x ).

Because our objective is a computationally efficient problem, we would usually employ also a normal prior model x ~ N (x ∗, Γ x ). Thus, we obtain the approximation for the posterior distribution (model):

1630



IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control ,

π(x |y) ∝ π(y | x )π(x )  1 2 ∝ exp  − ( L ν |x(y − Aλ0(x ) − ν ∗|x ) (3)  2 + Lx(x − x ∗)

2

)

),

T −1 T 1 where Γ − ν|x = L ν |x L ν |x , Γ x = L x L x , and ν ∗|x = ν ∗|x (x ). Often, the further approximations ν ∗|x ≈ e ∗ + ε∗ and Γν | x ≈ Γν are found to be adequate, and we also adopt them in this paper. In this paper, we assume that λ is related to the approximate absorbing boundary conditions and discretization of computational domain only, and that there are no other secondary unknowns. The statistics of the approximation errors are usually approximated by Monte Carlo simulation and by computing the associated sample mean and covariance. In the case of large dimensional measurements, which is the case also in this paper, the computation and storing of the matrix square root Lν is not a straightforward task. This topic will be discussed in Section IV. If the approximation errors are neglected, the posterior model is of the form

 1 2 2  π(x |y) ∝ exp  − ( σe−2 (y − Aλ0(x )) + Lx(x − x ∗) )  ,  2  (4) assuming zero mean and mutually uncorrelated measurement errors with equal variances σe2. In this paper, the reconstructions obtained as maximizers of (3) and (4), or, equivalently, the minimizers of the related posterior potentials, are compared against each other. IV. Computation of the Maximum a Posteriori Estimates A. Computing and Handling the Approximation Error Covariance In the approximation error approach, the prior model π(x ) that is used to compute the approximation errors is generally not the same as the one used in the inversion— that is, the model π (x) in the posterior model (3). For example, the model π(x ) can be based on an anatomical database or ad hoc simulated distribution, which is also the case in this paper. The model π(x ) is described in Section V. To compute the approximate second-order statistics of ε, we draw C samples x () from π(x ) and compute ε(ℓ) = Aλ(x ()) − Aλ0(x ()), ℓ = 1, …, C, the sample mean ε∗, and the sample covariance Γε. In principle, the matrix square root Lν contains KMsMr × KMsMr real values. By assuming realistic values for K, Ms, and Mr, we can see that Lν cannot be stored in

vol. 61, no. 10,

October

2014

computer memory without packing it. However, Γν = Γe + Γε, where Γe is typically diagonal, and rank (Γε) = C ≪ KMsMr. This gives a method to decompose the matrix in such a way that it can be stored using less memory. For the treatment of large dimensional measurements in the context of BAE, see [18]. Let

Vˆε =

1 (ε (1) − ε*, …, ε (C ) − ε*) = VεSεU Tε , C −1

where VεSεU T is singular value decomposition of Vˆε . We use the truncated SVD-based method described in [16], choose D ≤ C, and let S = diag (s1, s2, …, sD), where si is the ith singular value of Vˆε . In addition, let Vε = ⋅ ⋅ (Vε(,1) Vε(,2) … Vε(,⋅ D)), where Vε(,⋅ i) is the ith column of Vε. We get a truncated eigenvalue decomposition



Γ ε = V εΛ εV εT ,

 V ε1    V ε =    ∈ R KM rM s×D,  M  V ε s 

and Λε = S ε2 ∈ R D×D is a diagonal matrix. The use of the Sherman–Morrison–Woodbury identity [49] gives 1 −1 Γ− ν = (Γ e + Γ ε)



= (Γe + V εΛ εV εT )−1 T −1 1 −1 T −1 = Γe−1 − Γe−1V ε(Λ − ε + V ε Γ e V ε) V ε Γ e

(5)

= Γe−1 − Γe−1V εΣ εV εT Γe−1, 1 T −1 −1 where Σε = (Λ − ∈ R D×D. The matrices Γ −1 ε + V ε Γ e V ε) ν and Lν can be divided into Ms × Ms block structures, with each block corresponding to data measured while one transmitter sent the pulse. The definition of Lν and block matrix multiplication gives for the matrix block 1 (l,m) (Γ − = (LTν L ν )(l,m) ν )



Ms

=



∑(L(νn,l))T L(νn,m).

n =1

On the other hand, using (5) gives the matrix block 1 (l,m) (Γ − = (Γe−1 − Γe−1V εΣ εV εT Γe−1)(l,m) ν )



= (Γe−1)(l,m) − (Γe−1V εΣ εV εT Γe−1)(l,m) = (Γe−1)(l,m) − V νlΣ ε(V νm )T ,

where

 V ν1    V ν =    = Γe−1V ε ∈ R KM rM s×D.  M s  V ν 

Combination of these gives

koponen et al.: bayesian approximation error approach in full-wave ultrasound tomography Ms

∑(L(νn,l))T L(νn,m) = (Γe−1)(l,m) − VνlΣ ε(Vνm)T . (6)



m =1

with final conditions  Ψ n(q, K ) = 0, ∀q ∈ J q   ∂ Ψ (q, K ) = 0, ∀q ∈ J , n q  ∂t



If C is large, most standard SVD solvers use too much memory for a typical computer. However, the incremental SVD algorithm proposed in [50] can be used in such cases. B. Implementation of the Adjoint Method The MAP estimates xMAP are the minimizers of the posterior potentials in (3) and (4). Gauss–Newton and Newton–Raphson approaches require the Jacobian or the Hessian matrix to be computed for the determination of the search direction. However, in our case the dimension of parameters x is too large to compute these matrices, and thus we use an adjoint method to compute the steepest descent search direction. In the following, we discuss how the approach described in [12] can be adapted for the computation of the MAP estimates. We give the formulas only for the more complicated BAE version. The posterior potential X is of the form Ms

X =



∑X n ,

where Jq = {(i, j) : i, j ∈ {1, 2, …, N}}, with absorbing boundary conditions. Then, the update direction can be written in the form

Ms K   ∂Ψ n ∂u  −1 ∆x = −  Γ x (x − x *) + ∆t (t k ) n (t k )  . ∂t ∂t   n =1 k =1

∑∑

Here, zn are the weighted residual source fields

 Ms Ms   z n = −R  (L(νn,l ))T L(νn,m)(g m − Am(x ))  ,  m =1 l =1 

∑∑

where the operator R : R KM r → (Jq × Jk → R) converts the weighted residual signals into the source fields, and is defined by

(Rw s )(q, k) =

n =0

X0 =

1 L (x − x *) 2 x

1 Xn = 2

2

Ms

∑

L(νn,m)(g m m =1



− Am(x )) ,

xˆMAP = argmin(X ).



x

The Fréchet derivative is of the form 1 D xX(x ′) = 〈x ′, Γ − x (x − x *)〉 Ms

+

∑

K

x ′, ∆t

n =1

=

∑

k =1

1 〈x ′, Γ − x (x

− x *)〉

+ x ′, ∆t

∑∑

Ms K

n =1 k =1

∂Ψ n ∂u (t ) n (t ) ∂t k ∂t k



∂Ψ n ∂u (t ) n (t ) , ∂t k ∂t k N2

where 〈·,·〉 is standard inner product in R . Direction of the steepest descent is obtained by choosing x′ so that the value of the Fréchet derivative is minimized. Let Ψn­(q, k) be the solution of the adjoint equation

s,r,k ,

when r ∈ J r so that q r = q otherwise.

C. Algorithm

2

where n = 1, …, Ms and gm = ym − ν *m . Thus, we have



{w0,

The residual signals ws can be computed using (6).

where



1631

−∇ 2Ψ n + x

∂2 Ψ n = z n , (7) ∂t 2

Now, the algorithm to solve the inverse problem in BAE framework can be written. The algorithm is similar to the algorithm presented in [10], however, some details are different. The procedure to solve the inverse problem is as follows: 0) Construct the prior model π(x ), draw C samples, and compute the sample mean ν *, and covariance Γν. Choose a normal prior model π(x) = N (x *, Γ x ). This step only needs to be done once for each measurement geometry. 1) Measure y. 2) Set x0, and i = 0. 3) Compute the residual ri = y − A(xi) − ν *. 4) Compute the update direction Δx with the adjoint method. a) Solve Ψn from the adjoint equation (7). b) Set

∆x(q) = ∆t

Ms K

∑∑

n =1 k =1

∂Ψ n ∂u 1 (q, k) n (q, k) + Γ − x (z − x *). ∂t ∂t

5) Set x i +1 = x i − κ∆x , where κ ∈ R is found using a line search algorithm. 6) Increase i by 1. 7) If not converged, go to step 3.

1632

IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control ,

As the initial step (step 0) of the algorithm, the statistical parameters are defined for later usage. Note that the step is done only once; the computations are independent of the measurement data, and the results can be stored for later use. They are valid for the employed geometry and discretization as long as the employed prior model can be assumed to be feasible. In addition, because the step is done only once, the computational load of the step is not as critical as in other steps. In step 2, an initial guess is set. In our simulations, we used an initial guess of x0 = 1/c 02, where c0 is the speed of sound in water. It has also been proposed that initial guess can be adopted by some other method, for example, by ray-based reconstruction [10]. Computing the residual, in step 3, requires forward problem (1) to be solved with parameters xi. We utilize pressure field un later in step 4. Although, it is possible to save the pressure field u on every time step, it requires a lot of memory. However, it is possible to compute un in reversed time by saving un on the boundary of the computational domain, and saving both un(T) and (∂/∂t)un(T) on the whole computational domain [10], [51]. In step 4, the direction of the steepest descent is obtained by solving the adjoint equation. It is also possible to utilize other type of optimization schemes, for example, the conjugate gradient method. In step 5, the length of the step is defined. We use a golden section search [52] to define the optimal step length. In the current case, the speed of the line search algorithm is more important than the accuracy. We use two methods to speed up the line search; number of steps in the line search is limited, and we use an approximate cost function to minimize. The approximate cost function

ˆ x) = X(

∑ˆ ∑ L(νn,m)(Am(x) − g m),

n ∈J s m ∈J r

where Jˆs ⊂ Js contains only a few indices. In this work, the general purpose GPU (GPGPU) [53], [54] is utilized to speed up solving the wave equation and the adjoint problem in steps 3 and 4. There exist several articles about solving the full-wave inversion utilizing the GPGPU, for example, [10] and [55].

vol. 61, no. 10,

October

2014

V. Numerical Examples A. Simulation of Measurement Data We considered a numerical example which simulates detection of a breast tumor. The target is shown in Fig. 2. Parameters were chosen to imitate a human breast and the speed of sound in each sub-domain corresponds to an anatomically realistic value [56]. The outermost layer modeled skin with speed of sound 1580 m/s, and the middle layer modeled fat tissue with speed of sound 1480 m/s. For the innermost layer, the breast tissue, speed of sound was 1525 m/s. Speed of sound inside the star-like and circular inclusions simulating tumors varied from 1460 m/s to 1570 m/s. The breast was surrounded by water with speed of sound 1500 m/s. A rectangular computation domain Ω1 = 0.33 × 0.33 m was used, and Mr = 64 receivers and Me = 16 transmitters were located evenly on a circle with radius 0.1 m. Finite difference method in time domain (FDTD) with interpolated isotropic scheme (IIS) [57] with spatial discretization of 1858 × 1858 points and 4096 time steps was used to compute the signals. Absorbing layers near the boundaries of the computational domain were employed. As the measurement error model, we used the additive Gaussian noise e ~ N (0, σe2I ), where σe = 1, which corresponds to an error level of 2.6% of the maximum signal. B. Computation of Approximation Error Statistics Computing approximation error statistics in each test case utilized two different forward models; an accurate forward model Aλ(x ), and an approximate forward model Aλ0(x ). The approximate model varied between the test cases, but the same accurate forward model Aλ(x ) was utilized in both cases. The accurate model utilized FDTD with IIS with spatial discretization of 1794 × 1794 and 4096 time steps with computational domain Ω1. Absorbing layers near the boundaries of the computational domain were employed. Receivers and transmitters were located as they were when computing the simulated measurement signal.

Fig. 2. Speed of sound reconstructions of simulations in case 1. Images from left to right are the numerical model, reconstruction with the approximate forward solver utilizing CEM, reconstruction with the approximate forward solver utilizing BAE, and reconstruction with the accurate forward solver.

koponen et al.: bayesian approximation error approach in full-wave ultrasound tomography

1633

Fig. 3. Four draws form the prior model π(x ). Image area corresponds with computational domain Ω1, and thin lines on the left figure show the size of the computational domain Ω2 utilized in case 2.

The prior model π(x ) was constructed to model star-like smooth objects with smooth boundaries and three layers of variable thicknesses. The speed of sound inside each subregion was drawn from realistic probability distributions [56]. The outermost layer modeled skin, and its speed of sound was taken to be normally distributed with cskin ~ N (1590, 5 2). The middle layer modeled fat tissue with cfat ~ N (1470,1). For the innermost layer, we set cbreast ~ N (1530,10 2), and for the water subdomain, the speed of sound was cwater ~ N (1500, 0.12). We simulate uncertainty of choosing π(x ) by generating samples from models without tumors. No measurement noise was added to the computed signals. Four draws forming π(x ) are shown in Fig. 3. We computed C = 128 draws from π(x ) and D = 48 eigenvectors corresponding to the largest eigenvalues were used to compute the packed form of Γν. C. Case 1: Pure Discretization Errors In the first test case, the performance of BAE was examined for compensating errors generated by a coarse spatial and temporal discretization. The computation domain was the same as was used in the computation of the simulated measurements; that is, Ω1. As a reference, we also computed the reconstruction using dense spatial and temporal discretizations, and refer to this solution as an accurate solution. The estimates computed with

Fig. 4. Measured noiseless signals of accurate forward model, approximate forward model in case 1, and their difference.

the accurate solver used spatial discretization of 1794 × 1794 points and other parameters were the same as in the model used to simulate the data. The discretization is only slightly sparser than with the model that was used to compute the simulated measurements. The approximate solver Aλ0(x ) was discretized in a 466 × 466 grid with computation domain Ω1. To avoid the inverse crime, the IIS forward solver was not utilized here. Instead, the optimal nearly analytic discrete method (ONADM) [58] solver with 1024 time steps was used as the approximate forward solver. The approximation error samples were computed as explained in Sections IV-A and V-B. Examples of the measured noiseless signals that were computed with the accurate and approximate forward models, and their difference, are shown in Figs. 4 and 5. The difference signal is part of the approximation error vector for the particular target. The maximum error component here is around 2% of maximum signal value. The reconstructions are shown in Fig. 2 and cross-sections in Fig. 6. In addition, RMS errors and computing times are shown in Table I. The reconstruction time with the accurate solver is approximately 20 times higher than with the CEM and BAE. Reconstruction time with the BAE is only slightly higher than with CEM. The RMS error with the BAE is approximately the same as with the accurate solver, whereas the RMS error with the CEM model is about 40% larger than with the accurate model. Most importantly, of the three point-like inclusions around the larger one (shown in white), one can clearly be identified in the BAE estimate, whereas none can be in the CEM estimate. These inclusions correspond to increased wave speeds and potential tumors. One of the point-like inclusions is not visible even in the estimate obtained with the accurate solver. The starlike structure of the large white inclusion, which is typical for tumors, is clear in all reconstructions. In addition, we compared the influence of the number of samples C used to compute sample-based estimates for modeling error statistics. With C ≤ 40, reconstruction quality varied widely; however, even with small C values, the reconstruction result was in this case much better than with the method utilizing CEM.

1634

IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control ,

vol. 61, no. 10,

October

2014

Fig. 5. Simulated pressure field of accurate forward model (left), pressure field of approximate forward model in case 1 (middle), and their difference (right).

D. Case 2: Discretization and Absorbing Boundary Model Errors In the second case, reconstructions utilizing approximate model Aλ0(x ) were computed in the smaller computation domain Ω2 = 0.22 × 0.22 m; see Fig. 3. In addition, computationally lighter, but less accurate, first-order absorbing boundary conditions were used to cut the computation domain further. The boundaries of the computation domain were so close to the target and the receivers that reflections from the boundary cannot be taken as insignificant, as in can be seen in Figs. 7 and 8. In the computation of the second-order statistics of the approximation errors, Aλ(x ) was computed in Ω1 and Aλ0(x ) in Ω2. The approximation errors contained the contribution of errors that were due to both the numerical approximations and the reflections from the boundary of the computation domain. Examples of the measured noiseless signals that were computed with the accurate and approximate forward models, and their difference, are shown in Figs. 7 and 8. The difference signal is part of the approximation error vector for the particular target. The maximum error component here is around 40%. This is clearly larger than in the Case 1. The reconstructions are shown in Fig. 9, cross-sections in Fig. 10, and RMS errors and reconstruction computing times are presented in Table I. In this test case, conver-

Fig. 6. Cross-sections of the reconstructions in case 1. Figure is generated from the vertical line shown in the leftmost image of Fig. 2.

gence of the method utilizing CEM was random; in some steps the RMS error increased. However, the method utilizing BAE and the method utilizing with accurate forward model converged monotonically. At the same time, the posterior potential X decreased monotonically with all methods. The reconstruction time with the accurate solver was around 42 times higher than with the CEM and BAE. Reconstruction time with the BAE was only slightly higher than with CEM. The RMS error with the BAE was approximately same level as in the Case 1, approximately 8% higher than with the accurate solver, whereas the RMS error with the CEM model was about 90% larger than with accurate model. Overall, the BAE retained the tissue contrast better than the CEM. Of the three point-like inclusions around the larger one, two can be identified in the BAE estimate, whereas none can be identified in the CEM estimate. The larger inclusion near the bottom of the image is clearly visible in the accurate and BAE estimates, whereas it is quite faint in the CEM estimate. The same applies for the shape of the large star-like inclusion. Influence of number of samples C used to obtain sample based estimates for modeling error statistics was studied. However, in this case, RMS error with C = 2 was significantly higher than with a higher number of samples. Similar variation of RMS error with small number of samples can be seen as in test case 1. Again, the method with BAE

Fig. 7. Measured noiseless signals of accurate forward model, approximate forward model in case 2, and their difference.

koponen et al.: bayesian approximation error approach in full-wave ultrasound tomography

1635

Fig. 8. Simulated pressure field of accurate forward model (left), pressure field of approximate forward model in case 2 (middle), and their difference (right). The difference between the fields is significantly higher than in case 1.

Fig. 9. Speed of sound reconstructions of simulations in case 2. Images from left to right are the numerical model, reconstruction with the approximate forward solver utilizing CEM, reconstruction with the approximate forward solver utilizing BAE, and reconstruction with the accurate forward solver.

TABLE I. Parameters, Reconstruction Times, and Reconstruction Errors for the Conventional Error Model (CEM) and the Bayesian Approximation Error Model (BAE). Solver Measurements Accurate CEM BAE CEM BAE

K

N

N/λ

Stability factor

Ω

Scheme

4096 4096 1024 1024 1024 1024

1858 1794 466 466 306 306

52.5 50.7 13.2 13.2 13.0 13.0

0.52 0.50 0.52 0.52 0.51 0.51

Ω1 Ω1 Ω1 Ω1 Ω2 Ω2

IIS IIS ONADM ONADM ONADM ONADM

Reconstruction time 31 h 1h 1h

40 33 34 45 46

min min min min min

3 59 26 4 52

RMS error s s s s s

686 942 732 1293 738

Columns from left to right are forward solver name, number of discretization steps in temporal domain K, number of discretization points in one dimension in spatial domain N, number of discretization points per wavelength frequency of 150 kHz N/λ, stability factors cmax(Δt/Δx), the computational domains Ω, discretization schemes, reconstruction times, and the RMS error of reconstructions.

gave significantly higher quality reconstructions than the method with CEM, even with small C values. VI. Discussion

Fig. 10. Cross-sections of the reconstructions in case 2. Image is generated from the vertical line shown in the leftmost image of Fig. 9.

The use of approximate solvers for the ultrasound tomography problem of reconstructing the spatial distribution of the speed of sound was considered. With approximate solvers, the computation times can be reduced by at least an order of magnitude. The reconstructions with BAE retain the contrast better than the CEM, and allow for detection of small inclusions better. However, it was shown that the adoption of the Bayesian approximation error approach increases

1636

IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control ,

computational load only slightly when compared with the conventional error model (CEM). For a given geometry and transmitter–receiver configuration, the computationally heavy tasks in BAE approach have to be carried out only once; these computations consist of the solution of the forward problem only. All inversions are carried out with the fast approximate solvers.

References [1] A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm,” Ultrason. Imaging, vol. 6, no. 1, pp. 81–94, 1984. [2] K. A. Dines and R. J. Lytle, “Computerized geophysical tomography,” Proc. IEEE, vol. 67, no. 7, pp. 1065–1073, Jul. 1979. [3] R. Gordon, R. Bender, and G. T. Herman, “Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography,” J. Theor. Biol., vol. 29, no. 3, pp. 471–481, 1970. [4] A. Hormati, I. Jovanović, O. Roy, and M. Vetterli, “Robust ultrasound travel-time tomography using the bent ray model,” in SPIE Med. Imag., 2010, art. no. 76290I. [5] A. J. Devaney, “Geophysical diffraction tomography,” IEEE Trans. Geosci. Rem. Sens., vol. GE-22, no. 1, pp. 3–13, Jan. 1984. [6] J. F. Greenleaf, “Computerized tomography with ultrasound,” Proc. IEEE, vol. 71, no. 3, pp. 330–337, Mar. 1983. [7] N. Duric, P. Littrup, A. Babkin, D. Chambers, S. Azevedo, A. Kalinin, R. Pevzner, M. Tokarev, E. Holsapple, O. Rama, and R. Duncan, “Development of ultrasound tomography for breast imaging: Technical assessment,” Med. Phys., vol. 32, no. 5, pp. 1375–1386, 2005. [8] O. Gauthier, J. Virieux, and A. Tarantola, “Two-dimensional nonlinear inversion of seismic waveforms; Numerical results,” Geophysics, vol. 51, no. 7, pp. 1387–1403, 1986. [9] I. Jovanovic, “Inverse problems in acoustic tomography: Theory and applications,” Ph.D. thesis, School of Computer and Communication Sciences, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, 2008. [10] O. Roy, I. Jovanovic, A. Hormati, R. Parhizkar, and M. Vetterli, “Sound speed estimation using wave-based ultrasound tomography: Theory and GPU implementation,” in Society of Photo-Optical Instrumentation Engineers (SPIE) Conf. Series, 2010, vol. 7629, art. no. 76290J. [11] A. Tarantola, “Inversion of seismic reflection data in the acoustic approximation,” Geophysics, vol. 49, no. 8, pp. 1259–1266, 1984. [12] A. Fichtner, H. P. Bunge, and H. Igel, “The adjoint method in seismology: I. Theory,” Phys. Earth Planet. Inter., vol. 157, no. 1–2, pp. 86–104, 2006. [13] R.-E. Plessix, “A review of the adjoint-state method for computing the gradient of a functional with geophysical applications,” Geophys. J. Int., vol. 167, no. 2, pp. 495–503, 2006. [14] A. Tarantola, “Theoretical background for the inversion of seismic waveforms including elasticity and attenuation,” Pure Appl. Geophys., vol. 128, no. 1, pp. 365–399, 1988. [15] D. Calvetti and E. Somersalo, An Introduction to Bayesian Scientific Computing—Ten Lectures on Subjective Computing. New York, NY: Springer, 2007. [16] J. P. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems. New York, NY: Springer-Verlag, 2005. [17] A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation. Philadelphia, PA: SIAM, 2004. [18] J. P. Kaipio and V. Kolehmainen, “Approximate marginalization over modeling errors and uncertainties in inverse problems,” in P. Damien, N. Polson, and D. Stephens, Eds., Bayesian Theory and Applications, New York, NY: Oxford University Press, 2013, pp. 644–672. [19] J. Kaipio and E. Somersalo, “Statistical inverse problems: Discretization, model reduction and inverse crimes,” J. Comput. Appl. Math., vol. 198, no. 2, pp. 493–504, 2007. [20] S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors

vol. 61, no. 10,

October

2014

and model reduction with an application in optical diffusion tomography,” Inverse Probl., vol. 22, no. 1, pp. 175–195, 2006. [21] J. Heino, E. Somersalo, and J. P. Kaipio, “Compensation for geometric mismodelling by anisotropies in optical tomography,” Opt. Express, vol. 13, no. 1, pp. 296–308, 2005. [22] J. Heino and E. Somersalo, “A modelling error approach for the estimation of optical absorption in the presence of anisotropies,” Phys. Med. Biol., vol. 49, no. 20, pp. 4785–4798, 2004. [23] V. Kolehmainen, M. Schweiger, I. Nissilä, T. Tarvainen, S. R. Arridge, and J. P. Kaipio, “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A, vol. 26, no. 10, pp. 2257–2268, 2009. [24] M. Mozumder, T. Tarvainen, S. R. Arridge, J. Kaipio, and V. Kolehmainen, “Compensation of optode sensitivity and position errors in diffuse optical tomography using the approximation error approach,” Biomed. Opt. Express, vol. 4, no. 10, pp. 2015–2031, 2013. [25] T. Tarvainen, A. Pulkkinen, B. T. Cox, J. P. Kaipio, and S. R. Arridge, “Bayesian image reconstruction in quantitative photoacoustic tomography,” IEEE Trans. Med. Imaging, vol. 32, no. 12, pp. 2287–2298, 2013. [26] D. Calvetti, J. P. Kaipio, and E. Somersalo, “Aristotelian prior boundary conditions,” Int. J. Math., vol. 1, pp. 63–81, 2006. [27] A. Lehikoinen, S. Finsterle, A. Voutilainen, L. M. Heikkinen, M. Vauhkonen, and J. P. Kaipio, “Approximation errors and truncation of computational domains with application to geophysical tomography,” Inverse Probl. Imaging, vol. 1, no. 2, pp. 371–389, 2007. [28] A. Nissinen, L. M. Heikkinen, and J. P. Kaipio, “The Bayesian approximation error approach for electrical impedance tomography— Experimental results,” Meas. Sci. Technol., vol. 19, no. 1, art. no. 015501, 2008. [29] A. Nissinen, L. M. Heikkinen, V. Kolehmainen, and J. P. Kaipio, “Compensation of errors due to discretization, domain truncation and unknown contact impedances in electrical impedance tomography,” Meas. Sci. Technol., vol. 20, no. 10, art. no. 105504, 2009. [30] A. Nissinen, V. Kolehmainen, and J. P. Kaipio, “Compensation of modelling errors due to unknown domain boundary in electrical impedance tomography,” IEEE Trans. Med. Imaging, vol. 30, no. 2, pp. 231–242, 2011. [31] A. Nissinen, V. Kolehmainen, and J. P. Kaipio, “Reconstruction of domain boundary and conductivity in electrical impedance tomography using the approximation error approach,” Int. J. Uncertain. Quantif., vol. 1, no. 3, pp. 203–222, 2011. [32] J. M. J. Huttunen and J. P. Kaipio, “Approximation errors in nonstationary inverse problems,” Inverse Probl. Imaging, vol. 1, no. 1, pp. 77–93, 2007. [33] J. M. J. Huttunen and J. P. Kaipio, “Approximation error analysis in nonlinear state estimation with an application to state-space identification,” Inverse Probl., vol. 23, no. 5, pp. 2141–2157, 2007. [34] J. M. J. Huttunen and J. P. Kaipio, “Model reduction in state identification problems with an application to determination of thermal parameters,” Appl. Numer. Math., vol. 59, no. 5, pp. 877–890, 2009. [35] J. M. J. Huttunen, A. Lehikoinen, J. Hämäläinen, and J. P. Kaipio, “Importance filtering approach for the nonstationary approximation error method,” Inverse Probl., vol. 26, no. 12, art. no. 125003, 2009. [36] A. Lehikoinen, J. M. J. Huttunen, A. Voutilainen, S. Finsterle, M. B. Kowalsky, and J. P. Kaipio, “Dynamic inversion for hydrological process monitoring with electrical resistance tomography under model uncertainties,” Water Resour. Res., vol. 46, no. 4, art. no. W04513, 2010. [37] T. Lähivaara, N. F. Dudley Ward, T. Huttunen, J. Koponen, and J. P. Kaipio, “Estimation of aquifer dimensions from passive seismic signals with approximate wave propagation models,” Inverse Probl., vol. 30, no. 1, art. no. 015003, 2014. [38] T. Tarvainen, V. Kolehmainen, A. Pulkkinen, M. Vauhkonen, M. Schweiger, S. R. Arridge, and J. P. Kaipio, “Approximation error approach for compensating modelling errors between the radiative transfer equation and the diffusion approximation in diffuse optical tomography,” Inverse Probl., vol. 26, no. 1, art. no. 015005, 2010. [39] T. Tarvainen, V. Kolehmainen, J. P. Kaipio, and S. R. Arridge, “Corrections to linear methods for diffuse optical tomography using approximation error modelling,” Biomed. Opt. Express, vol. 1, no. 1, pp. 209–222, 2010. [40] V. Kolehmainen, T. Tarvainen, S. R. Arridge, and J. P. Kaipio, “Marginalization of uninteresting distributed parameters in inverse problems—Application to diffuse optical tomography,” Int. J.Uncertain. Quantif., vol. 1, no. 1, pp. 1–17, 2011.

koponen et al.: bayesian approximation error approach in full-wave ultrasound tomography [41] A. Pulkkinen, V. Kolehmainen, J. P. Kaipio, B. T. Cox, S. R. Arridge, and T. Tarvainen, “Approximate marginalization of unknown scattering in quantitative photoacoustic tomography,” Inverse Probl. Imaging, to be published. [42] M. Sambridge, “Geophysical inversion with a neighbourhood algorithm—I. Searching a parameter space,” Geophys. J. Int., vol. 138, no. 2, pp. 479–494, 1999. [43] U. Meier, A. Curtis, and J. Trampert, “Global crustal thickness from neural network inversion of surface wave data,” Geophys. J. Int., vol. 169, no. 2, pp. 706–722, 2007. [44] J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys., vol. 114, no. 2, pp. 185–200, 1994. [45] Q. H. Liu and J. Tao, “The perfectly matched layer for acoustic waves in absorptive media,” J. Acoust. Soc. Am., vol. 102, no. 4, pp. 2072–2082, 1997. [46] B. Engquist and A. Majda, “Absorbing boundary conditions for numerical simulation of waves,” Proc. Natl. Acad. Sci. USA, vol. 74, no. 5, pp. 1765–1766, 1977. [47] Y. Liu, L. Ding, and M. K. Sen, “Comparisons between the hybrid ABC and the PML method for 2D high-order finite-difference acoustic modeling,” in SEG Tech. Prog. Expanded Abstracts 2011, pp. 2952–2956. [48] J. F. Semblat, L. Lenti, and A. Gandomzadeh, “A simple multidirectional absorbing layer method to simulate elastic wave propagation in unbounded domains,” Int. J. Numer. Methods Eng., vol. 85, no. 12, pp. 1543–1563, 2011. [49] M. A. Woodbury, “Inverting modified matrices,” Princeton University, Statistical Research Group, Princeton, NJ, Memorandum Report, vol. 42, p. 4, 1950. [50] M. Brand, “Incremental singular value decomposition of uncertain data with missing values,” in Eur. Conf. Computer Vision, 2002, pp. 707–720. [51] M. Fink, “Time reversal of ultrasonic fields. I. Basic principles,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 39, no. 5, pp. 555–566, 1992. [52] J. Kiefer, “Sequential minimax search for a maximum,” Proc. Am. Math. Soc., vol. 4, pp. 502–506, 1953. [53] NVIDIA Corporation. About CUDA. https://developer.nvidia. com/about-cuda, 2012. [54] J. Sanders and E. Kandrot, CUDA by Example: An Introduction to General-Purpose GPU Programming. Reading, MA: Addison-Wesley Professional, 2010. [55] D. Michéa and D. Komatitsch, “Accelerating a three-dimensional finite-difference wave propagation code using GPU graphics cards,” Geophys. J. Int., vol. 182, no. 1, pp. 389–402, 2010. [56] E. Franceschini, S. Mensah, D. Amy, and J. P. Lefebvre, “A 2-D anatomic breast ductal computer phantom for ultrasonic imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 53, no. 7, pp. 1281–1288, 2006. [57] K. Kowalczyk and M. Van Walstijn, “Wideband and isotropic room acoustics simulation using 2-D interpolated FDTD schemes,” IEEE Trans. Audio Speech Lang. Process., vol. 18, no. 1, pp. 78–89, 2010. [58] D. Yang, M. Lu, R. Wu, and J. Peng, “An optimal nearly analytic discrete method for 2D acoustic and elastic wave equations,” Bull. Seismol. Soc. Am., vol. 94, no. 5, pp. 1982–1992, 2004.

1637

Janne Koponen received the M.Sc degree in information technology (simulation and optimization) from the University of Jyväskylä in 2007. Since 2010, he has worked as a junior researcher in the University of Eastern Finland. Since 2013, Janne has also worked as technology manager in Rocsole Ltd. His research interests include Bayesian methods with acoustic tomography, especially methods with industrial applications.

Tomi Huttunen received the M.Sc. and Ph.D. degrees in medical physics from the University of Kuopio, Finland, in 2001 and 2004, respectively. He has worked as a visiting researcher in the University of Delaware from 1999 to 2000; and in the Harvard Medical School in 2001 to 2002. Since 2007, Tomi Huttunen has been an adjunct professor in the University of Eastern Finland. He is a co-founder of Kuava Ltd. Currently, he is Kuava’s director of R\&D and has a part-time position in the University of Eastern Finland. His research interests are computational acoustics and ultrasonics.

Tanja Tarvainen received the M.Sc. and Ph.D. degrees from the University of Kuopio, Kuopio, Finland. She is currently Academy Research Fellow in the Department of Applied Physics, University of Eastern Finland, Finland, and part-time Research Associate at the Department of Computer Science, University College London, UK. Her research interests include Bayesian approaches to inverse problems, with applications especially in optical and acoustic tomographic methods.

Jari P. Kaipio received the M.Sc. and Ph.D. degrees from the University of Kuopio, Kuopio, Finland. He is currently a Professor of applied mathematics in the Department of Mathematics, University of Auckland, Auckland, New Zealand. He also leads the Inverse Problems Research Group, University of Eastern Finland, Kuopio, as a part-time Professor of computational physics. He has authored more than 140 papers, mostly in inverse problems, as well as the book Statistical and Computational Inverse Problems with Erkki Somersalo. His current research interests include inverse problems, especially Bayesian methods and nonstationary problems with various application areas. He serves on the editorial board of Inverse Problems, International Journal for Uncertainty Quantification, and Inverse Problems and Imaging. Prof Kaipio is a Fellow of the Institute of Physics, UK.