# Image Reconstruction of Two-Dimensional Highly Scattering Inhomogeneous Medium Using MAP-Based Estimation

Authors: Hong Qi; Yaobin Qiao; Shuangcheng Sun; Yuchen Yao; Liming Ruan
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A maximum a posteriori (MAP) estimation based on Bayesian framework is applied to image reconstruction of two-dimensional highly scattering inhomogeneous medium. The finite difference method (FDM) and conjugate gradient (CG) algorithm serve as the forward and inverse solving models, respectively. The generalized Gaussian Markov random field model (GGMRF) is treated as the regularization, and finally the influence of the measurement errors and initial distributions is investigated. Through the test cases, the MAP estimate algorithm is demonstrated to greatly improve the reconstruction results of the optical coefficients.


---

## Body

## 1. Introduction

The optical information reconstruction in nonhomogeneous medium exposed to collimated short-pulse irradiation has been a research focus in transient radiative transfer recently. It is widely used in many fields, such as nondestructive testing, optical tomography, infrared remote sensing, information processing, combustion diagnosis, biology, and medicine [1–7]. Among these fields, the diffuse optical tomography (DOT) is an emerging technology that can reconstruct the optical properties of internal biological tissues from the measured transient transmittance and reflectance signals at the boundary of tissues. The DOT is characterized as being noninvasive, portable, and producing real-time images of clinically relevant parameters. Commonly, the near-infrared laser pulse is used as the light source in DOT, which does no harm to the detected biological issues. The photon propagation in biological medium which is highly scattering always obeys diffusion equation. It offers the significant advantages of both flexibility to geometry and heterogeneity of tissue and is easy to be solved using numerical methods [8].Using the information carried by scattering photon, combined with the photon transmission model and the optimization algorithm, the information reconstruction can be converted to an inverse problem. The difficulties are how to choose radiative transfer models, detect transmission and reflected signals, and apply inverse optimization algorithms. Many research groups have done a lot of work in this field. Yamada and Hasegawa applied the finite element method (FEM) to solve the parabolic diffusion approximation equation and the transient transmission of pulse in 2D and 3D cylindrical media [9]. Schweiger et al. used the FEM to solve the diffusion approximation equation and the problems of light spread in biological tissue under different boundary conditions and light source [10, 11]. Guven et al. provided a computationally efficient iterative algorithm to simultaneously estimate the optical image and the unknown a priori model parameters [12]. More recently, Wang and Zuo proposed a variational model based on the maximum posterior probability (MAP) to reconstruct images degraded by Rician noise [13]. Shi et al. developed an efficient L1 regularization-based reconstruction algorithm to increase the computational speed with low memory consumption [14]. Zhang et al. proposed a new method to resolve the fluorescence molecular tomography inverse problem based on MAP estimation with structural priors in a Bayesian framework [15].However, to the authors’ best knowledge, there are few reports concerning the image reconstruction of biological tissues with the MAP estimation in the field of transient radiative transfer problems [12, 16]. In this work, based on the time-domain diffusion equation, the diffusion coefficient distribution of the biological tissues is reconstructed under the condition of uniform absorption coefficient firstly, and then the absorption coefficient and reduced scattering coefficient distribution are reconstructed simultaneously considering inhomogeneous absorption coefficient. The edges of the reconstruction image become clearer with the MAP estimation. Compared with the unconstrained method, the MAP estimation algorithm greatly improves the reconstruction quality and reduces the normalized root-mean square error (NRMSE).


## 2. Materials and Methods

### 2.1. The Principle of the Forward Problem

