Mathematical Problems in Engineering
Volume 2015 (2015), Article ID 412315, 9 pages
http://dx.doi.org/10.1155/2015/412315

Research Article

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

Hong Qi, Yaobin Qiao, Shuangcheng Sun, Yuchen Yao, and Liming Ruan

School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, China

Received 5 July 2014; Revised 11 November 2014; Accepted 15 November 2014

Academic Editor: Subhash Chandra Mishra

Copyright © 2015 Hong Qi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 by

where 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 by

where 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 as

where is the sensitivity of with respect to and can be obtained by differentiating (8) with respect to , which yields

The 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 rule

with

where 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 obtain

with

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.

Figure 1: The flowchart of the image reconstruction.

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  cm² with 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]

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⁻¹, the diffusion coefficient of background is  cm² ns⁻¹, and , , and are heterogeneous bodies whose diffusion coefficients are  cm² ns⁻¹,  cm² ns⁻¹, and  cm² ns⁻¹, respectively. The absorption coefficient of inclusions is the same to background. Initial value is uniform and the diffusion coefficient is  cm² ns⁻¹. The real distribution and reconstruction results for nonhomogeneous medium are shown in Figure 3.

(a)

(b)

Figure 3: The distribution of diffusion coefficient in nonhomogeneous medium: (a) the true distribution; (b) the reconstructed result (NRMSE = 0.17878).

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)

Figure 4: The reconstructed diffusion coefficients with MAP for different number of iterations (). (a) , NRMSE = 0.22287, (b) , NRMSE = 0.18365, (c) , NRMSE = 0.15461, and (d) , NRMSE = 0.13487.

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  cm² ns⁻¹,  cm² ns⁻¹,  cm² ns⁻¹, and  cm² ns⁻¹.

(a)

(b)

(c)

(d)

Figure 5: The reconstructed diffusion coefficients with MAP for different initial distributions. (a)  cm² ns⁻¹, (b)  cm² ns⁻¹, (c)  cm² ns⁻¹, and (d)  cm² ns⁻¹.

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⁻¹ and  cm⁻¹; , , and are heterogeneous bodies whose absorption coefficients are  cm⁻¹,  cm⁻¹, and  cm⁻¹, respectively, and reduced scattering coefficients are  cm⁻¹,  cm⁻¹, and  cm⁻¹, respectively. Initial values of the absorption coefficient and reduced scattering coefficient are uniform as  cm⁻¹ and  cm⁻¹. 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)

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.

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 as

where 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.

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

 cm² ns⁻¹
0
0.14971
0.15285

5
0.18290
0.18232

10
0.18663
0.18907

 cm² ns⁻¹
0
0.13487
0. 14928

5
0.17728
0.18274

10
0.18351
0.18815

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)

Figure 7: The reconstructed images of diffusion coefficient with 10% measurement error. (a) Initial guess  cm² ns⁻¹; (b) initial guess  cm² ns⁻¹.

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.

References

B. A. Brooksby, H. Dehghani, B. W. Poque, and K. D. Paulsen, “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities,” IEEE Journal on Selected Topics in Quantum Electronics, vol. 9, no. 2, pp. 199–209, 2003.

Z. C. Wang, J. L. Song, J. L. Chai, Q. Cheng, and H. C. Zhou, “The effect of BRDF surface on radiative heat transfer within a one-dimensional graded index medium,” International Journal of Thermal Sciences, vol. 77, pp. 116–125, 2014.

J. Boulanger and A. Charette, “Numerical developments for short-pulsed Near Infra-Red laser spectroscopy. Part I. Direct treatment,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 91, no. 2, pp. 189–209, 2005.

Z. Wang, Q. Cheng, and H. Zhou, “The DRESOR method for transient radiation transfer in 1-D graded index medium with pulse irradiation,” International Journal of Thermal Sciences, vol. 68, pp. 127–135, 2013.

Q. Cheng and H.-C. Zhou, “The DRESOR method for a collimated irradiation on an isotropically scattering layer,” Journal of Heat Transfer, vol. 129, no. 5, pp. 634–645, 2007.

