Image reconstruction of the absorption coefficients with l₁-norm minimization from photoacoustic measurements

Introduction

Biomedical imaging modalities, such as optical coherence tomography (1) and optical topography (2), have garnered attentions because of their abilities to provide structural and functional information in the tissue by taking advantages of the scattering and absorption of light while it propagates through biological media. Some promising clinical applications of optical imaging are being explored. Diffuse optical tomography (DOT) (3) has been developed for the non-invasive imaging of the cerebral blood concentration (4) and for the diagnoses of breast cancers (5). By applying the near-infrared spectroscopic technique, clinically useful information, such as oxygenation state of the blood, can be obtained by DOT. Fluorescent molecular imaging also exploits the light absorption, which can be used in small animals for pre-clinical tests (6).

In recent years, photoacoustic (PA) imaging (7-9) has been developed actively. PA imaging uses the thermoelastic wave, referred to as PA pressure from the photon absorbers, such as hemoglobin and exogenous contrast agents, excited by laser irradiation (10-12). The PA imaging can be applied to aid the diagnosis of cancers by providing the functional information about the increase of the blood volume caused by angiogenesis (13). PA pressure, which has a broad frequency band, is detected by the ultrasound transducer. The photon absorbers, which generate PA pressures, are imaged by the methods such as delay-and-sum projection or circular backprojection used in conventional medical ultrasound imaging (14). The image of the photon absorbers deep inside the biological medium can be obtained with high spatial resolution, compared with conventional optical imaging. Scattering of light by the tissues strongly attenuates the light intensity and decreases the spatial resolution of optical imaging. The ultrasound is also attenuated by the tissues, but the attenuation effect is relatively small.

The amplitude of PA pressure depends on the optical absorption coefficient of photon absorbers. So, the absorption coefficient can be estimated from the PA pressure by considering light propagation. The absorption coefficient depends on molar extinction coefficient and the concentration of photon absorber. By reconstructing the distribution of absorption coefficient, the concentration of targeted absorption coefficient is estimated (15-17). The image reconstruction of the absorption coefficient is carried out by minimizing the residual error between the measured and theoretically predicted PA signals. Therefore, the numerical calculation, such as finite element method (FEM) of the PA signal, is needed for the image reconstruction. Additionally, the amplitude of the PA pressure is not strictly linear to the absorption coefficient. As a result, the non-linear optimization scheme (18) is preferable in the image reconstruction. The precise prediction of the PA signals demands computing power. However, the non-linear optimization process also involves computational cost. Considering the actual clinical use, a quick and reliable image reconstruction is desired.

The linearization of the relation between the PA signals and absorption coefficient is effective to reduce the computational cost. However, the quality of the image, such as the spatial resolution and the quantification ability, may be compromised. In this study, a linearized PA image reconstruction is attempted in numerical simulation and phantom experiment. In the meantime, the l₁-norm of the reconstructed image is also minimized in the image reconstruction process in order to improve the localization of the targeted photon absorber. The latter has been used in the inverse problems and image processing in other fields (19-22). The quality of the reconstructed image with the l₁-norm minimization is discussed to compare with the conventional Tikhonov regularization.

Materials and methods

Image reconstruction algorithm

To reconstruct the distribution of the absorption coefficient from the detected PA pressures, the fundamental equations dealing with the generation and the propagation of the PA pressure, which is the photon diffusion equation (PDE) and the PA wave equation, were applied. The light is scattered and absorbed by the tissues in the biological media, and the propagation of the light is the radiative transfer of the photon energy. Therefore, the radiative transfer equation (RTE) rigorously describes the light propagation (23,24). The following PDE obtained from the diffusion approximation of the RTE (25) was used in this study to reduce the calculation cost,

where D=1/(3μ_s’) is the diffusion coefficient with the reduced scattering coefficient μ_s’, μ_a, the absorption coefficient, Φ, the fluence rate, q₀, the light source, and r, the position. The robin boundary condition is usually used with the PDE,

where n is the vector outer normal to the surface of the medium, and A is the parameter depending on the refractive index of the medium. Eq. [1] can be solved for Φ by FEM, and the PA source H, which is the absorbed photon energy, is calculated as H(r)=μ_a(r)Φ(r).

On the other hand, the propagation of the PA pressure in acoustically homogeneous medium is described by the PA wave equation (8,9,15),

where v is the speed of the PA pressure, t, time, p, the PA pressure, and Γ, the Grüneisen parameter representing the efficiency of the PA pressure generation. Eq. [3] also can be solved by the FEM.