To reconstruct the optical parameters of semitransparent media, a forward model is needed to simulate the photon propagation in semitransparent media and calculate transmission and reflected signals. The photon propagation in highly scattering medium obeys diffusion equation [12, 16]. The two-dimensional time-resolved diffusion equation is given as [10](1)∂G(r,t)∂t=∂∂xD(r)∂G(r,t)∂x+∂∂yD(r)∂G(r,t)∂y-cnμaG(r,t)+S(r,t),∀r∈Ω,where r denotes the location and t denotes the time. Ω indicates the medium. n is the refractive index. c is the speed of light in vacuum. G(r,t) is the incident radiation. S(r,t) is the radiative source term. μa is the absorption coefficient. D(r) is the diffusion coefficient, which is determined by (2)D(r)=c/n3μa(r)+μs′(r),where μs′=(1-g)μs is called reduced scattering coefficient and μs is the scattering coefficient. g is the scattering asymmetry parameter.The initial and boundary conditions are expressed as follows.Initial condition(3)Gr,0=0,∀r∈Ω∪∂Ω.Considering that the refractive indexes in the two sides of boundary are inconsistent, we employ the partial current boundary condition [17] in the four boundaries of two-dimensional media:(4)cnGr,t+2ADrn·∇Gr,t=0,∀r∈∂Ω,where ∂Ω indicates the boundary. n is the outward normal unit vector of the border. A can be obtained by(5)A=2/(1-R0)-1+cos⁡θc31-cos⁡θc2,where R0=n-12/n+12 and the critical angle θc=arcsin1/n.By replacing the spatial and temporal derivatives with the fully implicit scheme of finite difference approximations, (1) can be expressed as follows:(6)Gi,jk+1-Gi,jkΔt=Di+1,j+Di,jGi+1,jk+1-Gi,jk+12Δx2-Di-1,j+Di,jGi,jk+1-Gi-1,jk+12Δx2+Di,j+1+Di,jGi,j+1k+1-Gi,jk+12Δy2-Di,j-1+Di,jGi,jk+1-Gi,j-1k+12Δy2-cnμai,jGi,jk+1+12Si,jk+Si,jk+1.For square grids (Δx=Δy=Δ), the diffusion equation can be transformed into linear equations:(7)-(Ai+1,j+Ai,j)Gi+1,jk+1-(Ai-1,j+Ai,j)Gi-1,jk+1-(Ai,j+1+Ai,j)Gi,j+1k+1-(Ai,j-1+Ai,j)Gi,j-1k+1+cΔtn1+4Ai,j+Ai+1,j+Ai-1,j+Ai,j+1+Ai,j-1+cΔtnμai,jGi,jk+1=Gi,jk+Δt2Si,jk+Si,jk+1.Equation (7) can be abbreviated as(8)KGn+1=Gn+S-n+1/2.Equation (8) is solved by a stabilized biconjugate gradient method (Bi-CGSTAB) [18] in the present study.


### 2.2. The Principle of the Inverse Problem

