Blood perfusion estimation in heterogeneous tissue using BEM based algorithm

Engineering Analysis with Boundary Elements 39 (2014) 75–87

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Blood perfusion estimation in heterogeneous tissue using BEM based algorithm J. Iljaž n, L. Škerget Faculty of Mechanical Engineering, University of Maribor, Smetanova 17, SI-2000 Maribor, Slovenia

art ic l e i nf o

a b s t r a c t

Article history: Received 30 June 2012 Accepted 2 November 2013 Available online 30 November 2013

The estimation of space-dependent perfusion coefficient in homogeneous and non-homogeneous tissue has been investigated. While initial and Dirichlet boundary conditions are known, additional heat-flux measurement data is needed to render a unique solution. A numerical approach based on Boundary Element Method (BEM) combined with two different optimization routines and first-order Tikhonov regularization using L-Curve method has been developed. Efficiency of the algorithm, effect of initial guess, noise, perfusion distribution and non-homogeneous tissue on retrieving the perfusion coefficient has been studied on two test examples using exact as well as noisy data. Results show very good agreement with the true perfusion function under exact and low-noisy data, using Levenberg–Marquardt (LM) method combined with first-order regularization process. If the true perfusion function firstderivative is large, the function can be successfully retrieved only in the region near the boundary measurement, which is especially noticeable for a non-monotonic function in non-homogeneous tissue. This study represents the base for further research on the field of successful non-invasive blood perfusion determination in non-homogeneous tissue. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Bio-heat Inverse problem Boundary Element Method Optimization method Space-dependent perfusion Non-homogeneous tissue

1. Introduction Blood perfusion is an important clinical parameter for disease diagnostics, drug delivery, cancer hyperthermia, cryosurgery, etc. [1–4]. It is characterized as the blood flow rate through the capillary network, small arterioles, venules and extracellular space. For this reason the blood perfusion is multidirectional at the local, microscopic level, however from the global or continuum perspective it is directionless and defined as the volumetric flow rate per volume of tissue. It plays an important role in tissue physiology for the local transportation of oxygen, nutrition and waste products, as well it is the main attribute in thermoregulatory process. Therefore, its deviation from normal level indicates abnormal physiologic or pathologic condition of the tissue, like high vascularized tumour, tissue inflammation, thrombosis, vascular sclerosis, burn injuries. This deviation is also reflected in the temperature change of the tissue [1]. Thus, measurement techniques for blood perfusion represent a valuable medical diagnostic tool. Among various techniques like Magnetic Resonance Imaging (MRI), Positron Emission Tomography (PET), Laser Doppler Flowmetry (LDF) and Photoplethysmography (PPG), the

n

Corresponding author. Tel.: þ 386 2 220 7784; fax: þ386 2 220 7990. E-mail addresses: [email protected], [email protected] (J. Iljaž), [email protected] (L. Škerget). 0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.11.002

Thermal Clearance Techniques (TCT) [5] are showing a promising way due to their non-invasiveness, economical, quantitativeness, safe usage, as well as applicability to a wider range of clinical applications [6–10]. The TCT are based on the bio-heat equation, which correlates blood perfusion with tissue temperature, and known variables like tissue temperature, heat flux or power, which can be measured invasively or non-invasively. Since the measurement data are indirectly correlated with the unknown variables, in this case the blood perfusion, the problem is said to be inverse. Most important researches in blood perfusion determination [6–9,11–13] have been focused on techniques to measure temperature, heat flux or its response, combined with different approaches to estimate constant or time-dependent blood perfusion in homogeneous tissue non-invasively, based on Pennes bio-heat equation. These studies still lack in more realistic problem description (considering non-homogeneous tissue) and in evaluation of numerical and optimization methods together with problem limitations. Important steps to comprehend non-homogeneous blood perfusion have been presented by Partridge and Wrobel [2], Erhart et al. [14] and Trucu et al. [15] with temperature-dependent blood perfusion and Trucu et al. [16,17] with space-dependent blood perfusion using a numerical approach based on finite differences (FD) and NAG routine optimization. Studies of Trucu et al. [15–19] are mainly dealing with sufficient additional measurements (invasive and noninvasive) to render a unique solution of the inverse bio-heat problem

76

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

together with few numerical examples. Until now, tissue has been treated only as homogeneous with constant material properties; therefore this study presents a step further to estimate spacedependent blood perfusion in a non-homogeneous tissue for diagnostic purpose. This paper investigates the inverse bio-heat problem of estimating space-dependent blood perfusion in known non-homogeneous tissue. To solve the problem a numerical approach has been developed based on Boundary Element Method (BEM) and optimization algorithm. For optimization the Broyden–Fletcher–Goldfarb–Shanno (BFGS) and Levenberg–Marquardt (LM) method coupled with firstorder Tikhonov regularization combined with L-Curve method have been used. Only the Dirichlet and Neumann boundary data are taken into consideration. The reason lies in the possibility to measure the needed data non-invasively, which is the tendency of modern medicine. Trucu et al. [17] also discuss some other inverse problems to estimate space-dependent perfusion, as well as Rodrigues et al. [20] and Yang et al. [21] for general parabolic equations. These problems would need invasive measurements of temperature which can lead to discomfort, pain or even inflammation. Therefore, this type of problem has been left out. Blood perfusion and material properties have been treated as space-dependent, which can lead into treating different types of tissue like skin, muscle, fat, etc. For this reason, this study represents the base for further research on solving inverse bio-heat problems, as well as the base for development of new improved non-invasive techniques to measure space-dependent blood perfusion in real tissue for diagnostic purpose. The paper is organized as follows. In Section 2 the mathematical framework of the investigated inverse bio-heat problem is presented. The numerical approach based on BEM together with BFGS and LM regularized optimization method is given in Section 3. Section 4 discusses some examples used to test the proposed numerical approach and inverse problem limitations, using exact and noisy measurement data. In Section 5 the results and discussion of test examples are presented, followed by conclusions in Section 6.

2. Mathematical formulation The proposed numerical algorithm has been developed for solving an inverse bio-heat problem in two dimensions; for this reason all the on-going formulation is written in two-dimensional form. The base of the inverse problem represents a mathematical model for the bio-heat equation which correlates unknown variable such as blood perfusion with measurement data like temperature or heat flux. The solution of investigated inverse bio-heat problem in this study has been established on the Pennes bio-heat equation [22]:

 ∂T ∂ ∂T ∂ ∂T k þ k þ ωb ρb cb ðT  T a Þ þ I m ; ρc ¼ ð1Þ ∂t ∂x ∂x ∂y ∂y where ωb, ρb and cb are blood perfusion, density and specific heat of the blood, respectively, ρ, c and k are density, specific heat and thermal conductivity of the tissue, T is temperature of the tissue, Ta is arterial blood temperature and Im the metabolic heat generation. Pennes' model describes the tissue temperature as a function of blood perfusion, cell metabolic heat generation and heat exchange between the tissue and its surroundings. Through the years there were several proposals and improvements to describe bio-heat transfer better [23–25], but Pennes' model is still extensively used, due to its simplicity and good agreement with the experiments. The governing equation (1) has been treated and solved in its dimensionless form [26], written as 1 ∂θ ∂2 θ ∂2 θ 1 I ¼ þ  P θþ ; a ∂t ∂x2 ∂y2 k f k

ð2Þ

where θ, t, a, k and I are now dimensionless temperature, time, thermal diffusivity, thermal conductivity and metabolic heat generation, respectively. We omit the asterisk on the dimensionless variables like thermal diffusivity, time due to more clear writing. Now the blood perfusion is hidden inside the perfusion coefficient, Pf, which is defined as Pf ¼