Let us assume that the PA pressure is detected by the PA probe which irradiates the light and detects the PA pressure at the identical position, and that the PA pressure generated directly below the PA probe is detected. The PA pressures with T time samples are detected at the K positions. According to Eq. [2], p is linear to H, the PA pressure, which is m_k of the T-vector, detected at the kth position is formulated by discretizing the functions of t and r with FEM as,

where H_k is the S-vector which consists of the absorbed photon energy at S discrete positions below the PA probe, and L_k is the T×S-matrix which projects H_k to m_k.

H_k is linearly approximated as, , where is the energy absorbed by the background medium with the absorption coefficient , the S×N-matrix J_k consists of the differential coefficients of H_k about , and ∆μ_a is the N-vector of the changes in μ_a at the N discrete positions in the whole medium. To reconstruct the PA source, ∆μ_a was reconstructed in this study. By assuming the medium is large enough to say L_j≈L_k (j≠k), and by subtracting m_l form m_k, we obtain (26),

then ∆μ_a is reconstructed by solving the following minimization problem,

where ∆m is the CT-vector consists of ∆m_k,j with C of the number of the combinations of k and j. G is the CT×N matrix consists of G_k,j. f is a regularization function, and λ is the regularization parameter to adjust the effect of the regularization. When is given, the distribution of μ_a is calculated as , where is the reconstructed ∆μ_a.

In this study, we compared two regularization functions: l₁-norm and l₂-norm of the reconstructed ∆μ_a. The function f employing the l_p-norm is described as where ∆μ_a,i is the i th component of ∆μ_a (27).

Minimization of the l₁-norm generally provides sparse distribution of the reconstructed solution of the inverse problems expressed by the same manner of Eq. [6]. In this study, it was expected that the photon absorber with ∆μ_a larger than the background was localized. In the minimization process, ∆μ_a,i and f were reformulated as , respectively (26). The optimization was carried out with the non-linear conjugate gradient method.

On the other hand, minimization of the l₂-norm, which is usually referred to as Tikhonov regularization (28), provides the smooth solution by reducing the influence of noise. The solution with the l₂-norm minimization was obtained as,

where I is the identity matrix.

We tried the image reconstruction with following three regularizations: the l₁-norm minimization with λ at the corner of the L-curve (Reg.1), the l₂-norm minimization with λ at the plateau of the L-curve (Reg.2), and the l₂-norm minimization with λ at the corner of the L-curve (Reg.3). Reg.2 used unusual selection of λ. The L-curve can be regarded as the plot of f as a function of the squared error e in Eq. [6]. By selecting the position where |df⁄de| is minimized, it was expected that Reg.2 would provide the smoother distribution of ∆μ_a than Reg.3, because Reg.2 took larger λ than Reg.3.

Conditions of the numerical simulation

Figure 1A shows the geometrical conditions for the numerical simulation. The medium of the PA pressure was a 2D square with 50 mm side. The medium had the background optical properties of μ_s’ =0.8 mm^–1 and . The measurement of the PA pressure was conducted at eleven positions with an equal spacing of 2 mm on the surface of the medium. There existed a photon absorber with 2 mm side and μ_a =0.6, 1.1 or 1.7 mm^–1 on x=0 mm. The depth at which the photon absorber is located was y=5, 7 or 9 mm from the surface of the medium. The measured data ∆m used in the inversion process were calculated by the use of Eq. [5].

Figure 1 Geometrical conditions of (A) the numerical simulation and (B) the phantom experiment. PA, photoacoustic; ICG, indocyanine green.

The 2D image reconstruction was carried out. The positions of the irradiation and the detection were identical in the image reconstruction process. The matrices of L_k and H_k were calculated with FEM with Eqs. [1] and [2]. FEM employed 10,201 nodes and 20,000 triangular elements. Gaussian noise with the standard deviation (SD) of 10% of the maximum of ∆m was added. The image reconstruction was carried out on pixel basis as described in the literature (29). A single pixel consisted of 25 nodes. Three trials of the image reconstruction with Regs.1, 2 or 3 were carried out for each combination of μ_a and the depth of the photon absorber.

Conditions of the phantom experiment

Figure 1B shows the schema of the phantom experiments. The background of the phantom was an aqueous solution of the intralipid and indocyanine green (ICG). The optical properties of the phantom was adjusted to μ_s’ =0.8 mm^–1 and (30). The tube with an inner diameter of 1 mm and an outer diameter of 2 mm placed in the phantom as a photon absorber. The tube contained the intralipid and ICG with μ_s’ =0.8 mm^–1 and μ_a =0.6, 1.1 or 1.7 mm^–1. The depth of the photon absorber was y=5, 7 or 9 mm.