The reconstruction of unknown optical parameters is an ill-posed inverse problem. So the regularization technique should be used to make the solution well-behaved. For maximum a posteriori (MAP) estimation deduced from Bayesian theory can lead to Tikhonov regularization, the inversion of the distribution of the optical parameters is converted to the Bayesian estimation problems for image reconstruction by means of probability distribution models [15, 19]. The MAP estimation of the unknown image based on Bayesian theorem can be written as [20](9)α^MAP=argmax⁡αlog⁡P(M∣α)+log⁡P(α),where α denotes the retrieval optical parameters. M denotes the measured value. P(M∣α) is the conditional probability density and P(α) is the prior probability density of the image. The first term on the right is the agreement term and the second is the regularization term. In the inverse transient radiative problem, the objective function is usually given by(10)Fα=12∑s∑d∑nGs,dnα-Ms,dnMs,dn2+ψα,where α denotes the retrieval optical parameters. Gs,dn(α) is the estimated value of the dth detector at the nth moment when the sth light source is working. Ms,dn is the measured value of the dth detector at the nth moment when the sth light source is working. ψ(α) is the regularization term.To preserve the edges of the image, the generalized Gaussian Markov random field (GGMRF) model is used for image reconstruction. The regularization term can be defined as [21](11)ψ(α)=-1pσp∑{s,r}∈Nbs-rαs-αrp,where p is the image sharpness parameter. σ is the scale parameter. N is the set of all neighboring pixel pairs. r is the neighboring position of s. bs-r is the weight coefficient.The reconstruction of optical image is equal to minimization of the objective function, which needs to obtain the gradient of objective function with respect to the optical parameters. The adjoint differentiation (AD) method is an efficient method to calculate the gradient, which can reduce the consumption of computing time. The gradient of objective function with respect to the optical parameters can be expressed as(12)dFdαr=∑n∈T∑p∈ΩdFdGpn∂Gpn∂αr,r∈Ω,where ∂Gpn/∂αr is the sensitivity of Gpn with respect to α and can be obtained by differentiating (8) with respect to α, which yields(13)dKdαrGn+1+K∂Gn+1∂αr=0.The intensity G and matrix K can be calculated by the forward model. The solution of (18) is similar to the solution forward model using the Bi-CGSTAB. dF/dGpn is the derivative of objective function F with respect to intensity G and can be obtained recursively using the chain rule(14)dFdGpn=∑q∈ΩdFdGqn+1dGqn+1dGpn+∂F∂Gpnwith(15)dFdGPN=∂F∂GPN,N=max⁡n∈tsam,where tsam is the sampling time which means the time consumed in the measurement of transmission and reflected signals. ∂F/∂Gpn is obtained by partially differentiating objective function F with respect to Gpn:(16)∂F∂Gpn=Mpn-GpnMpn2,p∈∂Ω0,otherwise.Differentiating (8) with respect to Gpn, we can obtain(17)dGqn+1dGpn=K-1Ip,p,q∈Ωwith(18)Ip=1,q=p0,q≠p.The gradient of the objective function with respect to the optical parameters can be calculated by (12)–(18). After calculation of the gradient, CG method can be employed to solve inverse problem [12].The flowchart of the image reconstruction is shown in Figure1.Figure 1
The flowchart of the image reconstruction.


## 2.1. The Principle of the Forward Problem

To reconstruct the optical parameters of semitransparent media, a forward model is needed to simulate the photon propagation in semitransparent media and calculate transmission and reflected signals. The photon propagation in highly scattering medium obeys diffusion equation [12, 16]. The two-dimensional time-resolved diffusion equation is given as [10](1)∂G(r,t)∂t=∂∂xD(r)∂G(r,t)∂x+∂∂yD(r)∂G(r,t)∂y-cnμaG(r,t)+S(r,t),∀r∈Ω,where r denotes the location and t denotes the time. Ω indicates the medium. n is the refractive index. c is the speed of light in vacuum. G(r,t) is the incident radiation. S(r,t) is the radiative source term. μa is the absorption coefficient. D(r) is the diffusion coefficient, which is determined by (2)D(r)=c/n3μa(r)+μs′(r),where μs′=(1-g)μs is called reduced scattering coefficient and μs is the scattering coefficient. g is the scattering asymmetry parameter.The initial and boundary conditions are expressed as follows.Initial condition(3)Gr,0=0,∀r∈Ω∪∂Ω.Considering that the refractive indexes in the two sides of boundary are inconsistent, we employ the partial current boundary condition [17] in the four boundaries of two-dimensional media:(4)cnGr,t+2ADrn·∇Gr,t=0,∀r∈∂Ω,where ∂Ω indicates the boundary. n is the outward normal unit vector of the border. A can be obtained by(5)A=2/(1-R0)-1+cos⁡θc31-cos⁡θc2,where R0=n-12/n+12 and the critical angle θc=arcsin1/n.By replacing the spatial and temporal derivatives with the fully implicit scheme of finite difference approximations, (1) can be expressed as follows:(6)Gi,jk+1-Gi,jkΔt=Di+1,j+Di,jGi+1,jk+1-Gi,jk+12Δx2-Di-1,j+Di,jGi,jk+1-Gi-1,jk+12Δx2+Di,j+1+Di,jGi,j+1k+1-Gi,jk+12Δy2-Di,j-1+Di,jGi,jk+1-Gi,j-1k+12Δy2-cnμai,jGi,jk+1+12Si,jk+Si,jk+1.For square grids (Δx=Δy=Δ), the diffusion equation can be transformed into linear equations:(7)-(Ai+1,j+Ai,j)Gi+1,jk+1-(Ai-1,j+Ai,j)Gi-1,jk+1-(Ai,j+1+Ai,j)Gi,j+1k+1-(Ai,j-1+Ai,j)Gi,j-1k+1+cΔtn1+4Ai,j+Ai+1,j+Ai-1,j+Ai,j+1+Ai,j-1+cΔtnμai,jGi,jk+1=Gi,jk+Δt2Si,jk+Si,jk+1.Equation (7) can be abbreviated as(8)KGn+1=Gn+S-n+1/2.Equation (8) is solved by a stabilized biconjugate gradient method (Bi-CGSTAB) [18] in the present study.