ωb ρb c b L 2 k0

;

ð3Þ

where L and k0 represent characteristic length and thermal conductivity, respectively. Blood perfusion can be positive or equal to zero (ωb Z 0), therefore also the value of the perfusion coefficient has the same tendency (P f Z 0). Dimensionless temperature has been defined as θ ¼ ðT  T a Þ=ΔT, where ΔT ¼ ðI m0 L2 Þ=k0 and Im0 represent the characteristic metabolic heat generation. Non-dimensional time has been defined as t ¼ ða0 tÞ=L2 , where a0 ¼ k0 =ðρ0 c0 Þ represents the characteristic thermal diffusivity, ρ0 the characteristic density and c0 the characteristic specific heat. For the case of constant material properties, the characteristic variables are chosen in such manner such that k ¼ 1. Consequently, Eq. (2) can be rewritten as ∂θ ∂2 θ ∂2 θ ¼ þ  P f θ þ I; ∂t ∂x2 ∂y2

ð4Þ

as can be found in [26,17]. The above equation represents a second-order parabolic equation that can be also found in some other physical applications, like optical tomography, ground water flow [20,27]. Blood perfusion and tissue properties have been treated as space dependent variables; ωb ¼ ωb ðx; yÞ; a ¼ aðx; yÞ, where properties are treated as piecewise constant. For this reason the computational domain Ωðx; yÞ was divided into D subdomains ΩS ðx; yÞ, where constant material properties have been considered, aS ¼ const: To solve Eq. (2) for the whole computational domain we have to write it for each subdomain and merge the governing equations with compatibility θ1 ðr I ; tÞ ¼ θ2 ðr I ; tÞ and equilibrium  k1 ∂θ1 ðr I ; tÞ=∂νI ¼ k2 ∂θ2 ðr I ; tÞ=∂νI interface conditions, where indices 1 and 2 stand for the first and its adjacent or second subdomain, respectively, r I ¼ ðxI ; yI Þ is the position of the interface and νI its normal. Inverse bio-heat problem of estimating blood perfusion has now been transformed to a problem of estimating the perfusion coefficient based on the known initial and Dirichlet boundary conditions and additional measurement of heat flux, which is needed to render a unique solution. Therefore, perfusion coefficient and dimensionless temperature inside the computational domain are unknown functions, which we seek to retrieve. It is known that the problem is ill-posed and therefore challenging both theoretically and numerically. For this reason we restrict the blood perfusion dependency only in one direction, P f ¼ P f ðxÞ, for which we were able to find a uniqueness theorem stated below. The inverse bio-heat problem of space-dependent perfusion coefficient estimation can be now mathematically formulated as follows: Find the dimensionless temperature θ ¼ θðx; y; tÞ such that θ; ∂θ=∂x; ∂θ=∂y A Cð½0; L  ½0; H  ½0; 1ÞÞ, ∂θ=∂t; ∂2 θ=∂x2 A Cðð0; LÞ  ð0; HÞ  ð0; 1ÞÞ and the perfusion coefficient P f ¼ P f ðxÞ A Cð½0; LÞ, P f ðxÞ 4 0 in the computational domain ðx; y; tÞ A ½0; L  ½0; H ½0; 1Þ, satisfying system of Eqs. (2) written for each subdomain ΩS with constant properties; a; k A R þ , and coupled with compatibility and equilibrium interface conditions, subject to the initial condition:

θðx; y; 0Þ ¼ u0 ðx; yÞ;

ðx; yÞ A ½0; L  ½0; H;

ð5Þ

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

Dirichlet boundary conditions:

θð0; y; tÞ ¼ g0 ðtÞ;

ðy; tÞ A ½0; H  ½0; 1Þ;

ð6Þ

θðL; y; tÞ ¼ g1 ðtÞ;

ðy; tÞ A ½0; H  ½0; 1Þ;

ð7Þ

77

solution and integrating by parts, the boundary-domain integral form is obtained: Z Z cðξÞθðξ; tÞ þ θðR; tÞqn ðξ; RÞ dΓ S ¼ qðR; tÞθn ðξ; RÞ dΓ S ΓS

ΓS

Z

þ

Neumann boundary conditions: ∂θ=∂yðx; 0; tÞ ¼ h0 ðtÞ;

ðx; tÞ A ½0; L  ½0; 1Þ;

ð8Þ

∂θ=∂yðx; H; tÞ ¼ h2 ðtÞ;

ðx; tÞ A ½0; L  ½0; 1Þ;

ð9Þ

and additional Neumann data at the boundary x ¼ L: ∂θ=∂xðL; y; tÞ ¼ h1 ðtÞ;

ðy; tÞ A ½0; L  ½0; 1Þ:

ð10Þ

We formulate the theorem of uniqueness for the inverse problem given by Ramm [28] and Denisov [29, pp. 139–146] in two-dimensions, as follows. Theorem 1. Let u0 ¼ g 0 ¼ h0 ¼ h2 ¼ I ¼ 0 and g 1 ðtÞ a0 satisfy the following conditions; g 1 A C 2 ð½0; 1Þ, g 1 ð0Þ ¼ g ′1 ð0Þ ¼ 0 and there exists t 0 4 0 such that g 1 ðtÞ ¼ 0 for all t Z t 0 . Then the additional Neumann data h1 ðtÞ renders a unique solution ðP f ðxÞ; θðx; y; tÞÞ of the inverse problem (4)–(10). 3. Numerical approach

3.1. BEM The BEM is based on the integral transformation of the weighted governing equations with the fundamental solution. The reason for choosing this method is in the implicit treatment of Neumann boundary condition and its value (heat flux) for solving inverse bioheat problem. The normal derivative of the field function (temperature) on the boundary is directly incorporated into the algebraic system of equations, for which we do not need to use any approximation for its calculation from the field function. Therefore, it is defined and calculated more accurately. To derive the BEM for direct bio-heat problem, the Pennes bio-heat equation (2) written for subdomain ΩS has been treated as an inhomogeneous elliptic Poisson's equation as L½θðr; tÞ þbðr; tÞ ¼ 0;

E

cðξÞθðξ; tÞ þ ∑ fhg fθðtÞg ¼ ∑ fggtr fqðtÞgn tr

n

e¼1

e¼1

C

n

þ ∑ fsgtr fbðtÞg ;

ð13Þ

c¼1

where E and C are the number of boundary and domain elements, respectively, and n the number of discrete points inside the boundary or domain element. Elements of the vector {h}, {q} and {s} are R R n n n n consequently he ¼ Γ Se Φ qn ðξ; RÞ dΓ S , g ne ¼ Γ Se Φ θ ðξ; RÞ dΓ S and R n n n n sc ¼ ΩSc Φ θ ðξ; rÞ dΩS , where Φ represents the n-th interpolation function. The above discrete equation has been applied to all mesh nodal points of the subdomain (boundary and domain) from where the system of equations written in matrix form is obtained: ½HfθðtÞg ¼ ½GfqðtÞg þ ½SfbðtÞg;

ð14Þ

where [H], [G] and [S] are now fully populated matrixes, fθðtÞg, fqðtÞg and fbðtÞg are vectors containing the nodal value of field function, its normal derivative and the source term, respectively. The inhomogeneous term bðr; tÞ is composed of source term and time derivative and is defined as bðr; tÞ ¼ 

1 ∂θðr; tÞ 1 I  P f ðxÞθðr; tÞ þ : a ∂t k k