H. J. Niu, P. Guo, X. D. Song, and T. Jiang, “Improving depth resolution of diffuse optical tomography with an exponential adjustment method based on maximum singular value of layered sensitivity,” Chinese Optics Letters, vol. 6, no. 12, pp. 886–888, 2008.

F. Wang, J. Tan, and Z. Wang, “Heat transfer analysis of porous media receiver with different transport and thermophysical models using mixture as feeding gas,” Energy Conversion and Management, vol. 83, pp. 159–166, 2014.

J. Hao, L. Ruan, and H. Tan, “Effect of anisotropic scattering on radiative heat transfer in two-dimensional rectangular media,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 78, no. 2, pp. 151–161, 2003.

Y. Yamada and Y. Hasegawa, “Time-dependent FEM analysis of photon migration in biological tissues,” JSME International Journal Series B: Fluids and Thermal Engineering, vol. 39, no. 4, pp. 754–761, 1996.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Medical Physics, vol. 22, no. 11, pp. 1779–1792, 1995.

M. Schweiger, S. R. Arridge, and D. T. Delpy, “Application of the finite-element method for the forward and inverse models in optical tomography,” Journal of Mathematical Imaging and Vision, vol. 3, no. 3, pp. 263–283, 1993.

M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Physics in Medicine and Biology, vol. 50, no. 12, pp. 2837–2858, 2005.

Y. Wang and P. Zuo, “A variational model based on maximum posterior probability for restoration of MR images corrupted by Rician noise,” Journal of Jilin University (Science Edition), vol. 51, pp. 289–302, 2013.

J. Shi, B. Zhang, F. Liu, J. Luo, and J. Bai, “Efficient L1 regularization-based reconstruction for fluorescent molecular tomography using restarted nonlinear conjugate gradient,” Optics Letters, vol. 38, no. 18, pp. 3696–3699, 2013.

G. Zhang, X. Cao, B. Zhang, F. Liu, J. Luo, and J. Bai, “MAP estimation with structural priors for fluorescence molecular tomography,” Physics in Medicine and Biology, vol. 58, no. 2, pp. 351–372, 2013.

X. Intes, C. Maloux, M. Guven, B. Yazici, and B. Chance, “Diffuse optical tomography with physiological and spatial a priori constraints,” Physics in Medicine and Biology, vol. 49, no. 12, pp. N155–N163, 2004.

K. Xu, F. Gao, and H. Zhao, Biomedical Photonics, Science Press, Beijing, China, 2nd edition, 2011.

H. A. van der Vorst, “Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems,” Society for Industrial and Applied Mathematics. Journal on Scientific and Statistical Computing, vol. 13, no. 2, pp. 631–644, 1992.

P. Carvalho, A. Santos, A. Dourado et al., “Bayes information criterion for Tikhonov regularization with linear constraints: application to spectral data estimation,” in Proceedings of the 16th International Conference of Pattern Recognition Proceedings, vol. 1, pp. 696–700, 2002.

J. O. Berger, Statistical Decision Theory and Bayesian Analysis, Springer Series in Statistics, Springer, New York, NY, USA, 2nd edition, 1985.

S. S. Saquib, K. M. Hanson, and G. S. Cunningham, “Model-based image reconstruction from time-resolved diffusion data,” in Medical Imaging 1997: Image Processing, Proceedings of SPIE, pp. 369–380, Newport Beach, Calif, USA, February 1997.

J. C. Ye, K. J. Webb, C. A. Bouman, and R. P. Millane, “Optical diffusion tomography by iterative-coordinate-descent optimization in a Bayesian framework,” Journal of the Optical Society of America, vol. 16, no. 10, pp. 2400–2412, 1999.

B. Zhang, H. Qi, Y. T. Ren, S. C. Sun, and L. M. Ruan, “Application of homogenous continuous Ant Colony Optimization algorithm to inverse problem of one-dimensional coupled radiation and conduction heat transfer,” International Journal of Heat and Mass Transfer, vol. 66, pp. 507–516, 2013.