Effects of the approximations of light propagation on quantitative photoacoustic tomography using two-dimensional photon diffusion equation and linearization

Quantitative photoacoustic tomography (QPAT) employing a light propagation model will play an important role in medical diagnoses by quantifying the concentration of hemoglobin or a contrast agent. However, QPAT by the light propagation model with the three-dimensional (3D) radiative transfer equation (RTE) requires a huge computational load in the iterative forward calculations involved in the updating process to reconstruct the absorption coefficient. The approximations of the light propagation improve the efficiency of the image reconstruction for the QPAT. In this study, we compared the 3D/two-dimensional (2D) photon diffusion equation (PDE) approximating 3D RTE with the Monte Carlo simulation based on 3D RTE. Then, the errors in a 2D PDE-based linearized image reconstruction caused by the approximations were quantitatively demonstrated and discussed in the numerical simulations. It was clearly observed that the approximations affected the reconstructed absorption coefficient. The 2D PDE-based linearized algorithm succeeded in the image reconstruction of the region with a large absorption coefficient in the 3D phantom. The value reconstructed in the phantom experiment agreed with that in the numerical simulation, so that it was validated that the numerical simulation of the image reconstruction predicted the relationship between the true absorption coefficient of the target in the 3D medium and the reconstructed value with the 2D PDE-based linearized algorithm. Moreover, the the true absorption coefficient in 3D medium was estimated from the 2D reconstructed image on the basis of the prediction by the numerical simulation. The estimation was successful in the phantom experiment, although some limitations were revealed.


Introduction
Photoacoustic (PA) imaging [1][2][3] noninvasively provides physiological information such as hemodynamics in micro blood vessels. By employing the light propagation model described by the radiative transfer equation (RTE) or the photon diffusion equation (PDE) [4][5][6][7], the PA image of the light absorption coefficient can be reconstructed quantitatively. This technique is referred to as quantitative PA tomography (QPAT) [8]. The absorption coefficient of the photon absorber can reflect the conditions of biological tissues. Angiogenesis in cancer tissues increases the concentration of hemoglobin, which increases the absorption coefficient [9][10][11][12][13][14][15]. Some contrast agent molecules and nanoparticles increasing the light absorption coefficient are useful for finding diseased tissues [16,17] and mapping the sentinel lymph nodes [18,19]. Therefore, QPAT reconstructing the absorption coefficient will be useful in medical diagnoses.
Various QPAT algorithms have been proposed. Cox et al. proposed the 2D QPAT reconstruction of the absorption coefficient based on PDE by the efficient iterative method [20]. The model-based inversion scheme, which was based on the PDE with the d-Eddington approximation, successfully reconstructed an image from the PA pressure wave detected by using a Fabry-Perot polymer film [21]. Yuan et al. reconstructed the absorption coefficient using PDE and the finite element method (FEM)-based PA image reconstruction algorithm [22]. Saratoon et al. used RTE [23]. Pulkkinen and coworkers investigated the image reconstruction based on RTE with illumination from a single direction [24].
One of the technical difficulties of QPAT is the forward calculation with a RTE in 3D medium, which requires computational load and time. Forward calculation simulates the PA pressure wave from given optical properties by solving RTE, and the absorption coefficient is reconstructed by minimizing the error between the simulated and measured PA pressure waves with an iterative updating process. Although RTE describes the light propagation rigorously, the numerical calculation of RTE requires the discretization of the directions of the light propagation aside from the spatial positions. Thus, huge computational load and time are needed for image reconstruction. Therefore, QPAT has been discussed mainly by using 2D light propagation model and the PDE approximating RTE [8].
We have developed a 2D PDE-based linearized algorithm for QPAT to reconstruct the absorption coefficient [25][26][27][28][29]. The algorithm reconstructed the image from a small number of PA pressure waves measured at a single side of the media by using a probe for illumination and detection at an identical position. The algorithm took advantage of approximations of light propagation to reduce the computational load and time for use in actual quick medical diagnoses. The 2D PDE employed in the algorithm needed a much smaller number of variables in the numerical calculation of the light propagation than 3D RTE. The algorithm that did not required an iterative procedure reduced the calculation time by linearly approximating the relationship between the absorption coefficient of the photon absorber and the absorbed light energy, which generates the PA pressure wave. The approximations will make the image reconstruction algorithm work efficiently in clinical practice.
On the other hand, these approximations used in the algorithm surely cause errors in the reconstructed absorption coefficient. It is known that PDE is not valid in regions near the light source. Additionally, it is not appropriate to use PDE when the absorption coefficient is not sufficiently smaller than the scattering coefficient. Yao and coworkers demonstrated that the RTE-based image reconstruction was more accurate than the PDE-based image reconstruction in a 2D medium when PDE was invalid because of the large absorption coefficient and the anisotropy of the light scattering [30]. However, the effects of the approximations of light propagation on QPAT have not been investigated quantitatively and systematically. For the improvement of the algorithm and for the design and development of a new algorithm, quantitative and systematic investigations of the effects provide indispensable information. Moreover, by the information obtained by the investigation, it may be possible to estimate the true absorption coefficient from the approximated QPAT.
In this study, we calculated the light energies absorbed by the photon absorbers by FEM of 2D/3D PDE and that by 3D Monte Carlo simulation which approximated RTE more accurately than PDE. The differences among the calculations were investigated because the absorbed light energy generating the PA pressure wave was expected to affect the image reconstruction. Then, we discussed the errors caused by the approximations, (1) PDE, (2) 2D model, and (3) linearization, namely, in the 2D PDE-based linearized image reconstruction quantitatively and systematically. Effects of the other factors, such as the optical properties of the background, the size of imaging target, and noise, on the reconstructed image were also investigated. A phantom experiment was carried out to validate the numerical simulations. Additionally, the true absorption coefficient in a 3D phantom was estimated from the 2D reconstructed image on the basis of the prediction by numerical simulations.