ð15Þ

The time derivative term has been treated with second-order asymmetric finite difference scheme: ∂θðr; t i Þ 3θðr; t i Þ  4θðr; t i  1 Þ þ θðr; t i  2 Þ  ; ∂t 2dt

ð16Þ

where index i represents the time step and dt the time difference between two adjacent time steps. Therefore, term bðr; tÞ using fully implicit time scheme, θðr; tÞ ¼ θðr; t i Þ, can be written as bðr; t i Þ ¼ α1 θðr; t i Þ þ α2 θðr; t i  1 Þ þ α3 θðr; t i  2 Þ

or


 ∂2 θðr; tÞ ∂2 θðr; tÞ 1 ∂θðr; tÞ 1 I  P ¼ 0; þ þ  ðxÞ θ ðr; tÞ þ f a ∂t k k ∂x2 ∂y2

ð12Þ

where ξ ¼ ξðx; yÞ represents the position of the source point, Γ S is the boundary of subdomain ΩS , qðR; tÞ ¼ ∂θðR; tÞ=∂ν denotes the normal derivative of dimensionless temperature, R ¼ Rðx; yÞ is the position n vector of the boundary point, θ ðξ; rÞ ¼ 1=ð2π Þ  ln½1=dðξ; rÞ, n qn ðξ; RÞ ¼ ∂θ ðξ; RÞ=∂ν are the Laplace fundamental solution and its normal derivative in two-dimensions, respectively, where dðξ; rÞ represents the distance between the source and arbitrary point and cðξÞ is the geometrical free coefficient, which depends on the spatial position of the source point and has been determined using rigid body movement. To solve the integral equation (12) a numerical approach is used. Therefore, the computational subdomain ΩS , as well as its boundary Γ S , has to be discretized with boundary and domain elements. Integrals are then approximated with the sum of integrals over those tr elements. Introducing the interpolation functions fΦg for dimensional temperature and its normal derivative on the boundary and source inside the domain we obtain the discrete form of Eq. (12) as E

To solve the inverse problem a numerical approach based on the BEM and BFGS, as well as LM optimization method coupled with first-order Tikhonov regularization, has been used. The investigated problem has been transformed to an optimization one with defined objective function FðP f Þ, which we try to minimize. Consequently, the algorithm is composed of two steps. The first step represents solving the direct problem (4)–(9) using the BEM with known initial, Dirichlet and Neumann boundary conditions and assumed value of perfusion function, Pf, to determine the value of objective function, while the second step represents the minimization of the objective function in a manner that FðP f ;k þ 1 Þ o FðP f ;k Þ, where k ¼ 0; 1; 2; … represents an iteration step of the optimization process. In the following, a short description of the proposed algorithm is presented.

bðr; tÞθ ðξ; rÞ dΩS ; n

ΩS

þ α4 P f ðxÞθðr; t i Þ  α4 IðrÞ; ð11Þ

where LðÞ represents the Laplace differential operator, b is the inhomogeneous part, which incorporates all non-homogeneous effects like source term, as well as the time derivative term, and r ¼ rðx; yÞ the position of an arbitrary point in the subdomain. Multiplying Eq. (11) with the two-dimensional Laplace fundamental

ð17Þ

where α1 ¼ 3=ð2adtÞ, α2 ¼ 2=ðadtÞ, α3 ¼  1=ð2adtÞ and α4 ¼  1=k. Introducing Eq. (17) in (14) and using a fully implicit time scheme the system of equations can be rewritten as ð½H  α1 ½S  α4 fP f gtr ½SÞfθðt i Þg ¼ ½Gfqðt i Þg  α4 ½SfIg þ α2 ½Sfθðt i  1 Þg þ α3 ½Sfθðt i  2 Þg;

ð18Þ

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

78

which represents the linear system of equations that is solved using only one iteration for each time step. There is no need of using nonlinear iteration loop. Therefore, the computational time for direct problem solving is small, which leads to a more efficient and faster algorithm overall, due to the fact that optimization process needs to evaluate the objective function many times. For the direct bio-heat problem a continuous quadratic interpolation function (9 cell nodes) has been used for the field function θðr; tÞ and the source part I, a discontinuous constant interpolation function for the boundary value of ∂θðR; tÞ=∂ν and a continuous linear interpolation function (4 cell nodes) for the perfusion coefficient P f ðxÞ. Therefore, fP f g in Eq. (18) represents the vector of discrete values of the perfusion function at the mesh nodal points, which has also been used in the optimization process. The proposed numerical (18) has been applied to each subdomain: ΩS ; S ¼ 1; …; D, and then connected through the compatibility and equilibrium interface conditions described above to generate a global system of equations. The number of equations remains the same, only the number of unknowns is reduced by the number of interface unknowns [30,31]. Therefore, the final global system is overdetermined, which solves the direct bioheat problem in piecewise homogeneous tissue. Once the global initial and boundary conditions are applied, the system can be solved, where a LSQR solver has been used due to the overdetermination [32]. The presented BEM scheme using mixed elements has been developed to solve direct two-dimensional transient bio-heat problem in homogeneous, as well as piecewise homogeneous tissue, and to determine the objective function value that is needed for the optimization process. This BEM has also been analysed and tested on analytical examples to assure its correct implementation, solving and stability. Relatively coarse mesh is needed to attain the analytical solution, due to the used quadratic interpolation function, as well as because the direct bio-heat problem is well-posed. Here the BEM scheme analysis is not presented in detail, because solving the direct bio-heat problem is not the main aim of this paper.

3.2. Optimization method In the second step of the proposed numerical approach the first-order Tikhonov regularization method has been used, due to the fact that the inverse problem is ill-posed, especially in the presence of noisy data as it has been already pointed out by Trucu et al. [17]. For this reason, the regularized objective function F R ðfP f g; λÞ-R þ which we try to minimize has been defined as the difference between calculated Neumann boundary value at boundary x ¼ L with assumed perfusion function or its discrete value vector fP f g and the additional boundary measurement h1 ðtÞ:
2   1 m ∂θ   1 p ð P f ; ðL; H=2; t i ÞÞ h1 ðt i Þ þ λ ∑ ð½TfP f gj Þ2 F R ð Pf Þ ¼ ∑ 2 i ¼ 1 ∂x 2 j¼1   1   1 ¼ f ðfP f gÞtr f ð P f Þ þ λð½TfP f gÞtr ð½T P f Þ; 2 2

ð19Þ

where f ðfP f gÞ represents the residual vector, m the number of time measurements, λ the regularization parameter, p the number of mesh nodal points in the x-direction and [T] is a regularization matrix [33]. The first part of regularized objective function (19) represents the residual value or classical objective function, while the second one the used first-order Tikhonov regularization, which penalizes a high value of perfusion function first-derivative,

such that 2  1=Δx1 6 … ½T ¼ 4 0

1=Δx1

…

0

0

…

…

…

…

0

…

1=Δxn  1

1=Δxn  1

3 7 5;

ð20Þ

where Δxn ¼ xn þ 1  xn represents the distance between adjacent nodes in x-direction. Therefore, the vector fP f g had to be reorganized in the manner f…; P f ðxn Þ; P f ðxn þ 1 Þ; …g, where xn o xn þ 1 , to calculate the value of objection function (19). As can be seen the optimization problem depends on the mesh size p and can therefore be quite large. The defined objective function is non-linear, which makes the optimization problem much more difficult to solve. If m Z p the problem can also be seen as a non-linear least squares problem and for this reason two different methods have been used, BFGS and LM, which treat the problem in different ways. BFGS is a quasi-Newton method and treats the objective function in a classical manner. It finds a minimum of the objective function in a descendent way as fP f gk þ 1 ¼ fP f gk þ αfdgk -F R ðfP f gk þ 1 Þ o F R ðfP f gk Þ;