## 2.2. The Principle of the Inverse Problem

The reconstruction of unknown optical parameters is an ill-posed inverse problem. So the regularization technique should be used to make the solution well-behaved. For maximum a posteriori (MAP) estimation deduced from Bayesian theory can lead to Tikhonov regularization, the inversion of the distribution of the optical parameters is converted to the Bayesian estimation problems for image reconstruction by means of probability distribution models [15, 19]. The MAP estimation of the unknown image based on Bayesian theorem can be written as [20](9)α^MAP=argmax⁡αlog⁡P(M∣α)+log⁡P(α),where α denotes the retrieval optical parameters. M denotes the measured value. P(M∣α) is the conditional probability density and P(α) is the prior probability density of the image. The first term on the right is the agreement term and the second is the regularization term. In the inverse transient radiative problem, the objective function is usually given by(10)Fα=12∑s∑d∑nGs,dnα-Ms,dnMs,dn2+ψα,where α denotes the retrieval optical parameters. Gs,dn(α) is the estimated value of the dth detector at the nth moment when the sth light source is working. Ms,dn is the measured value of the dth detector at the nth moment when the sth light source is working. ψ(α) is the regularization term.To preserve the edges of the image, the generalized Gaussian Markov random field (GGMRF) model is used for image reconstruction. The regularization term can be defined as [21](11)ψ(α)=-1pσp∑{s,r}∈Nbs-rαs-αrp,where p is the image sharpness parameter. σ is the scale parameter. N is the set of all neighboring pixel pairs. r is the neighboring position of s. bs-r is the weight coefficient.The reconstruction of optical image is equal to minimization of the objective function, which needs to obtain the gradient of objective function with respect to the optical parameters. The adjoint differentiation (AD) method is an efficient method to calculate the gradient, which can reduce the consumption of computing time. The gradient of objective function with respect to the optical parameters can be expressed as(12)dFdαr=∑n∈T∑p∈ΩdFdGpn∂Gpn∂αr,r∈Ω,where ∂Gpn/∂αr is the sensitivity of Gpn with respect to α and can be obtained by differentiating (8) with respect to α, which yields(13)dKdαrGn+1+K∂Gn+1∂αr=0.The intensity G and matrix K can be calculated by the forward model. The solution of (18) is similar to the solution forward model using the Bi-CGSTAB. dF/dGpn is the derivative of objective function F with respect to intensity G and can be obtained recursively using the chain rule(14)dFdGpn=∑q∈ΩdFdGqn+1dGqn+1dGpn+∂F∂Gpnwith(15)dFdGPN=∂F∂GPN,N=max⁡n∈tsam,where tsam is the sampling time which means the time consumed in the measurement of transmission and reflected signals. ∂F/∂Gpn is obtained by partially differentiating objective function F with respect to Gpn:(16)∂F∂Gpn=Mpn-GpnMpn2,p∈∂Ω0,otherwise.Differentiating (8) with respect to Gpn, we can obtain(17)dGqn+1dGpn=K-1Ip,p,q∈Ωwith(18)Ip=1,q=p0,q≠p.The gradient of the objective function with respect to the optical parameters can be calculated by (12)–(18). After calculation of the gradient, CG method can be employed to solve inverse problem [12].The flowchart of the image reconstruction is shown in Figure1.Figure 1
The flowchart of the image reconstruction.