A tunable Ti:sapphire laser pumped by the second harmonic of a Q-switch Nd:YAG laser (LT-221 and LS-2134, Lotis Tii, Minsk, Belarus) was used for the laser irradiation. The laser light was introduced into the optical fiber. The optical fiber was installed in the cylindrical PA probe which had a ring shaped piezoelectric film P(VDF-TrFE) (KF piezo-film, Kureha Corp., Tokyo, Japan) on the edge of the PA probe (Figure 2). The edge of the optical fiber was placed at the center of the piezoelectric film so that the laser irradiation and the PA detection were conducted at the identical position. The energy of the light from the optical fiber was 4 mJ/pulse. The PA pressure was detected at the eleven positions aligned perpendicularly to the long axis of the tube of the photon absorber with an equal spacing of 2 mm. The measured data were acquired by the digital oscilloscope (DSO8104A, Agilent Technologies, Santa Clara, CA) via the amplifier (SA-220F5, NF Corp., Yokohama, Japan). The image reconstruction was carried out in 2D with the same manner of the numerical simulation.

Figure 2 The edge of the PA probe used in the phatom experiment. PA, photoacoustic.

Results

Numerical simulation

Figure 3 shows the reconstructed images using the regularization methods in the case of μ_a =1.1 mm^–1 and the depth of 7 mm of the photon absorber. Reg.1 reconstructed the localized μ_a distribution. The reconstructed value and the positions of the photon absorber were corrected in all cases with a combination of μ_a and depths. Reg.2 also reconstructed the photon absorber at the correct position, although the reconstructed value of μ_a was smaller and the distribution was slightly broader than the true one. The influence of the noise was not seen in the images reconstructed with Regs.1 and 2.

Figure 3 Reconstructed images of μ_a (mm^–1) with (A) Reg.1 (B) Reg.2 and (C) Reg.3 in the numerical simulation.

Reg.3 reconstructed the maximum value of μ_a at the correct position of the photon absorber. Although Reg.3 used Tikhonov regularization, the influence of the noise was seen at the deep positions in the medium. The reconstructed images in all cases with various μ_a and the depths of the photon absorber had similar characters depending on the regularization methods.

Figure 4 shows the signal-to-noise ratio (SNR) of the reconstructed μ_a of the photon absorber. SNR₁ was defined as SNR₁ =10 log₁₀ (∆μ_a,max /σ₁)², where ∆μ_a,max is the average of the maximum reconstructed ∆μ_a and σ₁ is SD of ∆μ_a,max estimated from the three trials with each of Regs1, 2 and 3. Figure 4A shows the averages of SNR₁ among the cases with the depths of 5, 7, 9 mm which were calculated for each of the cases with μ_a =0.6, 1.1 or 1.7 mm^–1. The averages of SNR₁ among the cases with μ_a =0.6, 1.1 or 1.7 mm^–1 calculated for the cases with the depths of 5, 7 or 9 mm are shown in Figure 4B. SNR₁ were in the range from 30 to 40 dB in all cases. This means that ∆μ_a,max was approximately 60 times larger than the error due to the measurement noise. The dependences of SNR₁ on the depth, ∆μ_a,max of the photon absorber or the regularization methods were not seen in Figure 4.

Figure 4 Averages of SNR₁ of the reconstructed images for (A) the cases with μ_a =0.6 mm^–1 (red), 1.1 mm^–1 (green) and 1.7 mm^–1 (blue), and (B) the cases with the depths of 5 mm (red), 7 mm (green) and 9 mm (blue).

On the other hand, Figure 5 shows SNR₂ defined as SNR₂ =10 log₁₀ (∆μ_a,max/σ₂)², where σ₂ was the total of the reconstructed |∆μ_a| except ∆μ_a,max. In the same manner as SNR₁, the averages of SNR₂ calculated for the cases with μ_a =0.6, 1.1 or 1.7 mm^–1 are shown in Figure 5A, and those calculated for the depth of 5, 7 or 9 mm are in Figure 5B. Reg.1 obtained SNR₂ of about 25 dB in all cases, while SNR₂ of Regs.2 and 3 were about –10 and –15 dB, respectively. ∆μ_a of the photon absorber reconstructed with Reg.1 was 20 times larger than the artifacts in the reconstructed image caused by the noisy measurement data, while the artifacts were larger than the reconstructed ∆μ_a of the photon absorber in the images reconstructed with Regs.2 and 3.