ð21Þ

where k represents the iterative step, α the step size and fdg the search direction. The search direction is calculated as fdgk ¼  ½Dk F ′R ðfP f gk Þ;

ð22Þ

where ½Dk represents an approximation of the inverse Hessian matrix and F ′R ðfP f gÞ represents the gradient of the objective function. In this manner the search direction calculation does not need equation solving and is therefore computationally less expensive. The approximation of the Hessian matrix ½Dk is calculated iteratively from the objective function gradients and initial approximation as ½Dk þ 1 ¼ ½Dk þ κ 1 fδgfδg  κ 2 ðfδgfvgtr þ fvgfδg Þ; tr

tr

ð23Þ

where fδg ¼ fP f gk þ 1  fP f gk represents an iteration step, vector fvg ¼ ½Dk fgg, where fgg represents the gradient change F ′R ðfP f gk þ 1 Þ  tr F ′R ðfP f gk Þ and coefficients κ 2 ¼ 1=ðfδg fggÞ and κ 1 ¼ κ 2 ð1 þ κ 2 fggtr fvgÞ [34]. As can be seen the second derivatives do not have to be evaluated, which is one of the advantages of this method and is therefore computationally less demanding. For initialization the identity matrix, ½D0 ¼ ½I, has been used. It can be shown that through the iteration process matrix ½Dk converges to the Hessian matrix F″ðfP f gÞ and for this reason quadratic convergence can be obtained when the solution is close to the local minimum. For this reason the number of iterative steps is small, which leads to faster computational time what was one of the reasons why to use it. After search direction establishment, the appropriate step size α needs to be defined for the purpose of reducing the number of iterative steps. The optimal step size in the search direction is one that minimizes the objective function the most. For this reason the exact line search with third-order polynomial approximation of F R ðαÞ ¼ F R ðfP f gk þ αfdgk Þ has been used, where fP f gk and fdgk are fixed [35]. The LM method treats the objective function in a different way, because it is designed to solve least-squares problems. The main concept of finding a better solution (21) is the same, only the calculation of search direction is different. Using the linear Taylor expansion of the residual vector f ðfP f g þ fdgÞ and regularized objective minimum condition F ′R ðfP f gÞ ¼ 0, the LM method incorporating first-order Tikhonov regularization is obtained as tr tr tr ð½Jtr k ½Jk þ μ½I þ λ½T ½TÞfdgk ¼  ½Jk f ðfP f gk Þ  λ½T ½TfP f gk ;

ð24Þ

where ½Jk ¼ ½JðfP f gk Þ is the Jacobian matrix, μ the damping parameter and [I] the identity matrix. As can be seen LM method is a hybrid between Gauss–Newton and Steepest descent controlled

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

by the damping parameter. The damping parameter also controls the step size fP f gk þ 1  fP f gk for which the step size has been chosen as one [36], α ¼ 1. If the step size α ¼ 1 does not minimize the objective function, F R ðfP f gk þ αfdgk Þ 4 FðfP f gk Þ, then the step size is reduced by αj þ 1 ¼ αj =3 and the calculation is repeated. Exact line search in this method has not been used. The damping parameter μ changes through the iteration process and is chosen to be large at the beginning for global convergence. The change is controlled with the gain ratio:

ρ¼

F R ðfP f gk Þ  F R ðfP f gk þ 1 Þ Zð0Þ  Zðαfdgk Þ

;

ð25Þ

where ZðÞ represents a linear Taylor expansion of the objective function F R ðfP f gÞ. A high gain ratio indicates good approximation with linear Taylor expansion and the damping parameter should be reduced in this case to improve the convergence speed. Otherwise, it represents a poor approximation and the damping parameter should be increased. The parameter is changed by the rule

μk þ 1 ¼ μk  maxf1=3; 1 ð2ρ  1Þ3 g:

ð26Þ

Owing to the damping parameter, the LM method is more stable than the classical Gauss-Newton method, which can produce unstable solution. The condition number of matrix ½Jtr ½J is reduced, which stabilizes the calculated search direction. The initial value of the damping parameter μ0 has been chosen to be μ0 ¼ 10  3  maxfaii g, where aii represents a diagonal element of matrix fJgtr 0 ½J0 [36]. The gradients necessary for both optimization methods have been calculated using first-order finite difference scheme, with pffiffiffiffiffi step size equal to machine accuracy, εR ; εR ¼ 10  16 . Three stopping criteria have been used for both optimization methods. The first one is the maximal number of iterative steps k 4kmax ¼ 50, the second one is the maximal gradient value J F′ðfP f gÞ J 1 o ε1 and the third one is minimal step size J fP f gk þ 1  fP f gk J r ε2 ð J fP f gk J þ ε2 Þ, where tolerances have been taken as ε1 ¼ 10  13 and ε2 ¼ 10  10 . If one of those criteria is fulfilled the optimization process is stopped. The regularization parameter λ has been determined according to the L-Curve method [33], which presents a novelty in this field. The main reason for choosing this heuristic technique is that there is no need of advanced knowledge of the noise, as for the discrepancy principle used by Trucu et al. [17]. The L-Curve method plots the logarithm of the norm of residual against the logarithm of the norm of the regularized derivative. Plotted function typically takes a characteristic L shape. This happens because the regularized derivative is a strictly decreasing function and the regularized residual is strictly increasing function of parameter λ. Therefore, the L-Curve method criterion for optimal λ selection is to pick the parameter corresponding to the corner of the curve which is defined as yL ðλÞ ¼ log ½ð½TfP f gÞtr ð½TfP f gÞ;

ð27Þ

xL ðλÞ ¼ log ½F R ðfP f gÞ  12 λð½TfP f gÞtr ð½TfP f gÞ:

ð28Þ

Due to the non-linearity of the problem, the L-Curve (27) and (28) has been evaluated numerically for a predefined set of regularized parameter values, λ ¼ f10  25 ; …; 10  5 g. The optimal value λ has been chosen according to the visual inspection of the plotted curve corner that is very distinct for the investigated inverse problem and for which the proposed approach does not fail. However, some precaution has to be made, because for some inverse problems the L-Curve method can fail, see Vogel [37].

79

4. Numerical examples The proposed numerical algorithm has been tested on some numerically generated test examples, where the base for them has been taken from Trucu et al. [17]. The efficiency of the BFGS and LM methods, as well as the solvability and the limitations of the inverse problem, under exact and noisy measurement data have been investigated. Optimization without regularization (λ ¼ 0) has been used only for exact measurement data that follows the numerical model exactly, just to analyse the solvability of the problem where there is no measurement and numerical error present. In real world both errors are present, numerical because of the numerical method and discretization and measurement because of the sensing error. Also to test the stability of the approach or solution independence, three different initial guesses for perfusion function have been used, close and far from the exact solution, where Trucu et al. [17] used only initial guess that is close. At the end, to study and analyse the effect of incorporating non-homogeneous tissue properties the testing has been done in two steps, where the first step represents tests examples for homogeneous tissue and the second one for piecewise homogeneous tissue. Two different test examples have been used, neglecting the metabolic heat generation I ¼ 0, on the computational domain: ðx; yÞ A ½0; 1  ½0; 0:1;