Light propagation in biological medium
Nanosecond pulse-laser light to excite the PA pressure wave is scattered and absorbed by tissues. The propagation of the light is described by the following RTE [4,23,24], s Á r/ðr; sÞ þ l a ðrÞ þ l s ðrÞ f g /ðr; sÞ Pðs; s 0 Þ/ðr; s 0 ÞdX 0 þ q 0 ðr; sÞ; ð1Þ where /ðr; sÞ is the radiance at the position r with the direction s. l a ðrÞ and l s ðrÞ are the absorption and scattering coefficients, respectively. Pðs; s 0 Þ is the scattering phase function, which describes the probability of the change in direction s 0 to s. q 0 represents the light source. The light propagation and the heating process to generate the PA pressure wave occur in a much shorter time range than the propagation of the PA pressure wave. Thus, the time dependence of the RTE was ignored. The solution of the RTE can be obtained by the finite difference method or FEM [5][6][7]. The MC simulation statistically approximates the solution of the RTE by tracking the positions and energies of the randomly moving photon packets [31]. In the MC simulation, the Henyey-Greenstein (H-G) phase function is used as Pðs; s 0 Þ. The anisotropy factor g, which is the average of s Á s 0 weighted with Pðs; s 0 Þ, is prescribed to determine the shape of the H-G phase function.
In this study, the following time-independent PDE, which is the approximation of the RTE [4], is used to formulate the reconstruction algorithm, where D ¼ 1=f3ðl 0 s þ l a Þg is the diffusion coefficient calculated with the reduced scattering coefficient l 0 s ¼ ð1 À gÞl s . UðrÞ represents the fluence rate integrating /ðr; sÞ with respect to s. It is known that the PDE is appropriate under the conditions where the medium is sufficiently thick and a long time has passed after the pulse light incidence. The scattering of light is treated as isotropic. The boundary condition, U þ A 2c n Á DrU ¼ 0, is given where n is the vector normal to the surface of the medium and A is the parameter depending on the internal reflection ratio. c is 1=p and 1 / 4 in 2D and 3D calculations, respectively [6]. The PDE was solved by FEM [32] in this study. A was set as unity by assuming that the refractive index matched at the boundary. The light source in arbitrary unit was placed at a depth of 1=l 0 s mm from the surface in the FEM calculation.
The PA source intensity depends on the absorbed energy, which is described as HðrÞ ¼ l a ðrÞUðrÞ, and can be calculated from the solution of Eq. (2). The PA pressure wave is generated by not only diagnostically interesting tissues, which are referred to as ''target'', such as cancer tissues with highly concentrated blood, but also the surrounding normal tissues, which are referred to as ''the background'' in this paper.

PA pressure wave propagation
H is converted into heat. Then thermal expansion generates the PA pressure wave. When the conditions of the thermal and stress confinements [2,3] are met, the PA pressure wave is efficiently generated. The propagation of the PA pressure wave in an acoustically homogeneous medium is described by the following wave equation [2,3]: where v is the speed of the PA pressure wave p,Ĉ, the Grüneisen parameter associated with the PA efficiency, and r, the position. The absorption and conversion to heat of light take a very short time compared with the propagation of the PA pressure wave. Therefore, H(r, t) in Eq. (3) is formulated as Hðr; tÞ ¼ l a UðrÞdðtÞ. 2D FEM was employed to simulate the PA pressure wave in this study.

2D PDE-based linearized image reconstruction algorithm
Let us assume that the illumination and detection are conducted at identical positions by a probe, that the probe illuminates the surface of the measured object, and that the probe is sensitive only to the PA pressure waves generated at the positions directly below the probe. Then the measured PA pressure wave which is linearly related to H below the probe at the k-th position is formulated as where m k 2 R T is the PA pressure waves measured at T discrete points of time by using the probe at the k-th measurement position (k ¼ 1; 2; . . .; K), and H k 2 R M is the absorbed energies at M discrete positions below the k-th measurement position. L k 2 R TÂM represents the contributions of H k to m k . The column vector of L k is obtained by calculating the measured PA pressure wave generated by an unit energy at each of the M positions. Equation (4) is rewritten as follows by the first-order Taylor expansion of H k about the background absorption coefficient l a , where H k is the energy absorbed by the background and Dl a 2 R N is the change in l a at all positions of the medium discretized into N positions. J k 2 R MÂN is the differential coefficients of H k at l a that consists of J k ii ¼ U i þ l a i Á oU i =ol a i and J k ij ¼ l a i Á oU i =ol a j as i 6 ¼ j, where i; j ¼ 1; 2; . . .; N. oU i =ol a j is calculated by FEM and the perturbation method based on Eq. (2) [32]. The reconstruction of Dl a is medically important to find targets, such as the cancer with angiogenesis and the sentinel lymph nodes with a contrast agent, which have larger absorption coefficients than the background. The set of differences between the PA pressure waves measured at two positions is used for the reconstruction to eliminate the PA pressure wave generated by the background and the trend of the baseline of the measured PA pressure wave. The trend is caused by pyroelectric effect in the probe and is commonly included in the PA pressure waves measured at different positions. By subtraction, the PA pressure wave from the background and the trend are removed. We assume that the medium is homogeneous acoustically and optically, which means that L k H k ' L l H l ; ðk 6 ¼ lÞ among the measurement positions. Then, we obtain Opt Rev (2017) 24:705-726 707 The reconstruction of Dl a is carried out by using the following forward equation: where Dm 2 R CT consists of Dm k;l , and the matrix G 2 R CTÂN has G k;l in a row. C ¼ KðK À 1Þ=2 is the total number of the combinations of k and l. Dl a is reconstructed by solving the optimization problem with the Tikhonov regularization [33] formulated as where k is a regularization parameter to adjust the effect of the regularization term on the reconstruction. Then the reconstructed image, c Dl a , is obtained as, where T indicates the transpose of the matrix. The reconstructed b l a was obtained as b l a ¼ c Dl a þ l a , if l a is given [27][28][29].
We assumed that the medium was a 2D square with a 50 mm side which had l 0 s ¼ 0:8 mm À1 and l a ¼ 0:0023 mm À1 [34] as the background. The probe measured the PA pressure waves at 11 positions of x ¼ À10; À8; . . .; 10 mm on the boundary of the medium at y ¼ 0 mm. The PA pressure waves generated from the regions with a width of 2 mm directly below the probe were detected. The 2D FEM calculation was carried out with 10,201 FEM nodes and 20,000 triangle elements. The nodes were placed with 0.5 mm spacing.
The image reconstruction was carried out on pixel basis, and the reconstructed image was displayed with linear interpolation. The medium was discretized into 625 pixels. 25 FEM nodes were contained in each pixel having 2 mm sides. l a and l 0 s were uniformly distributed in the single pixel in the calculation of the matrices G and J. The pixelbased algorithm reduced the matrix size and accelerated the image reconstruction.
The parameter k at the corner of the L-curve is usually selected on the basis of the conventional L-curve method [33]. In this study, k that minimized Dh=De ¼ jhðk u Þ À hðk uþ1 Þj=jeðk u Þ À eðk uþ1 Þj was used, where e is the residual error and h is the regularization function in Eq. (8). k u and k uþ1 were the neighboring numbers in a monotonically increasing numerical sequence of k. In other words, the solution of Eq. (9) with k at the flattest position on the L-curve was chosen as the reconstructed image. This criterion compromised on minimizing e, but reduced the artifact more than the conventional criterion choosing the corner of the L-curve [28].
The procedure of the image reconstruction on pixel basis with the 2D PDE-based linearized algorithm is summarized as follows: Preparation steps 1. Calculate the distribution H in each of M pixels below the k-th illuminating position (k ¼ 1; 2; . . .; K) by solving Eq. (2) with the background l a and l 0 s . 2. Calculate p from each of M pixels by solving Eq. (3) with H normalized to have unit energy for each pixel. Then, construct the matrix L k having p from each pixel in a column. 3. Calculate the matrix J k by the perturbation method based on Eq. (2). 4. Calculate G k;l ¼ L k J k À L l J l ; ðk 6 ¼ lÞ. Then, construct the matrix G having G k;l in a row.