Figure 5 Averages of SNR₂ of the reconstructed images for (A) the cases with μ_a =0.6 mm^–1 (red), 1.1 mm^–1 (green) and 1.7 mm^–1 (blue), and (B) the cases with the depths of 5 mm (red), 7 mm (green) and 9 mm (blue).

Phantom experiment

Figure 6 shows the reconstructed images in the case with μ_a =1.1 mm^–1 and the depth of 7 mm of the photon absorber of the phantom experiments. The reconstructed values were calibrated by the regression line fitting the give sets of the true μ_a of the photon absorber and the averages of the reconstructed values at the depths of 5, 7 and 9 mm. The maximum value of μ_a was reconstructed at the correct position for all combinations of μ_a and the depths of the photon absorber by each of the regularization methods. The reconstructed images demonstrated that Reg.1 reduced the influence of noise and localized μ_a distribution more strongly than Regs.2 and 3 as well as in the previous numerical simulation. The quality of the reconstructed image was improved by the l₁-norm minimization, although the reconstructed photon absorber had the area larger than true one. SNR₂ of the images reconstructed with Regs.1, 2 and 3 were –12.4, –23.0 and –31.3 dB, respectively. The phantom experiment showed that Regs. 1 can obtain larger SNR₂ than Regs.2 and 3.

Figure 6 Reconstructed images of μ_a (mm^–1) with (A) Reg.1 (B) Reg.2 and (C) Reg.3 in the phantom experiment.

μ_a of the photon absorber reconstructed with Reg.1 is plotted as a function of the true μ_a of the photon absorber for the numerical simulation and the phantom experiment in Figure 7A,B, respectively. In both of the numerical simulation and phantom experiment, the reconstructed μ_a of the photon absorber reflected the increase of the true μ_a with linearity on some level. Regs.2 and 3 also showed the similar results. However, μ_a,max depended on the depth of the photon absorber in the phantom experiment. The deeper position the photon absorber was placed at, the smaller the reconstructed μ_a of the photon absorber was, while μ_a,max did not depend on the depth of the photon absorber in the numerical simulation.

Figure 7 μ_a reconstructed with Reg.1 as a function of the true μ_a in (A) the numerical simulation and (B) the phantom experiment when the depth was 5 mm (circle symbol), 7 mm (square symbol) or 9 mm (triangle symbol).

Discussion

Artifacts appeared in the images reconstructed with Reg.3 in both the numerical simulation and the phantom experiment. The image reconstruction of μ_a from the PA pressure took the light propagation in the medium into account. At the deeper position, Φ became small owing to scattering and absorption of the light. To reconstruct μ_a of the photon absorber in various depths correctly, the image reconstruction algorithm compensated the rapid decrement of Φ. This function of the image reconstruction algorithm amplified the measurement noise observed at the late times. The time at which the PA pressure was detected corresponds to the distance from the detector to the photon absorber. So, the influence of the measurement noise observed at the late times appeared as artifacts at the deeper positions. The artifacts were reduced by Reg.2, because λ was larger, and the smoothing effect of the Tikhonov regularization of Reg.2 was stronger than that of Reg.3.

The photon absorber was reconstructed with large SNR₁ at the correct position by Regs.1, 2 and 3, regardless of μ_a and the depths of the photon absorber, although SNR₁ fluctuated owing to the measurement noise. SNR₂ quantitatively shows that Reg.1 strongly suppressed the influence of measurement noise on the reconstructed image, regardless of μ_a and the depths of the photon absorber. Reg.1 was superior to Regs.2 and 3 in the sense of SNR₂. Reg.2 is better to be used to reduce the influence of measurement noise than Reg.3 which used the conventional parameter selection for the Tikhonov regularization. SNR₂ in the phantom experiment was smaller than that in the numerical simulation, because the reconstructed large μ_a distributed larger than true one and some artifact appeared owing to noise.

Through the numerical simulation and the phantom experiment, Reg.1 obtained the reconstructed image effectively suppressing the influence of the measurement noise, because the reconstructed distribution of ∆μ_a needed to be sparse to minimize the l₁-norm. It is known that the regularization minimizing the l_p-norm with p smaller than unity generally obtains sparse solutions of inverse problems (19,20). It was demonstrated that the l₁-norm minimization provided the sparse solution robust to noise by the evaluations of the reconstructed images with SNR₁ and SNR₂. The characteristic of the l₁-norm minimization in obtaining the sparse solution can be useful to image small imaging targets, such as micro-blood vessels and the cancers at early stage labeled by the contrast dye in a small region. The robustness to noise of the image reconstruction with the l₁-norm minimization can be helpful for the image reconstruction of small imaging targets, because the PA pressure from the small imaging target is small and the SNR of the PA pressure is low.