t A ½0; 1;

ð29Þ

where L ¼1 and H ¼0.1. The domain is thin in the y-direction. The initial dimensionless temperature has been chosen, u0 ðx; yÞ ¼ 0; Dirichlet boundary conditions: 8 for t ¼ 0; > <0 2 g 0 ðtÞ ¼ 0; g 1 ðtÞ ¼ e1=ðt  tÞ for t A ð0; 1Þ; ð30Þ > :0 for t Z 1; and the Neumann boundary conditions: h0 ðtÞ ¼ 0;

h2 ðtÞ ¼ 0;

ð31Þ

which all satisfy the conditions of Theorem 1. As can be seen from (29), the time domain has been chosen as finite one, despite the Theorem 1 and problem formulation. The reason for that lies in the Dirichlet boundary condition g 1 ðtÞ, because after time t ¼ 1 it is zero and the additional Neumann boundary value h1 ðtÞ is small enough or close to zero and do not change significantly after that. If the time domain is prolonged the inverse solution does not improve. Problem occurs only when the time is long enough that the residual vectors are numerically zero, which leads to singular matrices in the optimization routine and its failure. Two different perfusion functions, which we attempt to recover, have been used and are considered as exact ones: P f 1 ðxÞ ¼ 1 þ x2 ;

ð32Þ

P f 2 ðxÞ ¼ 1 þ 4x  4x2 :

ð33Þ

Indices 1 and 2 stand for first and second test example, respectively, where the first perfusion function is a smooth monotonically increasing function with low value of first-derivative and has been taken from Trucu et al. [17]. The second test function has been chosen to be a non-monotonic one with higher spacedependent first-derivative, which represents the new test example, to test the first-order regularization and L-Curve method, as well as to analyse the solution stability and limitations of the investigated inverse problem. Because the blood perfusion has been chosen to be dependent on x-direction only and due to the selected computational domain (29) and boundary conditions (30) and (31), the test examples can be thought as one-dimensional.

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

80

numerical or measurement error and therefore the “inverse crime” is deliberately committed. The reason for committing it is just to analyse the solvability of the inverse problem when no error is present. For 1% and 5% noise levels the data does no longer follow the numerical model (that is used for solving inverse problem) and therefore, no “inverse crime” is committed. The 1% and 5% noise can be also viewed as the total error (measurement and numerical), where for the used discretization the numerical error is small. The first part of the analysis represents test examples using homogeneous tissue as a ¼ 1 and k ¼ 1, and for initial guess of P f 0 ðxÞ ¼ 1 repeats a test example found in Trucu et al. [17] considering 1% noisy data. As for the second part, which is the main focus of this paper, to estimate the space-dependent perfusion inside non-homogeneous tissue the following piecewise homogeneous thermal property has been considered:

To test the global stability and solution independence, three different initial guesses for perfusion function have been chosen, P f 0 ðxÞ ¼ 0, P f 0 ðxÞ ¼ 1 and P f 0 ðxÞ ¼ 2, where first and last one represent initial guess far from the exact solution and the second one a close one. Since the analytical solution of the problem (2), (5)–(9), (29)–(33) could not be determined, the Neumann boundary value h1 ðtÞ has been generated numerically. As mentioned, the problem has been solved using exact as well as noisy measurement data, where white noise has been added by h1 ðt i Þ ¼ h1 ðt i Þð1 þ εm  βÞ; n

ð34Þ

8 > < 2; kðx; yÞ ¼ aðx; yÞ ¼ 0:5; > : 1;

2.4

2.4

2

2

1.6

1.6

1.2

0.8

0.4

0.4

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

2.4

2.4

2

2

1.6

1.6

1.2

0.8

0.4

0.4

0.2

0.4

0.6 x

0.6

0.8

1

0.6

0.8

1

1.2

0.8

0

ð35Þ

x

Pf

Pf

x

0

ðx; yÞ A ½0:4; 0:8  ½0; 0:1; ðx; yÞ A ½0:8; 1  ½0; 0:1:

1.2

0.8

0

ðx; yÞ A ½0; 0:4  ½0; 0:1;

With this set of numerical examples we gain a complete analysis of the proposed numerical algorithm, as well as the limitations of the inverse problem regarding the effects of non-homogeneous tissue, noisy data, initial guess and perfusion distribution.

Pf

Pf

where εm represents the percentage of noise and β A ½  1; 1 a random number. To simulate the measurement data the 0, 1% and 5% levels of noise have been considered, where the first one represents exact measurement data. To solve the inverse bio-heat problem numerically, the following time and space domain discretization have been used. The space domain ½0; 1  ½0; 0:1 has been discretized with N d ¼ 90  2 uniform cells, where 90 cells have been used in x-direction and 2 cells in y-direction, and the time domain [0,1] with Nt ¼ 100 uniform time intervals [17]. The same discretization has also been used to obtain the additional Neumann boundary value or measurement data h1 ðtÞ. Therefore, for exact data there is no

0.8

1

0

0

0.2

0.4 x

Fig. 1. Estimated Pf1 (a, b) and Pf2 (c, d) using different initial guesses and exact measurement data without regularization: (a, c) represent the solution using LM method and (b, d) the solution using BFGS method.

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

5. Results and discussion 5.1. Homogeneous tissue Results of numerical test examples for homogeneous tissue are presented first. Fig. 1 shows the obtained perfusion coefficient for exact Neumann data without regularization, using different initial guesses. It can be seen that the solution of inverse problem strongly depends on the initial guess and the method used, even for exact data. As it can be seen, if the initial guess is close to the exact perfusion function that we try to recover, the estimated perfusion will be in a good agreement with the exact perfusion, otherwise an oscillatory solution is obtained. The LM method produces a less oscillatory solution and describes perfusion function a little better than BFGS. The reason why both methods do not give better results for the initial guess of zero lies in the small sensitivity of perfusion coefficient near the boundary x ¼ 0. The gradient of the objective function is nearly 103 times smaller than the gradient of the objective function near measurement boundary. This can be seen from the plot of the Jacobian matrix elements for positions x ¼ 0:05 and x ¼ 0:95, which are presented in Fig. 2. It can be also seen that the time domain is suitably long, because the gradients after time t ¼ 1 becomes very small or close to zero. Thus, the oscillatory behaviour of the solution is a result of

81

averaging around the exact solution. As can be seen the inverse problem is ill-posed and the exact solution cannot be obtained even for exact measurement data. The convergence properties of both methods are shown in Fig. 3 for initial guess close to the exact solution, which yields better solution. For the same iteration step, the LM method obtains solution with lower values of the objective function, which means a better solution that can be also confirmed from Fig. 1. The BFGS method terminates prematurely and could not achieve such low value of objective function, because of a small step size or search direction, which did not improve the objective function. This can also be seen from the convergence tendency, as the objective function does not reduce for several iteration steps. The better LM method also starts with the steepest descent direction, due to a large parameter μ, which helps to stabilize the solution during the iteration process. The solution of inverse bio-heat problem strongly depends on the initial guess when no regularization is used, as can be seen from Fig. 1, and can be highly oscillatory if the initial guess is far from the exact solution, which indicates the ill-posedness of the problem. Although an oscillatory behaviour of the solution is obtained, the tendency of the exact functions is well described. To stabilize the oscillations, first-order Tikhonov regularization has been used. Therefore, Fig. 4 shows the estimated perfusion

0.0002

6E-07

0.00015

⎨J⎬5

⎨J⎬95

4E-07 0.0001

2E-07 5E-05

0 0

0.2

0.4