Reconstruction steps
1. Measure m k at K positions. 2. Calculate Dm k;l ¼ m k À m l ; ðk 6 ¼ lÞ. Then construct the vector Dm having Dm k;l in a row.

Display c
Dl a ðkÞ found in step 4.

Effects of the approximations on the absorbed light energy
The light energies absorbed by the targets, H targ , were calculated with the light propagation models to investigate the differences caused by the approximations. The differences in H targ caused by the approximations were expected to cause the errors in the image reconstruction, because the PA pressure wave simulated by the forward calculation depended on H targ . Figure 1 shows the geometrical conditions of the calculations. The backgrounds of both 2D and 3D media had l a ¼ 0:0023 mm À1 and l 0 s ¼ 0:8 mm À1 . We assumed that the illuminating laser light had the wavelength in the nearinfrared (NIR) range in this study. The NIR light is penetrated deeply into the biological media, because the hemoglobin and water that are contained substantially in the living body do not absorb the NIR light considerably, although the main absorber of the NIR light in the living body is hemoglobin. The optical properties of the biological tissues have been studied in the research field of diffuse optical imaging. It was reported that healthy or diseased breast tissues had l a from 0.001 mm À1 (0.01 cm À1 ) to 0.02 mm À1 (0.2 cm À1 ) and l 0 s from 0.5 mm À1 (5 cm À1 ) to 1.0 mm À1 (10 cm À1 ) in the NIR range [10,35]. Diseased tissues tended to have a larger l a than normal tissues. Breast tissue diagnosis is one of the important applications of QPAT and diffuse optical tomography (DOT). Many numerical studies of DOT employed the optical properties in the NIR range [36,37]. The optical phantom that simulated tissues had l a ¼ 0:001 mm À1 and l 0 s ¼ 0:8 mm À1 [38]. The optical properties used mainly in this numerical simulation were selected on the basis of those previous publications.
The center of the target was placed at the depths y ¼ 1, 3, 5, 7 and 9 mm at x ¼ 0. The 2D medium had the squared target with 2 mm sides, that is, the size of the single pixel of the reconstructed image. The 3D medium had the cylindrical target with height and diameter of 2 mm. The target had l a in the range from 0.0023 to 1.0 mm À1 . The light sources were placed above the target.
H targ was calculated by the 2D PDE, 3D PDE and 3D MC simulations with g ¼ 0.66 and 0.9. g ¼ 0.9 was assumed to simulate the biological medium, which scatters the light forward strongly. g ¼ 0.66 simulated the phantom made of an aqueous solution of intralipid. For 2D PDE, the FEM described in the previous section was used. For 3D PDE, 75,626 nodes and 345,600 tetrahedral elements were employed. For the 3D MC simulation, the Monte Carlo for Multi-Layer media code [31] was used, which was modified for the medium to have the target. The medium was discretized with the pixels with 0.5 mm sides. 10 6 photon packets were tracked. The absorbed light energies calculated by 3D PDE and MC simulations were divided by the light source intensities to be comparable.
The effects of the approximations on H targ were examined on the basis of the simulations. Firstly, H targ with PDE was compared with that with MC simulation. Secondly, H targ with 3D PDE was compared with that with 2D PDE.
Thirdly, H targ values with and without linearization were compared.

Effects of the background on the absorbed light energy
Additionally, the effects of the background optical properties on H targ were examined by considering a variety of tissues and wavelengths. H targ at a depth of 5 mm was calculated for l 0 s ¼ 0.2, 0.4, 0.8, 2, and 4 mm À1 with l a ¼ 0.0023 mm À1 . The l a of the target varied as 0.0023, 0.01, 0.1, and 1 mm À1 . H targ was also calculated with l a ¼ 0.02 mm À1 with l 0 s ¼ 0.8 mm À1 to be compared with the cases with l a ¼ 0.0023 mm À1 . The l a of the target took the values from 0.02 to 1 mm À1 . MC simulation with g ¼ 0.9 and 3D PDE were used for the calculations.

Effects of the approximations on the image reconstruction
The image reconstructions from the PA pressure waves simulated with the light propagation models were attempted to investigate the error caused by the approximations used in the 2D PDE-based linearized algorithm. The geometrical conditions were identical to that in the previous section. The absorption coefficient of the target was varied from 0.0023 to 1.0 mm À1 . The simulated distribution of the absorbed energy was used to calculate the PA pressure waves by FEM. The distribution of the absorbed energy on the x-y plane calculated with the 3D PDE was linearly interpolated to be matched with the 2D FEM nodes. The pixels of the MC simulation were overlaid on the 2D FEM nodes. Gaussian noise having the standard deviation (SD) of 1% of the maximum of Dm was added to Dm.
The errors in c Dl a of the target reconstructed with the approximations were examined. Firstly, c Dl a reconstructed from PA pressure waves simulated with 3D PDE and MC simulation were compared. Secondly, c Dl a reconstructed from PA pressure waves simulated with 3D PDE was compared with that with 2D PDE. Thirdly, c Dl a was compared with the linearized value of c Dl a .

Effects of the background and others on the image reconstruction
Firstly, the effects of the background optical properties on the b l a of the target at a depth of 5 mm were examined. The image reconstructions were carried out for the targets with l a ¼ 0.01, 0.1, and 1.0 mm À1 at a depth of 5 mm with the background of l 0 s ¼ 0.2, 0.8, and 2 mm À1 and l a ¼ 0.0023  mm À1 . The image reconstructions were also conducted for l a ¼ 0.02 mm À1 with l 0 s ¼ 0.8 mm À1 . The target had l a from 0.04 to 1.0 mm À1 .
Secondly, the effect of the target size was investigated. The target size does not always fit the pixel size. The image reconstructions using the 2-mm pixel for the 2D target having side lengths of 0.5 and 1 mm at a depth of 5 mm were compared with the image reconstruction with the target fitting the pixel with a side length of 2 mm.
Thirdly, to investigate the effect of noise, the image reconstructions were carried out with Gaussian noise having SDs of 10 and 40% of the maximum of Dm to be compared with the image reconstruction with the SD of 1%. The 2D target that had 2-mm sides and l a ¼ 0.01 mm À1 was placed at a depth of 5 mm.
The numerical simulations and the image reconstructions were performed using computer with a 3.3 GHz CPU and 128 GB of RAM. The FEM calculations and the image reconstructions were performed using Matlab (The Math-Works Inc., Natick, MA, USA).