## 3. Results and Discussions

The computational model which contains three inhomogeneous media is shown in Figure2. The diffusion coefficient in the whole region Ω is unknown. There are 16 optical fibers distributed uniformly over the boundary. When a short pulsed light source irradiates from one of the fibers, the other 15 fibers receive time-resolved signal at the same time. Move the source around the medium to obtain more detector readings. The refractive index is fixed as n=1.4, sampling time tsm=0.5 ns, and residual ε=0.5×10-4. The forward calculation model is 2×2 cm2 with 21×21 grids and 50 time steps.Figure 2
The schematic of the model with fibers.To compare the overall accuracy of the reconstructions, the NRMSE is introduced and defined as [22](19)NRMSE=∑rNα^r-α~r2∑rNα~r21/2,where α^r is the retrieval value of the optical parameter at the rth pixel. α~r is the true value of the optical parameter at the rth pixel. N is the total number of the pixels in the medium.


### 3.1. Image Reconstruction of Diffusion Coefficient

Assume that the absorption coefficient is constant asμa=0.05 cm−1, the diffusion coefficient of background is D=1.5 cm2 ns−1, and A, B, and C are heterogeneous bodies whose diffusion coefficients are 3.0 cm2 ns−1, 3.0 cm2 ns−1, and 1.0 cm2 ns−1, respectively. The absorption coefficient of inclusions is the same to background. Initial value is uniform and the diffusion coefficient is D0=1.0 cm2 ns−1. The real distribution and reconstruction results for nonhomogeneous medium are shown in Figure 3.Figure 3
The distribution of diffusion coefficient in nonhomogeneous medium: (a) the true distribution; (b) the reconstructed result (NRMSE = 0.17878).
(a)
(b)Figure3(b) shows that the locations and distributions of diffusion coefficient of three heterogeneous bodies are well reconstructed. The conjugate gradient method is used in the study of inverse problems. Since calculating the gradient of diffusion coefficient is needed, there are some step changes appearing in the retrieved results.Although the reconstruction peak of the diffusion coefficient is close to the real value, the shapes of three inclusions have been blurred in the computing process. So the MAP method is introduced to reconstruct the edge better. The variation of reconstructed images and NRMSE with number of iterations of the optimization algorithm with MAP is shown in Figure4. The shape becomes more accurate and the value of NRMSE becomes smaller with the increasing of iterations.Figure 4
The reconstructed diffusion coefficients with MAP for different number of iterations (k). (a) k=10, NRMSE = 0.22287, (b) k=50, NRMSE = 0.18365, (c) k=100, NRMSE = 0.15461, and (d) k=200, NRMSE = 0.13487.
(a)
(b)
(c)
(d)Considering the influence of initial value, we reconstruct the optical image with different initial guess. The retrieved results are shown in Figure5 when the initial values of the diffusion coefficient are D0=0.5 cm2 ns−1, D0=1.0 cm2 ns−1, D0=1.5 cm2 ns−1, and D0=2.0 cm2 ns−1.Figure 5
The reconstructed diffusion coefficients with MAP for different initial distributions. (a)D0=0.5 cm2 ns−1, (b) D0=1.0 cm2 ns−1, (c) D0=1.5 cm2 ns−1, and (d) D0=2.0 cm2 ns−1.
(a)
(b)
(c)
(d)As shown in Figure5, the reconstruction results with MAP are close to the real distribution of diffusion coefficient when the distribution of initial guess is close to the background medium. Compared with the reconstruction image without MAP as shown in Figure 3(b), it can be seen from Figure 5(b) that the shape of the inclusions reconstructed by the MAP estimation algorithm is closer to the original distribution and the image is smoother in the entire domain, although the retrieval values of the diffusion coefficient of A and B are smaller than the true values.The predicted values obtained without consideration of MAP algorithm are more correct than that with the MAP algorithm, but the reconstructed results with MAP can predict the shape correctly that is more important for optical image reconstruction. The regularization term based on MAP is the least squares function of optical parameters in neighbor points, as shown in (11); the minimization of objective function will diminish the difference of optical parameters in neighbor points. So the reconstructed images with MAP algorithm become smoother and the predicted magnitude shows deviations from the true magnitude.