τ

0.6

0.8

0

1

0

0.2

0.4

τ

0.6

0.8

1

Fig. 2. Plot of Jacobian matrix elements value for position x ¼ 0:95 (a) and x ¼ 0:05 (b) at the final iteration.

10

-6

10

-8

10

-10

10

-12

10

-14

FR, ∇FR, μ

FR, ∇FR, μ

10-4

10

-4

10

-6

10

-8

10-10 10

-12

10

-14

10-16

10

-16

10-18

10

-18

10

-20

10

-20

0

10

20

30 k

40

50

0

10

20

30

40

50

k

Fig. 3. Convergence properties of BFGS and LM methods for both examples using exact data and initial guess P f 0 ðxÞ ¼ 1: (a) represents the first test case Pf1 and (b) the second test case Pf2.

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

2.4

2.4

2

2

1.6

1.6

Pf

Pf

82

1.2

1.2

0.8

0.8

0.4

0.4

0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

x

2.4

2.4

2

2

1.6

1.6

Pf

Pf

x

1.2

1.2

0.8

0.8

0.4

0.4

0

0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x

x

Fig. 4. Obtained perfusion coefficient for both test examples using regularization method on noise-free data: (a) and (c) represent the solution using LM method for first Pf1 and second Pf2 test function, respectively, and (b, d) the solution using BFGS method.

10

4

10

3

LM BFGS

LM BFGS

λ=1E-14

yL(λ)

yL(λ)

103

102

10

2

λ=1E-20

λ=1E-14 λ=1E-21 1

10

-23

10

-21

10

-19

10

-17

10

-15

10

-13

10

-11

10

-9

10

xL(λ)

-7

10 -21 10

10

-19

10

-17

10

-15

10

-13

10

-11

10

-9

xL(λ)

Fig. 5. Determination of the optimal regularization parameter λ using the L-Curve method for both test examples using first-order Tikhonov regularization, exact data and initial guess P f 0 ðxÞ ¼ 1: (a) represents the first test example Pf1 and (b) the second test example Pf2.

function using regularization process with noise-free data for both test examples and different initial guesses. The results are better and less oscillatory. Using the BFGS method provides a very good retrieval for the smooth monotonic function (32) and initial guess close to the exact solution. Otherwise, the oscillatory solution is obtained, which still depends on the initial guess. However, the LM

method incorporating the first-order regularization gave far better results, with very good agreement to the exact profile for monotonic, as well as non-monotonic test function and the solution does not depend on the initial guess as before. For a monotonic function and initial guess close to the exact solution the exact profile has been obtained. Some deviation from the exact solution

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

can be seen in the region near the boundary x ¼ 0, especially for non-monotonic function. For similar reasons as before, small gradients of the objective function or sensitivity coefficients and additional first-order regularization try to smooth out the unknown function in the first-order sense. For the perfusion function with small first-derivative in the region of small sensitivity, the function can be retrieved almost exactly as for the first test function, while for the second one the deviation occurs, because of the higher function first-derivative. The reason why the LM method is globally better and more stable than BFGS using regularization is hidden in the calculation of the search direction. Regularization for the LM method is included directly, through Eq. (24), which produces steps that are smooth in the first-order sense, while regularization for the BFGS method is incorporated indirectly through the gradient of the regularized objective function and penalizes large step sizes producing high first-order solution. For this reason the convergence of BFGS is very slow, while for the LM method stayed the same. The value of the optimal regularization parameter λ has been determined using the L-Curve method that is shown in Fig. 5, which depicts the obtained regularization parameter for noise-free data and initial guess close to the exact solution. As mentioned, the optimal value has been obtained by visual inspection of the plotted curve for a predefined set of regularized parameter values.

2.4

For both test cases and different initial guess the corner of the L-Curve function is distinctive and the parameter estimation is easy. From the figure it can also be seen that the LM method produces solutions with smaller values of the objective function, indicating better and less oscillatory solutions, than the BFGS method. The estimation of perfusion coefficient has been also done for 1% and 5% noisy measurement data. As seen from precedent analysis, the inverse bio-heat problem is ill-posed, for which the regularization process has to be used. While the LM method yields results independent of the initial guess, the BFGS method produces good results only when initial guess is close to the exact solution. Therefore, the results of inverse problem will be shown only for the initial guess P f 0 ðxÞ ¼ 1. Figs. 6 and 7 show perfusion coefficient for both test examples using 1% and 5% noisy data, respectively. The solution using BFGS method is acceptable only for the smooth monotonic function with small amount of noise. The reason for this lies in the first-order smoothness of the exact perfusion function. Also the solution using LM method is showing a good agreement with the first exact function. A small deviation occurs in the region x A ½0; 0:3 because of the previously stated reasons, which for noisy data becomes even more expressed. The perfusion coefficient in the region of small sensitivity, due to noisy data and first-order regularization,

2.4

exact LM λ=8E-10 BFGS λ=1E-9 Trucu [17]

2

2 1.6

Pf

Pf

1.6 1.2

1.2

0.8

0.8

0.4

0.4

0

83

0

0.2

0.4

0.6

x

0.8

0

1

exact LM λ=1E-10 BFGS λ=5E-11

0

0.2

0.4

x

0.6

0.8

1

Fig. 6. Estimation of perfusion function for both test examples using 1% noisy data and close initial guess P f 0 ðxÞ ¼ 1: (a) represents the first test function Pf1 and (b) the second test function Pf2.

2.4

exact LM λ=1E-8 BFGS λ=1E-8

2

2

1.6

1.6

Pf

Pf

2.4

1.2

1.2

0.8

0.8

0.4

0.4

0

0

0.2

0.4

x

0.6

0.8

1

0 0

exact LM λ=5E-10 BFGS λ=1E-9

0.2

0.4

x

0.6

0.8

1

Fig. 7. Estimation of perfusion function for both test examples using 5% noisy data and close initial guess P f 0 ðxÞ ¼ 1: (a) represents the solution of first test function Pf1 and (b) the solution of second test example.

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

84

cannot be retrieved anymore and therefore information is lost, as can be seen distinctly for the non-monotonic perfusion function. However, for a small amount of noise, the perfusion distribution is still in good agreement with the exact one for the region near the boundary measurement x A ½0:4; 1. Higher noise level makes perfusion estimation even more difficult, especially for the nonmonotonic function, as can be seen from Fig. 7. The deviation is more noticeable for the non-monotonic function with high value of first-derivative, while for a smooth monotonic function the reconstruction is still possible. However, the LM method proves to be more superior, which produces solutions that are more in agreement with the exact solution and independent of initial guess. Also the convergence is a lot faster, as presented in Fig. 8 for 1% noisy data, finishing within 5 iteration steps. The BFGS produces very slow convergence and also poorer solutions. The convergence properties of both methods are the same in 5% noisy data and for this reason are not presented. The regularization parameter has been obtained in the same way as for noise-free data and is presented in Fig. 9 for 1% noisy data, which is similar for 5% noisy data. As can be seen the corner of the L-Curve function is still distinctive, which makes the visual determination of the optimal regularization parameter still possible. From Fig. 9 we can also notice that the solutions for first