PAT system
The phantom experiments were carried out to validate the numerical results in the previous section, and to explore the possibility of estimating the true l a of the target in the 3D medium from the 2D PDE-based linearized image reconstruction. Figure 2 shows the PAT system and the structures of phantoms used in the phantom experiment [29]. A tunable Ti:sapphire laser pumped by the second harmonic of a Q-switched Nd:YAG laser (LT-2211, LS-2134, Lotis Tii, Minsk, Belarus) was used to produce excitation light pulses with a width of about 10 ns and a repetition frequency of 15 Hz at a wavelength of 756 nm. The light pulse was coupled into a multimode optical fiber (M41L02, Thorlabs, Newton, NJ, USA) with a core diameter of 0.6 mm. The optical fiber was introduced into a cylindrical probe that had a ring-shaped piezoelectric film P(VDF-TrFE) (KF piezo-film, Kureha Corp., Tokyo, Japan) on the detection surface with a diameter of 3 mm [39]. The edge of the optical fiber was located at the center of the detecting surface. Therefore, the illumination and detection are conducted at identical positions. The energy of the illuminating light was 4 mJ/pulse. The signal was amplified by the amplifier (SA-220F5, NF Corporation, Yokohama, Japan) and recorded with a digital oscilloscope (DSO7104A, Agilent Technologies, Santa Clara, CA, USA). The sampling rate was 20 MHz. The measured PA pressure waves were averaged over 256 light pulses online. The PAT images were reconstructed by a personal computer from the measured PA pressure wave.

Phantom
A liquid phantom of an aqueous solution of the intralipid and indocyanine green (ICG) was measured. The optical properties of the phantom were adjusted to l 0 s ¼ 0:8 mm À1 and l a ¼ 0:0023 mm À1 [40]. The acoustic coupling polymer gel (SONAGEL, TAKIRON, Tokyo, Japan) with a thickness of 10 mm was used as a spacer between the probe and the liquid phantom so that the PA signal measured in the period from 0 s to 6.7 ls was eliminated and was not used for the image reconstruction. The PA pressure wave measured in the period was severely contaminated by noise owing to the illumination adjacent to the detecting surface of the probe. We tried image reconstruction for 2 cases  with (i) a single target and (ii) multiple targets. In case (i), the target, which was diluted intralipid and ICG contained in a tube with an inner diameter of 1 mm, had the optical properties of l 0 s ¼ 0:8 mm À1 and l a ¼ 0.6, 1.1, 1.7, or 2.3 mm À1 . Therefore, the pixel-sized region including the tube had l a ¼ 0.12 0.22, 0.33 and 0.45 mm À1 on average, because of the ratio of the cross-sectional area of the tube, 0.79 mm 2 , to the area of the pixel, 4.0 mm 2 . The circular cross section of the target was located at y ¼ 5, 7, or 9 mm and x ¼ 0 mm. In case (ii), 3 targets were placed at ðx ½mm; y ½mmÞ ¼ ðÀ6; 5Þ; ð0; 7Þ and (6,9). Each region with the target tube had l 0 s ¼ 0:8 cm À1 and l a ¼ 0:22 mm À1 .

Estimation of l a in 3D medium
The c Dl a of the target reconstructed in case (i) was compared with the c Dl a reconstructed from the PA pressure wave, which was simulated with MC simulation (g ¼

0.66). For the comparison, c
Dl a with MC simulation in the range of the true l a from 0.12 to 0.45 mm À1 was fitted by a quadratic curve as a function of the true Dl a and the depth y of the target as c Dl a ¼ f ðDl a ; yÞ. It was expected that f was able to predict the c Dl a of the target in the 3D phantom reconstructed by the 2D PDE-based linearized algorithm.
In case (ii), it was expected that the 2D PDE-based linearized image reconstruction caused the errors in the reconstruction of the Dl a of the targets in the 3D phantom. On the basis of the quadratic function f calculated in case (i), the true Dl a values of the targets were estimated from c Dl a as Dl a ¼ f À1 ð c Dl a ; yÞ.
3 Results and discussion 3.1 Numerical simulations Figure 3 shows H targ ðmodelÞ, where ''model'' in the parentheses indicates the light propagation models, namely, 2D PDE, 3D PDE, and MC simulation (g ¼ 0.9 and 0.66) and their categories, which calculated H targ . The parameters in the simulation, namely, g and the depth of the target, y, are also noted in the parentheses hereafter if necessary. H targ was linearly related to the PA pressure wave, and the differences in H targ could cause the errors in the image reconstruction. Hence, on the basis of the results shown in Fig. 3, the differences in H targ calculated with the approximations were investigated.