### 3.2. Image Reconstruction of the Absorption Coefficient and Reduced Scattering Coefficient

Considering that both the absorption coefficient and scattering coefficient are inhomogeneous in the real biological tissues, the distribution of the absorption and reduced scattering coefficients was reconstructed simultaneously.Assume that the absorption coefficient and reduced scattering coefficient of background areμa=0.05 cm−1 and μs′=4.5 cm−1; A, B, and C are heterogeneous bodies whose absorption coefficients are 0.02 cm−1, 0.02 cm−1, and 0.1 cm−1, respectively, and reduced scattering coefficients are 3.0 cm−1, 3.0 cm−1, and 1.0 cm−1, respectively. Initial values of the absorption coefficient and reduced scattering coefficient are uniform as μa0=0.05 cm−1 and μs0′=4.5 cm−1. The distribution of the absorption coefficient and reduced scattering coefficient is reconstructed with MAP estimation algorithm. The real and reconstruction results for nonhomogeneous medium are shown in Figure 6.Figure 6
The distribution of absorption coefficient and reduced scattering coefficient in nonhomogeneous medium: (a) the true distribution of absorption coefficient, (b) the true distribution of reduced scattering coefficient, (c) the reconstructed distribution of absorption coefficient, and (d) the reconstructed distribution of reduced scattering coefficient.
(a)
(b)
(c)
(d)Figures6(c) and 6(d) show that the locations and distributions of the absorption coefficient and reduced scattering coefficient are well reconstructed simultaneously with MAP estimation algorithm. The algorithm proves to work well for the reconstruction of both absorption coefficient and scattering coefficient.


### 3.3. The Influence of the Measurement Error and Initial Distribution

To further illustrate the performance of the MAP estimation algorithm, random errors are considered. Measured values with random errors are obtained by adding normal distributed errors to the exact measured values, as shown in what follows [23]:(20)Ms,dn=M~s,dn+randnσs,dn,where Ms,dn and M~s,dn are the measured value and the exact measured value of the dth detector at the nth moment when the sth light source is working, respectively. randn is a standard normal distributed random number. The standard deviation of measured value σs,dn, for a measurement error of γ% at 99% confidence, is determined as(21)σs,dn=M~s,dn×γ%2.576,where 2.576 arises from the fact that 99% of a normally distributed population is contained within ±2.576 standard deviation of the mean.The influence of MAP is also investigated. Taking the model with three inclusions as example, the NRMSE values of the reconstruction results with the MAP estimation algorithm and without MAP estimation algorithm for different measurement errors and different initial distribution are shown in Table1.Table 1
The NRMSE values of the reconstruction results with the MAP estimation and without MAP for different measurement errors and different initial distributions.
Initial guess
γ%
With MAP
Without MAP
D
0
=
1.0 cm2 ns−1
0
0.14971
0.15285
5
0.18290
0.18232
10
0.18663
0.18907
D
0
=
1.5 cm2 ns−1
0
0.13487
0. 14928
5
0.17728
0.18274
10
0.18351
0.18815It can be seen from Table1 that the image quality of the reconstruction results deteriorates slightly with increasing measurement errors. The MAP estimation algorithm can improve the quality of reconstructed images with the decreasing NRMSE values. The initial distribution also affects the reconstruction results slightly. Even though the value of NRMSE is comparatively larger which is mainly caused by the deviation between the actual value and reconstructed value, the shape and location of inclusions are retrieved clearly in the reconstructed images (Figure 7) with measurement errors. From the viewpoint of real image reconstruction, the present algorithm has good robustness.Figure 7
The reconstructed images of diffusion coefficient with 10% measurement error. (a) Initial guessD0=1.0 cm2 ns−1; (b) initial guess D0=1.5 cm2 ns−1.
(a)
(b)