test examples using both methods are very similar in the sense of objective function value and perfusion function first-derivative as can be also seen from Fig. 6. However, for BFGS method we have to take the initial guess that is close to the exact solution which is hard to do when the exact solution is not known. For the first test example and 1% noisy data the solution has also been compared with the solution of Trucu et al. [17] and is shown in Fig. 6. The BFGS provides the best estimate solution for test function (32) followed by LM and by the NAG routine of Trucu et al. [17]. 5.2. Piecewise homogeneous tissue The second part treats perfusion function estimation in piecewise homogeneous tissue, composed of three different homogeneous layers described in (35), under the same boundary conditions and test examples. The previous analysis for homogeneous tissue shows that due to ill-posedness of the problem the regularization process has to be used even if exact measurement data is used, and that BFGS method produces good results only when initial guess is close to the exact solution. Therefore, the solution for piecewise homogeneous tissue will be represented only for the initial guess, P f 0 ðxÞ ¼ 1, using first-order Tikhonov regularization.

10-4

10

-5

10

-5

10

-6

10

-6

10

-7

10

-7

10

-8

10

-8

10

-9

10

-9

10

-10

10

-11

FR, ∇FR, μ

FR, ∇FR, μ

10-4

0

10

20

30

40

50

10

-10

10

-11

0

10

20

k

30

40

50

k

10

5

10

4

10

3

10

2

10

1

10

0

10

-1

10

-2

10

-3

106

LM BFGS

10

LM BFGS

4

λ=5E-11

λ=1E-9 λ=8E-10

yL(λ)

yL(λ)

Fig. 8. Convergence of optimization methods for both test examples, close initial guess P f 0 ðxÞ ¼ 1 and 1% noisy data: (a) represents the first test example and (b) the second test example.

102 λ=1E-10 10

10

10

-7

10 xL(λ)

-6

0

-2

10

-6

10

-5

xL(λ)

Fig. 9. L-Curve method for determination of optimal λ for both test examples using close initial guess P f 0 ðxÞ ¼ 1 and 1% noisy data: (a) represents the first test example and (b) the second test example.

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

Fig. 10 shows the estimated perfusion function for exact measurement data using regularization process. The conclusions are the same as for the homogeneous tissue. The BFGS method produces good solution only for the smooth monotonic function,

while LM method produces good solution for both test examples regardless of the initial guess. The solution using BFGS method is now showing even more oscillatory behaviour than for homogeneous tissue. The reason lies in the second low conductive layer,

2.4

exact LM λ=5E-20 BFGS λ=1E-14

2

2

1.6

1.6

Pf

Pf

2.4

1.2

1.2

0.8

0.8

0.4

0.4

0

0

0.2

0.4

0.6

85

0.8

0

1

exact LM λ=5E-20 BFGS λ=1E-14

0

0.2

0.4

x

0.6

0.8

1

x

Fig. 10. Estimated perfusion function for both test examples using regularization process, exact measurement data and close initial guess P f 0 ðxÞ ¼ 1: (a) represents the first test and (b) the second test example.

2.4

exact LM λ=2E-10 BFGS λ=1E-9

2

2

1.6

1.6

Pf

Pf

2.4

1.2

1.2

0.8

0.8

0.4

0.4

0

0

0.2

0.4

0.6

0.8

0

1

exact LM λ=2E-10 BFGS λ=5E-10

0

0.2

0.4

x

0.6

0.8

1

x

Fig. 11. Obtained perfusion function for both test cases using 1% noisy data and close initial guess: (a) represents the first test example Pf1 and (b) the second test example Pf2.

2.4

exact LM λ=1E-8 BFGS λ=1E-8

2

2

1.6

1.6

Pf

Pf

2.4

1.2

1.2

0.8

0.8

0.4

0.4

0

0

0.2

0.4

0.6

x

0.8

1

0

exact LM λ=1E-10 BFGS λ=1E-8

0

0.2

0.4

0.6

0.8

1

x

Fig. 12. Estimated perfusion function for both test cases using 5% noisy data and close initial guess: (a) represents the solution of first test example Pf1 and (b) the solution of second test example Pf2.

86

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87

which blocks the penetration of the heat wave through the domain. Therefore, the information of the blood perfusion inside of the domain cannot be gathered, which reflects in even smaller sensitivity of perfusion for the region x A ½0; 0:4. This is especially visible for the second test example. On the contrary the LM method is showing less sensitivity to the piecewise homogeneous tissue and for this reason describes the perfusion function better. Consequently the piecewise homogeneous tissue makes the estimation of perfusion function more difficult especially after the low conducting tissue. This is clearly visible for the second nonmonotonic function in the region of x A ½0; 0:3. Convergence of both methods for piecewise homogeneous tissue is similar to the convergence for homogeneous tissue and is therefore not presented here, as well as the L-Curve method or determination of optimal regularization parameter. We also tested the reconstruction of perfusion function under noisy measurement data. Fig. 11 shows obtained perfusion function for both test examples using 1% noise level. The conclusion is similar as for the homogeneous tissue. The information about perfusion is completely lost in the region of small sensitivity x A ½0; 0:4 for both test functions. Seemingly the BFGS method provides a better solution for first test function in this region, however this is due to the close initial guess, P f 0 ðxÞ ¼ 1, which is not changed during the optimization process. This is especially visible for the second test function and therefore the information is considered to be lost for BFGS method as well. For this reason the LM method still shows superior performance, because it is independent of the initial guess and can retrieve the perfusion function better, particularly in the region near boundary measurement x A ½0:5; 1. If the noise level is large, the region of information loss is also increased, which can be also seen in Fig. 12 for 5% level of noise. Therefore, the successful reconstruction can be made only in the region of boundary measurement, x A ½0:5; 1, for reasonably smooth monotonic functions, while for non-monotonic function the information can be lost due to the high first-order regularization. For noisy measurement data the convergence properties of both optimization methods is similar to the convergence of homogeneous tissue already described above (Fig. 8) and are therefore not presented again. The LM method is still faster from BFGS method and converges in 5 iteration steps. As well the L-Curve function for noisy data is omitted here, because it is very similar to the inverse problem of homogeneous tissue and the optimal values of regularization parameter are given in Figs. 11 and 12.

6. Conclusion The inverse bio-heat problem of space-dependent perfusion coefficient estimation inside homogeneous, as well as piecewise homogeneous tissue, has been solved. The approach is based on the Pennes bio-heat equation and has been solved numerically and investigated in the presence of known initial, Dirichlet and Neumann boundary conditions, together with additional exact and noisy Neumann boundary measurement. To solve the inverse problem, a BEM based optimization algorithm has been developed. The first part represents a direct solver based on BEM using mixed interpolation functions, to solve direct bio-heat problems with assumed perfusion, coupled with the second step representing optimization routines. Two different optimization methods have been used, BFGS quasi-Newton and Levenberg-Marquardt, for minimizing the firstorder Tikhonov regularization. Therefore, the unknown perfusion coefficient function has been found as a global minimizer of the prescribed objective function. The methodology has been tested on different test examples using exact and noisy data, where the efficiency of the optimization methods, effect of different initial