Effects of the approximations on the absorbed light energy
(1) Effect of PDE Firstly, the difference between MC simulation and 3D PDE was investigated. Figure 4 shows the absorbed light energies in the cylindrical regions with a diameter and height of 2 mm with l a ¼ 0:0023 mm À1 at depths of y ¼ 1, 3, 5, 7 and 9 mm in the homogeneous medium. PDE caused error in the region at y ¼ 1 mm where the illuminating position was closely located. H targ (3D PDE, 1 mm) was 67% of H targ (MC, 1 mm). PDE approximates the light propagation well when the light is scattered several times, and scattering randomizes the directions of light propagation [41]. Therefore, PDE did not approximate the MC simulation in the position close to the illuminating position where the laser light conserved the directionality. This simulation indicated that PDE was effective at the positions deeper than 3 mm at least. When the depth was larger than 3 mm, 3D PDE and MC simulations agreed well. The difference between MC simulations with g ¼ 0.66 and 0.9 was small. On the other hand, when the cylindrical target with l a ¼ 0:01 mm À1 existed at y ¼ 5 mm, H targ (3D PDE, 5 mm) was larger than H targ (MC, 5 mm), although the absorbed light energies with 3D PDE coincided with those with MC simulations at the depths of y ¼ 3, 7 and 9 mm as shown in Fig. 4.
H targ ð3D PDEÞ was about 1.4 times as large as H targ ðMCÞ when the target with l a of 0.01 mm À1 was at y ¼ 5 mm, as shown in Fig. 5a. The difference in H targ tended to increase with l a until H targ stopped increasing with large l a . The ratio of H targ ð3D PDEÞ to H targ ðMCÞ is plotted in Fig. 5b. The ratio increased until l a reached approximately 0.1 mm À1 for the depths from 3 mm to 9 mm. It was observed that the difference between PDE and MC simulations is enlarged when l a becomes comparable to l 0 s in Fig. 5b. The ratio depended on the depth of the target. The deeper the target existed, the larger the ratio. The maximum ratio in the simulation was 1.7 with the target having l a ¼ 0.1 mm À1 at the depth of 9 mm. The error of PDE increased as the depth of the target and l a increased.
As observed in Figs. 3 and 5a, especially when the target had l a ! 0:1 mm À1 , the increasing rate of H targ became small. That means that the PA pressure wave did not change considerably with the change in l a larger than 0.1 mm À1 . In this circumstance, it is difficult to precisely identify Dl a from the PA pressure wave, which was contaminated by noise. As a result, the c Dl a of the target can be varied largely by small noise. When the target has a large l a and the signal-to-noise ratio of the PA pressure wave is small, a large error in c Dl a can occur. This reconstruction Opt Rev (2017) 24:705-726 711 error is referred to in the discussion of the results of the phantom experiment. The result shown in Fig. 5b indicated that in the target with a large l a , the light propagation was different from diffusion, and it became difficult for PDE to approximate the light propagation precisely. The error caused by PDE with a homogeneous medium having large l a was discussed previously [41]. It was reported that U calculated with PDE became larger than that with MC simulation as l a increased. The difference between PDE and RTE was demonstrated in the 2D medium in the literature [30]. It is said that the photon flux calculated with 2D PDE was smaller than that with 2D RTE in the target region with a large absorption coefficient. These results indicated that H targ (3D PDE) became larger than H targ (MC) and agreed with the numerical simulation of this study.
The differences between H targ (3D PDE) and H targ (MC) shown in Figs. 4 and 5 originated from the equations that were used to calculate the light propagation. MC simulation can approximate RTE better than PDE. Therefore, the differences can be explained by observing the differences between PDE and RTE. Although Green's functions of RTE have been studied [42,43], the use of the equation is complicated. To understand how the differences between H targ (3D PDE) and H targ (MC) occurred, we tried to compare 1D light propagation models with RTE and PDE. A one-dimensional transport theory can be used to develop a three-dimensional theory [43]. They are deeply related to each other. We assumed that light propagated into a 1D target with unit net flux as shown in Fig. 6a. The fluence rate in the target, U RTE , calculated with 1D RTE was compared with U PDE with 1D PDE (see Appendix 1). Because the absolute value of the power of the exponential term of U RTE in Eq. (16) is smaller than that of U PDE in Eq. (21), U PDE decreases more rapidly than U RTE . As a result, H targ (PDE) is different from H targ (RTE). Figure 6b plots H targ (PDE)/H targ (RTE) as a function of ðl a =l 0 s Þ. It was demonstrated that H targ (PDE) ! H targ (RTE). When ðl a =l 0 s Þ was small, PDE approximated RTE. The difference between H targ (PDE) and H targ (RTE) became larger as ðl a =l 0 s Þ increased. It is well known that PDE cannot approximate RTE when ðl a =l 0 s Þ is large. The result demonstrated this fact evidently. When ðl a =l 0 s Þ [ 1, most of the light was absorbed within the target. Therefore, the difference became small. A similar mechanism can explain the reason why H targ (3D PDE) [ H targ (MC) at the depths from 3 to 9 mm, as shown in Figs. 4 and 5. In the 3D case, the light propagated through paths around the target. Thus, the light energy absorbed by the tissues surrounding the target was not affected by the target so much when the target had a small l a .
(2) Effect of 2D model Secondly, the difference between 2D and 3D light propagation models was investigated. To compare the dependence of H targ on the depth y, the ratio, RðmodelÞ ¼ H targ ðmodel; yÞ=H targ ðmodel; 3 mmÞ, was calculated with l a ¼ 0.01 mm À1 of the target. RðmodelÞ was plotted for the depths y ¼ 1, 3, 5, 7 and 9 mm in Fig. 7a. In both 2D and 3D models, the absorbed light energy decreased as the depth increased. However, the decreasing rates were different. Rð2D PDEÞ was approximately 0. more rapidly than H targ ð2D PDEÞ, because the light propagated in the direction orthogonal to the x-y plane, and the light that reached the target decreased in 3D models more than in the 2D model. Rð2D PDEÞ=Rð3D PDEÞ is shown in Fig. 7b, which indicated that with 2D PDE, the calculated ratio of H targ was approximately 1.6, 2.0 and 2.3 times as large as that with 3D PDE at the depths of 5, 7 and 9 mm, respectively. The ratio did not depend so much on l a , although the ratio became slightly small when l a was larger than 0.1 mm À1 .
To demonstrate that the ratio Rð2D PDEÞ=Rð3D PDEÞ shown in Fig. 7 is caused by the difference between 2D and 3D PDEs, their solutions were compared (see Appendix 2). U 2D ðrÞ and U 3D ðrÞ were calculated with 2D and 3D PDEs, respectively. The media had l a ¼ 0:0023 mm À1 and l 0 s ¼ 0:8 mm À1 . Then, the ratio of U 2D ðrÞ to U 3D ðrÞ was calculated for r ¼1, 3, 5, 7, and 9 mm (Fig. 8). It was found that the ratio shown in Fig. 8 coincided with the ratio in Fig. 7, although the error between Figs. 7 and 8 caused by the pixel size was observed. It was demonstrated that the ratio depended on the difference between 2D PDE and 3D PDE. U and H were linearly increased with the l a of the target to some extent while l a 0:1 mm À1 as shown in Fig. 3. Therefore, the ratio took an almost constant value in the range of l a 0:1 mm À1 . The ratio changed nonlinearly with l a [ 0:1 mm À1 in Fig. 7b.
(3) Effect of linearization Thirdly, the effect of linearizing the relationship between l a and H targ was investigated. Figure 9a shows H targ (2D PDE, 5 mm) and the linearized value as a function of l a . The absorbed light energy did not increase linearly to l a as shown in Fig. 9a. The error between the actual and the linearized absorbed light energy became large as l a became large. The relative error between the actual and linearized H targ ðmodelÞ was plotted in Fig. 9b. The relative errors at the depths of 1, 3, 5, 7 and 9 mm were averaged. The relative error was about 0.01 and smaller when the target had l a smaller than 0.01 mm À1 . The relative error became larger than 0.1 when l a ¼ 0:1 mm À1 . The linearization expanded the error in the calculation of H targ as l a increased. Figure 10a plots H targ (MC, 5 mm) with l 0 s ¼ 0.2, 0.4, 0.8, 2, and 4 mm À1 . As l 0 s increased, H targ (MC, 5 mm) decreased, because the light that reached the target decreased owing to scattering. Similarly to the result shown in Fig. 9, H targ began to stop increasing when the l a of the target became larger than 0.1 mm À1 . The tendency became more significant with the increase in l 0 s . The nonlinear character of H targ strongly appeared when l 0 s was larger. The ratio H targ (3D PDE)=H targ (MC) is shown in Fig. 10b. 3D PDE approximated MC simulation well while 0:4 l 0 s 2, although the ratio was larger than unity because H targ (3D PDE) [ H targ (MC), as demonstrated previously. Since PDE underestimated the light that reached the target owing to strong scattering, the ratio became smaller when l 0 s ¼ 4 mm À1 . It was reported that some tissues have l 0 s ! 0:1 mm À1 especially for l 0 s . The error between 3D PDE and MC simulation became larger especially when l 0 s ¼ 0:2 mm À1 and l a ¼ 0:1 and 1 mm À1 of the target. This indicated that PDE could not approximate RTE when l 0 s was small and l a was large. As shown in Fig. 10b, the error was large when the light was propagated in the medium with a small l 0 s and absorbed suddenly by the region with a large l a . Therefore, it is expected that the error in the PDE calculation is significant when the light is absorbed after propagating in the void-like region which hardly scatters and absorbs light. The error can be enlarged as l 0 s becomes smaller. Figure 11a compares H targ (3D PDE) and H targ (MC) with l a ¼ 0.02 mm À1 . The relationship between H targ (3D PDE) and H targ (MC) was similar to that shown in Fig. 5a. H targ (3D PDE) became larger than H targ (MC) as l a of the target increased. The ratio H targ (3D PDE)=H targ (MC) with l a ¼ 0.02 mm À1 was compared with that with l a ¼ 0.0023 mm À1 in Fig. 11b. The ratio increased while the target had l a 0:2 mm À1 . Then the ratio decreased when the target had l a ! 0:2 mm À1 . The ratio with l a ¼ 0.02 mm À1 changed similarly to that with l a ¼ 0.0023 mm À1 . l a depends on the concentration of the chromophores absorbing light such as hemoglobin, lipid and water [44]. Therefore, l a can become larger than 0.0023 mm À1 by depending on the tissue conditions and wavelength used in the PA measurement. While l a is around 0.01 mm À1 , PDE can approximate RTE with the errors shown in Fig. 11. However, it is difficult for PDE to approximate RTE when ð l a = l 0 s Þ is considerably larger. The numerical simulations demonstrated that the approximations of the light propagation caused the errors in the calculation of the absorbed light energy. Light propagation in the scattering media such as biological tissues obeys 3D RTE and the relationship between l a and H targ is nonlinear by nature. Therefore, it is impossible to examine the effects of the approximations using PDE, 2D model and linearization with phantom experiments. The effects of the approximations can occur only in the forward calculation used for the image reconstruction algorithm, and can be examined only by the numerical simulations. The approximations affect the image reconstruction through the forward calculation by causing the mismatch between the calculations and the real phenomena. The PA pressure wave depends on H targ . Therefore, the mismatches in H targ between the forward calculation with the approximations and the real experiments affect the simulation of  Figure 12 shows the c Dl a ðmodelÞ of the target reconstructed by the 2D PDE-based linearized algorithm from the PA pressure wave calculated with various Dl a values and depths. The ''model'' in the parentheses indicates the light propagation model, 2D PDE, 3D PDE, MC simulation, and their categories which were used to calculate the PA pressure wave. The parameters in the simulation, namely, g and the depth of the target, y, are also noted in the parentheses if necessary hereafter. Figure 13 shows the images reconstructed from the PA pressure waves simulated with 2D PDE, 3D PDE, and MC simulations (g ¼ 0:9 and 0.66). The target with l a ¼ 0.01 mm À1 was placed at the depth of 5 mm. The maximum of c Dl a was reconstructed at the target position correctly in each case of the numerical simulations of the image reconstruction.

Effects of the approximations on the image reconstruction
On the basis of the reconstructed c Dl a (model), the errors caused by the 2D PDE-based linearized image reconstruction are discussed next.
(0) Performance of the algorithm The 2D PDE-based linearized algorithm reconstructed c Dl a ð2D PDEÞ as approximately 72% of the true Dl a when the target had Dl a 0:0077 mm À1 (l a 0:01 mm À1 ). The percentage decreased to 53% when the target had Dl a in the range from 0.0177 to 0.0977 (0:02 l a 0:1 mm À1 ). The 2D PDE-based linearized algorithm reconstructed the image stably from the noisy PA pressure wave. This character of the algorithm is important for obtaining the solution of the inverse problem. Generally, measurement noise and the mismatches between the forward calculation and the measurement often expand the errors in image reconstruction solving ill-posed inverse problems. When the regularization method did not function sufficiently, the artifacts caused by the measurement noise and modeling error appeared in the image reconstructed by the 2D PDEbased linearized algorithm.
The image reconstruction by the 2D PDE-based linearized algorithm took 2.8 s. The 2D PDE-based linearized algorithm reconstructed the image fast enough to be used in clinical practice. That was achieved by the approximations and the pixel based algorithm. In an advanced gradientbased 2D image reconstruction with RTE, it took approximately 16 h to obtain the reconstructed image with 10,000 computational meshes [23].  According to the simulation results, it is better to use PDE in the image reconstruction to reconstruct the target at the depth larger than 3 mm. When the imaging object is not sufficiently thick, it is not appropriate to use PDE. larger than c Dl a ðMCÞ at the depth of 5 mm. This tendency in Fig. 14a was similar to that in Fig. 5a. Figure 14b shows the ratio of c Dl a ð3D PDEÞ to c Dl a (MC, g ¼ 0.66) as a function of Dl a of the target. The ratio indicated that c Dl a ð3D PDEÞ was 1.6 times as large as c Dl a (MC, g ¼ 0.66) on average at depths larger than 3 mm. The ratio approximately coincided with H targ (3D PDE)/H targ (MC) shown in Fig. 5b. Thus, it can be said that c Dl a was affected directly by the error in H targ with the PDE calculation shown in Fig. 5b. From the comparison of Figs. 5b and 14b, it was observed that c Dl a by the 2D PDE-based linearized algorithm was changed to a similar extent with the change in H targ related to the amplitude of the PA pressure wave. smaller than c Dl a (3D PDE) at depths from 3 to 9 mm, because the MC simulation calculated the amplitude of the PA pressure wave smaller than 3D PDE. Similarly, it can be anticipated that a 3D PDE-based algorithm reconstructs c Dl a (MC) smaller than the true Dl a , because 3D PDE expects a larger amplitude of the PA pressure wave than MC simulations for identical l a of the target. This error can be evaluated from the result shown in Fig. 5b. In real image reconstruction, which was approximated by MC simulation more accurately than by 3D PDE, the same error will occur in a PDE-based image reconstruction.
(2) Effect of 2D model Secondly, the difference between  plots R r (model). R r (3D models) decreased as the depth increased. In quantitative image reconstruction, the c Dl a of the targets having an identical Dl a should be reconstructed as an identical value, regardless of the depths of the targets.

However, c
Dl a (3D models) varied with the depth. R r ð3D modelsÞ was 0.55, 0.44 and 0.35 at depths of 5, 7 and 9 mm, respectively, while R r (2D PDE) was almost unity. Figure 15b shows the ratio of R r (2D PDE) to R r (3D PDE). The ratio was approximately 1.5, 2.0 and 2.5 for y ¼ 5, 7 and 9 mm, respectively, although some fluctuations that could be caused by noise added to the measurement and the selection of k was observed. The ratio coincided with Rð2D PDEÞ=Rð3D PDEÞ shown in Fig. 7b. This means that the error shown in Fig. 7b was preserved after the image reconstruction. The difference in the dependence Dl a (3D models) reconstructed from the PA pressure wave generated by the target at 5 mm in 3D medium. The PA pressure waves were simulated with 3D PDE (circle symbol) and MC simulations with g ¼ 0:9 (square symbol) and 0.6 (triangle symbol), and b ratio of the reconstructed values from the PA pressure wave simulated with 3D PDE to that with MC simulation (g ¼ 0:66). The ratio is plotted for the depths of the targets of 1 mm (circle symbol), 3 mm (square symbol), 5 mm (triangle symbol), 7 mm (times symbol) and 9 mm (cross symbol) of the depth between 2D and 3D light propagations cannot be overcome by the 2D PDE-based linearized image reconstruction. In the image reconstruction of a real 3D medium, it was expected that the Dl a of the target at a deep position was underestimated by the 2D image reconstruction.
As shown in Fig. 15a, R r ð2D PDE; 1 mmÞ was 0.7, while R r was almost unity at the other depths. The error in c Dl a (2D PDE, 1 mm) was relatively larger than the errors at y ! 3 mm. Recall that L k in the formulation of the forward equation Eq. (7) had the PA pressure waves generated by the unit absorbed light energy. L k was determined prior to the image reconstruction. Meanwhile, in the pixel at y ¼ 1 mm, which was close to the illuminating position, the shape of the PA pressure wave generated in the pixel was changed by the change in the distribution of the absorbed light energy in the pixel with Dl a . That caused the mismatch between the PA pressure wave and the forward calculation at y ¼ 1 mm. This mismatch was relatively small in the deeper pixel, because the distribution of the absorbed light energy was smoothed by the diffusive propagation of the light. Therefore, the error in c Dl a (2D PDE, 1 mm) became larger. The mismatch aggravated the overestimation of c Dl a (MC, 1 mm), in addition to the error caused by the use of PDE. The reconstruction pixel should be sufficiently small to improve c Dl a . Additionally, the approximations such as the d-Eddington method to reduce the error in H near the illuminating position may be effective [21,45].
(3) Effect of linearization Thirdly, the effect of the linearization in the image reconstruction was investigated.  Figure 16a is similar to Fig. 9a. Therefore, it was clear that the linearization directly caused the error in the image reconstruction. The true value of H targ was smaller than the linearized value especially when l a ! 0.1 mm À1 as shown in Fig. 9a. Then the measured amplitude of the PA pressure wave also became smaller than the amplitude expected by the linearized algorithm.
Since the linearized algorithm reconstructed c Dl a in response to the amplitude of the PA pressure wave, c Dl a was smaller than the true value. The error between the true Dl a and c Dl a by the linearized algorithm expanded as the Dl a of the target increased. The relative error between c Dl a (model) and the linearized value are plotted in Fig. 16b. The relative errors at y ¼ 1, 3, 5, 7 and 9 mm were averaged. The relative error became larger than 0.1 when the Dl a of the target became larger than 0.1 mm À1 , similarly to the result in Fig. 9b. The error caused by the linearization was also observed in c Dl a (MC). Therefore, the error by the linearized algorithm was expected in the image reconstruction for the real medium.

Effects of the background and others on the image reconstruction
(1) Effect of background Figure 17a shows the c Dl a (MC) of the target at a depth of 5 mm for l 0 s ¼ 0.2, 0.8 and 2 mm À1 . c Dl a with l 0 s ¼ 0.2 mm À1 was larger than that with l 0 s ¼ 0.8 and 2 mm À1 . This coincides with Fig. 10a that plots H  small l 0 s . On the other hand, the difference between the results of l 0 s ¼ 0.8 and 2 mm À1 was relatively small. It was observed that the difference became larger as the l a of the target increased. Because the nonlinearity of H targ was enhanced by increases in the l 0 s and l a of the target as shown in Fig. 10a, c Dl a with l 0 s ¼ 2 mm À1 became smaller than that with l 0 s ¼ 0.8 mm À1 with the increase in the l a of the target. The ratio of c Dl a (3D PDE) to c Dl a (MC) is plotted in Fig. 17b. Similarly to H targ (3D PDE)/H targ (MC) shown in Fig. 10b, c Dl a (3D PDE) was much larger than c Dl a (MC) with l 0 s ¼ 0.2 mm À1 and the l a of the target of 1.0 mm À1 . PDE could not approximate RTE when ðl a =l 0 s Þ became large. The 2D PDE-based image reconstruction was affected by the error between PDE and MC simulation.
c Dl a (3D PDE) and c Dl a (MC) with l a ¼ 0.02 mm À1 were plotted in Fig. 18a. The result was similar to that with l a ¼ 0.0023 mm À1 shown in Fig. 14a. c Dl a (3D PDE) was larger than c Dl a (MC), because the PA pressure wave generated by H targ (3D PDE) was larger than that by H targ (MC).

The ratio c
Dl a (3D PDE)= c Dl a (MC) was plotted for l a ¼ 0.0023 and 0.02 mm À1 . Although the ratio with l a ¼ 0.02 mm À1 was smaller than that with l a ¼ 0.0023 mm À1 , the ratio was approximately 1.5, and the corresponding tendency was observed in both cases. Some biological tissues can have l a around 0.02 mm À1 in a certain wavelength [44]. The effects of the approximations with l a ¼ 0.02 mm À1 were similar to that with l a ¼ 0.0023 mm À1 discussed in Sect. 3.1.3.
(2) Effect of target size c Dl a (2D PDE) in a 2-mm pixel at a depth of 5 mm was varied with the true target size as shown in Fig. 19. This is because the absorbed light energy in the pixel was varied by the target size. As the target size decreased, c Dl a (2D PDE) became smaller. The ratio plotted in Fig. 19 coincided with the ratio of the side length of the target to that of the target with 2-mm side. Therefore, it seemed that c Dl a (2D PDE) depended on the side length of the target. The amplitude of the PA pressure wave from the target can predominantly depend on U at the edge of the target from which the light enters the target, because U is the largest there. Therefore, the amplitude of the PA pressure wave depended on the side length of the target in this 2D case. Because the side length of the target was regarded as 2 mm of the pixel in the reconstruction, c Dl a (2D PDE) with the 1-mm target was half as large as c Dl a (2D PDE) was normalized to have a unit peak value for each SD. The width of the rise and decay of c Dl a indicating the existence of the target became larger as SD increased. The full widths of half maximum with SDs of 10 and 40% measured in Fig. 20 were 1.2 and 1.4 times as large as that with SD of 1%, respectively. When the noise was large (SD was 40%), the spatially oscillating artifact was observed in the profile. The peak values of c Dl a (2D PDE) were 0.006, 0.0041, and 0.0038 mm À1 with SDs of 1, 10 and 40%, respectively, when the target had l a ¼ 0.01 mm À1 . In the case of the target with l a ¼ 0.1 mm À1 , the peak values changed to 0.053, 0.036, and 0.034 mm À1 with SDs of 1, 10 and 40%, respectively. The ratio of the peak value with l a of 0.01 mm À1 to that with l a of 0.1 mm À1 was approximately 8.8 constantly for each SD. The noise decreased the spatial resolution and the peak value of c Dl a . When the measured PA pressure wave had a small signalto-noise ratio, a large k was selected to suppress the effect of noise on the reconstructed image. As a result, the distribution of c Dl a is excessively smoothed. Spatial resolution can be improved by reducing noise. The pixel size should be adjusted, depending on the required spatial resolution. The measurement system using focusing ultrasound transducers can improve the spatial resolution in PA imaging generally.
The numerical simulation of the image reconstruction demonstrated that the 2D PDE-based linearized algorithm was able to reconstruct the target in the 3D medium. The errors caused by the approximations in the calculation of H targ directly affected the image reconstruction.

Phantom experiments
We attempted the conventional phantom experiment of the image reconstruction in a similar manner to the previous publications to validate the numerical simulations in this

Moreover, c
Dl a (MC) obtained from 2D image reconstruction is used to predict the reconstructed values of the target in the 3D phantom experiment in this section. Figure 21 shows the distribution of c Dl a reconstructed from the measurement data obtained by using the phantom in case (i) with the target regions having l a ¼ 0.12 mm À1 at depths of 5, 7 and 9 mm. The maximum value of c Dl a , which indicated the existence of the target, was reconstructed at the correct position. The 2D PDE-based linearized algorithm can reconstruct the target region in the 3D medium from the real measurement data successfully. As demonstrated in Sect. 3.1.4 (3), noise decreased the spatial resolution. Similarly to the numerical simulation, some spatially oscillating artifacts were observed around the reconstructed target region. The artifact often appears in the solutions of the inverse problems employing truncated singular value decomposition and Tikhonov regularization techniques. The c Dl a of the target region in case (i) was compared with the predicted values by the use of f formulated by the simulation of the image reconstruction from the PA pressure wave calculated with MC simulation (g ¼ 0.66) in Fig. 22. When the target region with l a ¼ 0:45 mm À1 existed at 9 mm, the error between c Dl a and the predicted value was large. This was caused by the large l a of the target and the small signal-to-noise ratio of the PA pressure wave. l a ¼ 0:45 mm À1 of the target region was in the range where the increasing rate of H targ became small as shown in Fig. 3. Additionally, in this phantom experiment, the PA pressure wave from the target at 9 mm had a small signal-to-noise ratio. In this condition, as mentioned in the paragraph on the effect of PDE on H targ , the large error in c Dl a could be caused by the small changes in the PA pressure wave due to measurement noise and error.

Case (i) with single target
Although some mismatches in c Dl a due to noise were found, The c Dl a of the target region of the phantom agreed well with the values predicted using f. This means that the numerical simulation succeeded in the prediction of the values reconstructed by the 2D PDE-based linearized algorithm in the real image reconstruction with a 3D medium. Thus, it seemed possible to obtain Dl a from c Dl a by the use of f, which related c Dl a to the true Dl a .  128 predicted using f. The c Dl a of the target at a depth of 7 mm was 0.12, which was slightly larger than the value of 0.10 predicted using f. Measurement noise and error in the experimental condition caused the mismatch. As expected, the c Dl a of the target region at a depth of 9 mm was the smallest among the three targets. However, c Dl a of 0.049 was much smaller than the predicted value of 0.081. Figure 23b shows the c Dl a and Dl a of the target regions. The Dl a value of the target regions at the depths of 5 and 7 mm were successfully estimated as 0.215 and 0.265 mm À1 , respectively. The errors in Dl a at the depths of 5 and 7 mm were 3 and 20% of the true Dl a , respectively. When a real cancer lesion with high concentrations of hemoglobin or contrast agents had l a larger than the background normal tissues as reported in the literature [9][10][11][12][13][14][15], the error may not disturb the diagnosis. The literature showed that the l a of cancer tissue was 2 -4 times as large as the l a of normal tissue. In such cases, the errors in the estimations at the depths of 5 and 7 mm seemed not critical. The Dl a of the target region, which was 0.11 mm À1 , was much smaller than the true value at the depth of 9 mm, because of the large mismatch between c Dl a and the predicted values. Besides the condition with a large l a of the target and a small signal-to-noise ratio, the regularization caused the large error in c Dl a at the deep position in the case (ii).
When the multiple targets existed in the medium, c Dl a of the target located at a deeper position was underestimated by the regularization technique to smoothen the distribution of c Dl a and to reduce the amount of artifacts owing to the measurement noise. Tikhonov regularization techniques minimizing the norm of the distribution of c Dl a tended to suppress the c Dl a of the target in deeper position, the contribution of which to the measured PA pressure wave was small, although the effect was important to solve the inverse problem stably.

Effects of the experimental conditions
There were some experimental conditions which can cause the mismatches between the numerical simulations and the phantom experiment. The difference in the shape of the target between the MC simulation and the phantom experiment must cause the mismatch in the comparison shown in Fig. 13. However, the numerical simulation agreed well with the phantom experiment in this study. The limited directional angle of the probe measuring the PA pressure wave could reduce the error. The region from which the probe detected the PA pressure wave was so limited that the measured PA pressure wave was not affected by the difference. Additionally, the algorithm based on the 2-mm pixel was not sensitive to the small structure of a thin tube. The effect of smoothing the absorbed light energy distribution by the diffusive light propagation might also reduce the error caused by the difference in the shape of the target.
Although in this study we focused on the effects of the approximations in the light propagation models, it was expected that the modeling of the propagation of the PA pressure waves affected the image reconstruction. The phantom used in this study was a 3D medium, so that the 2D FEM calculation of the PA pressure wave must cause the error in the image reconstruction. However, the phantom had a rectangular parallel-piped shape and the tubes were horizontal to the surface of phantom. This 2D-like geometry of the phantom might alleviate the error in the 2D FEM calculation of the PA pressure wave. The phantom experiments demonstrated that the 2D PDE-based linearized algorithm can image the target in a real 3D medium, and validated the numerical simulation reconstructing c Dl a with the MC simulations in the previous section. Additionally, the method of obtaining the Dl a of the target in the 3D medium based on the numerical simulation of the image reconstruction functioned in this phantom experiment, although some limitations became apparent.

Conclusions
The effects of the approximations of the light propagation on the calculation of the absorbed light energy and the image reconstruction for quantitative photoacoustic image reconstruction were reported in this paper.
It was numerically demonstrated that the approximations such as the use of PDE, 2D model and linearization caused some errors in the calculation of the absorbed light energy. The errors directly affected the image reconstruction taking full advantage of the approximations. The error caused by the use of PDE became larger near the illuminating position. The 2D image reconstruction caused the depth dependence of the reconstructed absorption coefficient in a 3D medium. The use of linearization underestimated the absorption coefficient when the absorption coefficient of the target became large. In the phantom experiment, it was validated that the numerical simulation of the 2D PDE-based linearized image reconstruction can predict the reconstructed value in a 3D medium. Then the true value of the absorption coefficient was estimated on basis of the numerical simulation of the image reconstruction. Although some limitations of the estimation became apparent, the estimation was partly successful.
Precise modeling and measurement of the photoacoustic phenomenon and a certain nonlinear image reconstruction are needed for precise quantitative diagnoses. However, the lack of promptness and simplicity in the PA measurement and the image reconstruction may impair the advantages of optical and ultrasound techniques. The approximations should be appropriately utilized while considering the precision required in medical diagnoses.