## 3.1. Image Reconstruction of Diffusion Coefficient

Assume that the absorption coefficient is constant asμa=0.05 cm−1, the diffusion coefficient of background is D=1.5 cm2 ns−1, and A, B, and C are heterogeneous bodies whose diffusion coefficients are 3.0 cm2 ns−1, 3.0 cm2 ns−1, and 1.0 cm2 ns−1, respectively. The absorption coefficient of inclusions is the same to background. Initial value is uniform and the diffusion coefficient is D0=1.0 cm2 ns−1. The real distribution and reconstruction results for nonhomogeneous medium are shown in Figure 3.Figure 3
The distribution of diffusion coefficient in nonhomogeneous medium: (a) the true distribution; (b) the reconstructed result (NRMSE = 0.17878).
(a)
(b)Figure3(b) shows that the locations and distributions of diffusion coefficient of three heterogeneous bodies are well reconstructed. The conjugate gradient method is used in the study of inverse problems. Since calculating the gradient of diffusion coefficient is needed, there are some step changes appearing in the retrieved results.Although the reconstruction peak of the diffusion coefficient is close to the real value, the shapes of three inclusions have been blurred in the computing process. So the MAP method is introduced to reconstruct the edge better. The variation of reconstructed images and NRMSE with number of iterations of the optimization algorithm with MAP is shown in Figure4. The shape becomes more accurate and the value of NRMSE becomes smaller with the increasing of iterations.Figure 4
The reconstructed diffusion coefficients with MAP for different number of iterations (k). (a) k=10, NRMSE = 0.22287, (b) k=50, NRMSE = 0.18365, (c) k=100, NRMSE = 0.15461, and (d) k=200, NRMSE = 0.13487.
(a)
(b)
(c)
(d)Considering the influence of initial value, we reconstruct the optical image with different initial guess. The retrieved results are shown in Figure5 when the initial values of the diffusion coefficient are D0=0.5 cm2 ns−1, D0=1.0 cm2 ns−1, D0=1.5 cm2 ns−1, and D0=2.0 cm2 ns−1.Figure 5
The reconstructed diffusion coefficients with MAP for different initial distributions. (a)D0=0.5 cm2 ns−1, (b) D0=1.0 cm2 ns−1, (c) D0=1.5 cm2 ns−1, and (d) D0=2.0 cm2 ns−1.
(a)
(b)
(c)
(d)As shown in Figure5, the reconstruction results with MAP are close to the real distribution of diffusion coefficient when the distribution of initial guess is close to the background medium. Compared with the reconstruction image without MAP as shown in Figure 3(b), it can be seen from Figure 5(b) that the shape of the inclusions reconstructed by the MAP estimation algorithm is closer to the original distribution and the image is smoother in the entire domain, although the retrieval values of the diffusion coefficient of A and B are smaller than the true values.The predicted values obtained without consideration of MAP algorithm are more correct than that with the MAP algorithm, but the reconstructed results with MAP can predict the shape correctly that is more important for optical image reconstruction. The regularization term based on MAP is the least squares function of optical parameters in neighbor points, as shown in (11); the minimization of objective function will diminish the difference of optical parameters in neighbor points. So the reconstructed images with MAP algorithm become smoother and the predicted magnitude shows deviations from the true magnitude.


## 3.2. Image Reconstruction of the Absorption Coefficient and Reduced Scattering Coefficient