guess, as well as the effect of non-homogeneous properties and perfusion distribution have been studied. Results show that the solution depends on the initial guess for the case of exact data using no regularization and good solutions can be obtained only for initial guess close to the exact one, which also shows that the inverse problem is ill-posed. First-order regularization helps to stabilize the solution and makes solutions using the LM method, independent of initial guess for exact, as well as noisy data, with very good agreement obtained for smooth monotonic perfusion function. In the case of non-monotonic distribution, the information in the region of small sensitivity coefficients is lost, especially in the presence of noisy data. Therefore, the unknown function can be efficiently retrieved only in the region near boundary measurement, with low noise levels. For higher noise level, the reconstruction is possible only for smooth monotonic functions in first-order sense. Problem difficulties of estimating perfusion coefficient only increase when considering non-homogeneous properties of the tissue, especially if a layer of low conducting tissue occurs, because it blocks the heat wave propagation into the domain and therefore the information is not retrieved. Therefore, the perfusion can be successfully retrieved only in the region near the boundary measurement or before the low conducting tissue and at low level of noise. The optimal regularization parameter has been determined using the L-Curve method, which is distinctive for exact as well as noisy data and makes the determination by visual inspection possible. According to the obtained results, the LM method is shown to be more superior and is more appropriate for treating problems under exact and noisy measurement data. It is stable, less sensitive to piecewise homogeneous tissue and has fast convergence. Therefore, the proposed BEM based algorithm together with LM first-order Tikhonov regularization and L-Curve method produces an efficient algorithm for perfusion estimation in homogeneous and piecewise homogeneous tissue. Further work will be focused on including higher-order regularization and improved L-Curve method, trying to retrieve the perfusion coefficient in regions of small sensitivity, as well as upgrading the problem to multidimensions.

References [1] Liu J, Xu LX. Boundary information based diagnostics on the thermal states of biological bodies. Int J Heat Mass Transf 2000;43(16):2827–39. [2] Partridge P, Wrobel L. A coupled dual reciprocity BEM/genetic algorithm for identification of blood perfusion parameters. Int J Numer Methods Heat Fluid Flow 2009;19(1):25–38. [3] Thiebaut C, Lemonnier D. Three-dimensional modelling and optimisation of thermal fields induced in a human body during hyperthermia. Int J Therm Sci 2002;41(6):500–8. [4] Zhao G, Zhang H, Guo X, Luo D, Gao D. Effect of blood flow and metabolism on multidimensional heat transfer during cryosurgery. Med Eng Phys 2007;29 (2):205–15. [5] Xu LX, Anderson GT. Techniques for measuring blood flow in the microvascular circulation. In: Biomechanical systems techniques and applications, Vol. 4: Biofluid methods in vascular and pulmonary systems. CRC Press; 2000 (Chap. 4). [6] Liu J, Zhou Y, Deng Z. Sinusoidal heating method to noninvasively measure tissue perfusion. IEEE Trans Biomed Eng 2002;49(8):867–77. [7] Yue K, Zhang X, Zuo Y. Noninvasive method for simultaneously measuring the thermophysical properties and blood perfusion in cylindrically shaped living tissues. Cell Biochem Biophys 2008;50(1):41–51. [8] Ricketts P, Mudaliar A, Ellis B, Pullins C, Meyers L, Lanz O, et al. Non-invasive blood perfusion measurements using a combined temperature and heat flux surface probe. Int J Heat Mass Transf 2008;51(23):5740–8. [9] Loulou T, Scott E. An inverse heat conduction problem with heat flux measurements. Int J Numer Methods Eng 2006;67(11):1587–616. [10] Sandberg M, Zhang Q, Styf J, Gerdle B, Lindberg L-G. Non-invasive monitoring of muscle blood perfusion by photoplethysmography: evaluation of a new application. Acta Physiol Scand 2005;183(4):335–43. [11] Liu J, Xu LX. Estimation of blood perfusion using phase shift in temperature response to sinusoidal heating at the skin surface. IEEE Trans Biomed Eng 1999;46(9):1037–43. [12] Deng Z, Liu J. Parametric studies on the phase shift method to measure the blood perfusion of biological bodies. Med Eng Phys 2000;22(10):693–702.

J. Iljaž, L. Škerget / Engineering Analysis with Boundary Elements 39 (2014) 75–87 [13] Upreti S, Jeje A. A noninvasive technique to determine peripheral blood flow and heat generation in a human limb. Chem Eng Sci 2004;59(21):4415–23. [14] Erhart K, Divo E, Kassab A. An evolutionary-based inverse approach for the identification of non-linear heat generation rates in living tissues using a localized meshless method. Int J Numer Methods Heat Fluid Flow 2008;18 (3/4):401–14. [15] Trucu D, Ingham D, Lesnic D. Inverse temperature-dependent perfusion coefficient reconstruction. Int J Non-Linear Mech 2010;45:542–9. [16] Trucu D, Ingham D, Lesnic D. Inverse space-dependent perfusion coefficient identification. J Phys Conf Ser 2008;135:012098 (8 pp). [17] Trucu D, Ingham D, Lesnic D. Space-dependent perfusion coefficient identification in the transient bio-heat equation. J Eng Math 2010;67:307–15. [18] Trucu D, Ingham D, Lesnic D. Inverse time-dependent perfusion coefficient identification. J Phys Conf Ser 2008;124:012050 (28 pp). [19] Trucu D, Ingham D, Lesnic D. An inverse coefficient identification problem for the bio-heat equation. Inverse Probl Sci Eng 2009;17(1):65–83. [20] Rodrigues F, Orlande H, Dulikravich G. Simultaneous estimation of spatially dependent diffusion coefficient and source term in a nonlinear 1D diffusion problem. Math Comput Simul 2004;66:409–24. [21] Yang L, Yu J-N, Deng Z-C. An inverse problem of identifying the coefficient of parabolic equation. Appl Math Model 2008;32(10):1984–95. [22] Pennes H. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol 1948;1(2):93–122. [23] Arkin H, Xu LX, Holmes KR. Recent developments in modeling heat transfer in blood perfused tissues. IEEE Trans Biomed Eng 1994;41(2):97–107. [24] Khaled A-RA, Vafai K. The role of porous media in modeling flow and heat transfer in biological tissues. Int J Heat Mass Transf 2003;46(26):4989–5003. [25] Nakayama A, Kuwahara F. A general bioheat transfer model based on the theory of porous media. Int J Heat Mass Transf 2008;51(11):3190–9.

87

[26] Ren Z, Liu J, Wang C, Jiang P. Boundary element method (BEM) for solving normal or inverse bio-heat transfer problem of biological bodies with complex shape. J Therm Sci 1995;4:117–24. [27] Tadi M, Klibanov M, Cai W. An inverse method for parabolic equations based on quasireversibility. Comput Math Appl 2002;43:927–41. [28] Ramm A. An inverse problem for the heat equation. J Math Anal Appl 2001;264(2):691–7. [29] Denisov A. Elements of the theory of inverse problems. Utrecht, Netherlands: VSP; 1999. [30] Popov V, Power H, Škerget L. Domain decomposition techniques for boundary elements. Application to fluid flow. Southampton, Boston: WIT Press; 2007. [31] Ramšak M, Škerget L, Hriberšek M, Žunič Z. A multidomain boundary element method for unsteady laminar flow using stream functionvorticity equations. Eng Anal Bound Elem 2005;29(1):1–14. [32] Ramšak M, Škerget L. Iterative least squares methods for solving overdetermined matrix for mixed boundary elements. Z Angew Math Mech Berl 2000;80(Suppl3):S657–8. [33] Aster R, Borchers B, Thurber C. Parameter estimation and inverse problems. London: Elsevier; 2012. [34] Frandsen P, Jonasson K, Nielsen H, Tingleff O. Unconstrained optimization, informatics and mathematical modelling. Technical University of Denmark, DTU. [35] Vanderplaats G. Numerical optimization techniques for engineering design: with applications. New York: McGraw-Hill; 1984. [36] Madsen K, Nielsen H, Tingleff O. Methods for non-linear least squares problems. Informatics and Mathematical Modelling, Technical University of Denmark, DTU. [37] Vogel C. Non-convergence of the L-curve regularization parameter selection method. Inverse Probl 1996;12:535–47.