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.
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 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]where denotes the location and denotes the time. indicates the medium. is the refractive index. is the speed of light in vacuum. is the incident radiation. is the radiative source term. is the absorption coefficient. is the diffusion coefficient, which is determined by where is called reduced scattering coefficient and is the scattering coefficient. is the scattering asymmetry parameter.
The initial and boundary conditions are expressed as follows.
Initial condition
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:where indicates the boundary. is the outward normal unit vector of the border. can be obtained bywhere and the critical angle .
By replacing the spatial and temporal derivatives with the fully implicit scheme of finite difference approximations, (1) can be expressed as follows:
For square grids (), the diffusion equation can be transformed into linear equations:
Equation (7) can be abbreviated as
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]where denotes the retrieval optical parameters. denotes the measured value. is the conditional probability density and 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 bywhere denotes the retrieval optical parameters. is the estimated value of the th detector at the th moment when the th light source is working. is the measured value of the th detector at the th moment when the th 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]where is the image sharpness parameter. is the scale parameter. is the set of all neighboring pixel pairs. is the neighboring position of . 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 aswhere is the sensitivity of with respect to and can be obtained by differentiating (8) with respect to , which yieldsThe intensity and matrix can be calculated by the forward model. The solution of (18) is similar to the solution forward model using the Bi-CGSTAB. is the derivative of objective function with respect to intensity and can be obtained recursively using the chain rulewithwhere is the sampling time which means the time consumed in the measurement of transmission and reflected signals. is obtained by partially differentiating objective function with respect to :
Differentiating (8) with respect to , we can obtainwith
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 Figure 1.
3. Results and Discussions
The computational model which contains three inhomogeneous media is shown in Figure 2. 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 , sampling time ns, and residual . The forward calculation model is cm2 with grids and 50 time steps.
To compare the overall accuracy of the reconstructions, the NRMSE is introduced and defined as [22]where is the retrieval value of the optical parameter at the th pixel. is the true value of the optical parameter at the th pixel. 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 cm−1, the diffusion coefficient of background is cm2 ns−1, and , , and are heterogeneous bodies whose diffusion coefficients are cm2 ns−1, cm2 ns−1, and cm2 ns−1, respectively. The absorption coefficient of inclusions is the same to background. Initial value is uniform and the diffusion coefficient is cm2 ns−1. The real distribution and reconstruction results for nonhomogeneous medium are shown in Figure 3.
(a)
(b)
Figure 3(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 Figure 4. The shape becomes more accurate and the value of NRMSE becomes smaller with the increasing of iterations.
(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 Figure 5 when the initial values of the diffusion coefficient are cm2 ns−1, cm2 ns−1, cm2 ns−1, and cm2 ns−1.
(a)
(b)
(c)
(d)
As shown in Figure 5, 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 and 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 cm−1 and cm−1; , , and are heterogeneous bodies whose absorption coefficients are cm−1, cm−1, and cm−1, respectively, and reduced scattering coefficients are cm−1, cm−1, and cm−1, respectively. Initial values of the absorption coefficient and reduced scattering coefficient are uniform as cm−1 and 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.
(a)
(b)
(c)
(d)
Figures 6(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]:where and are the measured value and the exact measured value of the th detector at the th moment when the th light source is working, respectively. is a standard normal distributed random number. The standard deviation of measured value , for a measurement error of % at 99% confidence, is determined aswhere 2.576 arises from the fact that 99% of a normally distributed population is contained within 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 Table 1.
It can be seen from Table 1 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.
(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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The supports of this work by the National Natural Science Foundation of China (no. 51476043), the Major National Scientific Instruments and Equipment Development Special Foundation of China (no. 51327803), and the Technological Innovation Talent Research Special Foundation of Harbin (no. 2014RFQXJ047) are gratefully acknowledged.