Considering that both the absorption coefficient and scattering coefficient are inhomogeneous in the real biological tissues, the distribution of the absorption and reduced scattering coefficients was reconstructed simultaneously.Assume that the absorption coefficient and reduced scattering coefficient of background areμa=0.05 cm−1 and μs′=4.5 cm−1; A, B, and C are heterogeneous bodies whose absorption coefficients are 0.02 cm−1, 0.02 cm−1, and 0.1 cm−1, respectively, and reduced scattering coefficients are 3.0 cm−1, 3.0 cm−1, and 1.0 cm−1, respectively. Initial values of the absorption coefficient and reduced scattering coefficient are uniform as μa0=0.05 cm−1 and μs0′=4.5 cm−1. The distribution of the absorption coefficient and reduced scattering coefficient is reconstructed with MAP estimation algorithm. The real and reconstruction results for nonhomogeneous medium are shown in Figure 6.Figure 6
The distribution of absorption coefficient and reduced scattering coefficient in nonhomogeneous medium: (a) the true distribution of absorption coefficient, (b) the true distribution of reduced scattering coefficient, (c) the reconstructed distribution of absorption coefficient, and (d) the reconstructed distribution of reduced scattering coefficient.
(a)
(b)
(c)
(d)Figures6(c) and 6(d) show that the locations and distributions of the absorption coefficient and reduced scattering coefficient are well reconstructed simultaneously with MAP estimation algorithm. The algorithm proves to work well for the reconstruction of both absorption coefficient and scattering coefficient.


## 3.3. The Influence of the Measurement Error and Initial Distribution

To further illustrate the performance of the MAP estimation algorithm, random errors are considered. Measured values with random errors are obtained by adding normal distributed errors to the exact measured values, as shown in what follows [23]:(20)Ms,dn=M~s,dn+randnσs,dn,where Ms,dn and M~s,dn are the measured value and the exact measured value of the dth detector at the nth moment when the sth light source is working, respectively. randn is a standard normal distributed random number. The standard deviation of measured value σs,dn, for a measurement error of γ% at 99% confidence, is determined as(21)σs,dn=M~s,dn×γ%2.576,where 2.576 arises from the fact that 99% of a normally distributed population is contained within ±2.576 standard deviation of the mean.The influence of MAP is also investigated. Taking the model with three inclusions as example, the NRMSE values of the reconstruction results with the MAP estimation algorithm and without MAP estimation algorithm for different measurement errors and different initial distribution are shown in Table1.Table 1
The NRMSE values of the reconstruction results with the MAP estimation and without MAP for different measurement errors and different initial distributions.
Initial guess
γ%
With MAP
Without MAP
D
0
=
1.0 cm2 ns−1
0
0.14971
0.15285
5
0.18290
0.18232
10
0.18663
0.18907
D
0
=
1.5 cm2 ns−1
0
0.13487
0. 14928
5
0.17728
0.18274
10
0.18351
0.18815It can be seen from Table1 that the image quality of the reconstruction results deteriorates slightly with increasing measurement errors. The MAP estimation algorithm can improve the quality of reconstructed images with the decreasing NRMSE values. The initial distribution also affects the reconstruction results slightly. Even though the value of NRMSE is comparatively larger which is mainly caused by the deviation between the actual value and reconstructed value, the shape and location of inclusions are retrieved clearly in the reconstructed images (Figure 7) with measurement errors. From the viewpoint of real image reconstruction, the present algorithm has good robustness.Figure 7
The reconstructed images of diffusion coefficient with 10% measurement error. (a) Initial guessD0=1.0 cm2 ns−1; (b) initial guess D0=1.5 cm2 ns−1.
(a)
(b)


## 4. Conclusion

In conclusion, a MAP estimation based on Bayesian framework algorithm is applied to image reconstruction of two-dimensional highly scattering inhomogeneous medium. The time-resolved diffusion equation is solved using the FDM method, and for inverse solution of diffusion equation CG algorithm is used, in which the gradient of the objective function is calculated by the adjoint differentiation. GGMRF model can protect image edges. From the sense of spatial position of the defects, the reconstructed distribution of the optical parameters with the aid of MAP estimation algorithm is found to be more consistent with that of the original distribution than those with the unconstrained method. The image quality of the reconstruction results deteriorates with an increasing measurement error and the initial distribution also has a slight effect on the reconstruction results.



---

*Source: 412315-2015-03-25.xml*