Some side effects of the powerful regularization should be noticed. Regs.1 and 2 may suppress the smaller PA signals from the photon absorber embedded deeply inside the medium. Additionally, the strong localization effect of Reg.1 can reconstruct the distribution of the photon absorber smaller than true ones when the true photon absorbers are broadly distributed. The PA image reconstructed with a certain regularization method should be carefully interpreted for the medical diagnosis by considering the characteristics of the regularization method.

Figure 6A shows that Reg.1 quantitatively reconstructed μ_a of the photon absorber when the image reconstruction algorithm fit the physical model of the PA measurement. The difference between Figure 6A,B can occur because of the 2D image reconstruction used in this study. In the phantom experiment, the light and the PA pressure propagated in larger volume of 3D medium than that of 2D medium. Then, the excitation light and PA pressure from deeper position in 3D medium was attenuated more strongly than that predicted in 2D medium.

However, μ_a at the depth of 5 mm was twice as large as that of 9 mm, while the PA pressure from the depth of 5 mm was about 5-fold of that of 9 mm in the phantom experiment. This means that the 2D image reconstruction compensated the attenuation of the light and PA pressure to some extent, although the 2D image reconstruction cannot recover them perfectly.

Conclusions

In this study, the effect of the l₁-norm minimization on the image reconstruction from the PA measurement was compared to the Tikhonov regularization. By the numerical simulation and the phantom experiment, we demonstrated that the l₁-norm minimization localized the distribution of the absorption coefficient and clearly imaged the photon absorber in the medium. The influence of noise on the reconstructed image was sensible by using the Tikhonov regularization in this study. The image reconstruction with the regularization minimizing the l₁-norm is useful for the PA measurement when the distribution of the photon absorber is sparse and the signal-to-noise ratio of the detected PA pressure is low. The l₁-norm minimization can be used for the nonlinear image reconstruction algorithm.

Acknowledgements

This work was partially supported by JST collaborative Research Based on Industrial Demand (In vivo Molecular Imaging: Towards Biophotonics Innovations in Medicine) and JSPS KAKENHI (Grant Number 24760314). The authors thank Hiroaki Ishihara for his assistance in the phantom experiment.

Authors’ contribution: S.O. designed overall study and wrote the paper. S.O. constructed the image reconstruction algorithm and carried out the numerical simulations. S.O., T.H. and T.K. constructed the experimental setup, collected and analyzed the data. S.O., T.H., T.K. and M.I. discussed and edited the paper. M.I. supervised this study.

Disclosure: The authors declare no conflict of interest.

References

1. Yin X, Chao JR, Wang RK. User-guided segmentation for volumetric retinal optical coherence tomography images. J Biomed Opt 2014;19:086020. [PubMed]
2. Wang S, Shibahara N, Kuramashi D, Okawa S, Kakuta N, Okada E, Maki A, Yamada Y. Effect of Spatial Variation of Skull and Cerebrospinal Fluid Layer on Optical Mapping of Brain Activities. Opt Rev 2010;17:410-20.
3. Gibson AP, Hebden JC, Arridge SR. Recent advances in diffuse optical imaging. Phys Med Biol 2005;50:R1-43. [PubMed]
4. Fukuzawa R, Okawa S, Matsuhashi S, et al. Reduction of image artifacts induced by change in the optode coupling in time-resolved diffuse optical tomography. J Biomed Opt 2011;16:116022. [PubMed]
5. Yates T, Hebden JC, Gibson A, Everdell N, Arridge SR, Douek M. Optical tomography of the breast using a multi-channel time-resolved imager. Phys Med Biol 2005;50:2503-17. [PubMed]
6. Okawa S, Yano A, Uchida K, Mitsui Y, Yoshida M, Takekoshi M, Marjono A, Gao F, Hoshi Y, Kida I, Masamoto K, Yamada Y. Phantom and mouse experiments of time-domain fluorescence tomography using total light approach. Biomed Opt Express 2013;4:635-51. [PubMed]
7. Wang LV, Hu S. Photoacoustic tomography: in vivo imaging from organelles to organs. Science 2012;335:1458-62. [PubMed]
8. Li C, Wang LV. Photoacoustic tomography and sensing in biomedicine. Phys Med Biol 2009;54:R59-97. [PubMed]
9. Xu M, Wang LV. Photoacoustic imaging in biomedicine. Rev Sci Instrum 2006;77:041101.
10. Hu S, Wang LV. Photoacoustic imaging and characterization of the microvasculature. J Biomed Opt 2010;15:011101. [PubMed]
11. Erpelding TN, Kim C, Pramanik M, Jankovic L, Maslov K, Guo Z, Margenthaler JA, Pashley MD, Wang LV. Sentinel lymph nodes in the rat: noninvasive photoacoustic and US imaging with a clinical US system. Radiology 2010;256:102-10. [PubMed]
12. Wang X, Pang Y, Ku G, Xie X, Stoica G, Wang LV. Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain. Nat Biotechnol 2003;21:803-6. [PubMed]
13. Heijblom M, Piras D, Xia W, van Hespen JC, Klaase JM, van den Engh FM, van Leeuwen TG, Steenbergen W, Manohar S. Visualizing breast cancer using the Twente photoacoustic mammoscope: what do we learn from twelve new patient measurements? Opt Express 2012;20:11582-97. [PubMed]
14. Sperl JI, Zell K, Menzenbach P, Haisch C, Ketzer S, Marquart M, Koenig H, Vogel MW. Photoacoustic image reconstruction: a quantitative analysis. Proc SPIE 2007;6631:663103.
15. Cox B, Laufer JG, Arridge SR, Beard PC. Quantitative spectroscopic photoacoustic imaging: a review. J Biomed Opt 2012;17:061202. [PubMed]
16. Laufer J, Cox B, Zhang E, Beard P. Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme. Appl Opt 2010;49:1219-33. [PubMed]
17. Yuan Z, Wang Q, Jiang H. Reconstruction of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton method. Opt Express 2007;15:18076-81. [PubMed]
18. Saratoon T, Tarvainen T, Cox BT, Arridge SR. A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation. Inverse Problems 2013;29:075006.
19. Okawa S, Hoshi Y, Yamada Y. Improvement of image quality of time-domain diffuse optical tomography with l sparsity regularization. Biomed Opt Express 2011;2:3334-48. [PubMed]
20. Okawa S, Ikehara T, Oda I, Yamada Y. Reconstruction of localized fluorescent target from multi-view continuous-wave surface images of small animal with lp sparsity regularization. Biomed Opt Express 2014;5:1839-60. [PubMed]
21. Xu P, Tian Y, Chen H, Yao D. Lp norm iterative sparse solution for EEG source Localization. IEEE Trans Biomed Eng 2007;54:400-9. [PubMed]
22. Chang CH, Ji JX. Improving multi-channel compressed sensing MRI with reweighted l 1 minimization. Quant Imaging Med Surg 2014;4:19-23. [PubMed]
23. Tarvainen T, Kolehmainen V, Arridge SR, Kaipio JP. Image reconstruction in diffuse optical tomography using the coupled radiative transport-diffusion model. J Quant Spectrosc Radiat Transfer 2011;112:2600-8.
24. Fujii H, Okawa S, Yamada Y, Hoshi Y. Hybrid model of light propagation in random media based on the time-dependent radiative transfer and diffusion equations. J Quant Spectrosc Radiat Transfer 2014;147:145-54.
25. Arridge SR. Optical tomography in medical imaging. Inverse Problems 1999;15:R41-93.
26. Okawa S, Hirasawa T, Kushibiki T, Ishihara M. Reconstruction of the optical properties of inhomogeneous medium from photoacoustic signal with lp sparsity regularization. Proc SPIE 2013;8581:858135.
27. He Z, Cichocki A, Zdunek R, Xie S. Improved FOCCUS Method With Conjugate Gradient Iterations. IEEE Trans Signal Process 2009;57:399-404.
28. Vogel CR. Computational Methods for Inverse Problems (Frontiers in Applied Mathematics), Philadelphia, 2002.
29. Okawa S, Hirasawa T, Kushibiki T, Ishihara M. Numerical Evaluation of Linearized Image Reconstruction Based on Finite Element Method for Biomedical Photoacoustic Imaging. Opt Rev 2013;20:442-51.
30. van Staveren HJ, Moes CJ, van Marie J, Prahl SA, van Gemert MJ. Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm. Appl Opt 1991;30:4507-14. [PubMed]