Translator Disclaimer
29 January 2018 Deep-tissue temperature mapping by multi-illumination photoacoustic tomography aided by a diffusion optical model: a numerical study
Author Affiliations +
Temperature mapping during thermotherapy can help precisely control the heating process, both temporally and spatially, to efficiently kill the tumor cells and prevent the healthy tissues from heating damage. Photoacoustic tomography (PAT) has been used for noninvasive temperature mapping with high sensitivity, based on the linear correlation between the tissue’s Grüneisen parameter and temperature. However, limited by the tissue’s unknown optical properties and thus the optical fluence at depths beyond the optical diffusion limit, the reported PAT thermometry usually takes a ratiometric measurement at different temperatures and thus cannot provide absolute measurements. Moreover, ratiometric measurement over time at different temperatures has to assume that the tissue’s optical properties do not change with temperatures, which is usually not valid due to the temperature-induced hemodynamic changes. We propose an optical-diffusion-model-enhanced PAT temperature mapping that can obtain the absolute temperature distribution in deep tissue, without the need of multiple measurements at different temperatures. Based on the initial acoustic pressure reconstructed from multi-illumination photoacoustic signals, both the local optical fluence and the optical parameters including absorption and scattering coefficients are first estimated by the optical-diffusion model, then the temperature distribution is obtained from the reconstructed Grüneisen parameters. We have developed a mathematic model for the multi-illumination PAT of absolute temperatures, and our two-dimensional numerical simulations have shown the feasibility of this new method. The proposed absolute temperature mapping method may set the technical foundation for better temperature control in deep tissue in thermotherapy.



Monitoring the temperature distribution in deep tissues during thermotherapy of cancers is critically important. With the temperature map, the heating process can be precisely controlled to efficiently kill the tumor cells and minimize collateral damage to healthy tissues. Currently, magnetic resonance thermometry is the major clinical tool for temperature mapping in thermotherapy, especially for the high-intensity focused ultrasound treatment. However, magnetic resonance thermography is an expensive technology and suffers from low resolution and low speed. For noninvasive real-time temperature mapping, photoacoustic (PA) tomography (PAT) provides an alternative solution due to its inherent sensitivity to temperature.13 In PAT, the initial acoustic pressure generated by the pulsed optical excitation is proportional to the Grüneisen parameter, the optical absorption coefficient, and the optical fluence.4 The Grüneisen parameter Γ is related to the specific heat capacity cp, thermal expansion coefficient β, and speed of sound vs. Γ is a temperature-dependent, mainly due to the temperature dependence of β, which leads to the temperature mapping capability of PAT. Although vs in tissues is also temperature-dependent (0.1% to 0.3% increase per degree temperature rise), its effect on Γ is usually negligible compared with that of β (5% increase per degree temperature rise).5,6

Several PA methods have been developed to measure temperatures.1,3,711 In these methods, it is generally assumed that the optical fluence is independent of temperature and the temperature has a linear relationship with the Grüneisen parameter. Thus, by measuring the changes in PA signals at different temperatures over time, the ratiometric change in the PA signals as a function of temperatures can be obtained by fitting the data to a linear model.2,8,9,12 These ratiometric methods can obtain an accuracy of 0.1°C in relative temperature mapping and have been used in photothermal therapy,10,13 cryoablation of prostate tissue,14 and development of gold nanorods.7 However, these ratiometric methods are based on the photoacoustic microscopy systems and mostly work only in superficial depths. Moreover, they cannot provide absolute temperatures due to the unknown tissue optical properties and thus optical fluence. We previously proposed another temperature measurement method using the dual-temperature dependences of Grüneisen parameter and the speed of sound in tissue.11 The temperature irrelevant factors, which are difficult to correct for in deep tissue, such as the optical fluence and the optical absorption coefficient, are eliminated by taking ratiometric measurements at two different temperatures. Although this method realizes the absolute temperature measurement in deep tissue, it assumes a relatively homogeneous temperature distribution, which may not be valid for thermotherapy in which the temperature elevation is usually highly confined within a small tissue volume.

All of the above methods avoid directly quantifying the tissue’s optical parameters and thus optical fluence, which leads to the necessity of either calibration or ratiometric measurement under the strong assumption that the local optical fluence and optical parameters do not change with the temperature. Such an assumption may not be valid in clinical settings, largely due to the temperature-induced hemodynamic changes, including increased blood perfusion and oxygenation in response to the temperature elevation.1517 Since the initial acoustic pressure can be reconstructed by the measured PA signals, the absolute Grüneisen parameter distribution can be recovered once the optical fluence and optical parameters are reconstructed. Previous studies have introduced various optical models, such as the Monte Carlo simulation,18 radiative transfer equation (RTE),1921 and optical diffusion equation, to quantify the optical variables above.22,23 Generally, two assumptions are used in the previous studies to simplify the reconstruction problem: homogeneous Grüneisen parameter and homogeneous optical properties. Several studies focusing on reconstructing the optical parameters have been proposed,1820,2224 but temperature-dependent Grüneisen parameter is usually not considered. Other studies have focused on reconstructing optical parameters and Grüneisen parameter simultaneously, but prior information is required in these methods.2529 Due to the strong underdetermined nature of the inverse problem, all the unknown variables, including absorption and scattering coefficients, optical fluence, and Grüneisen parameter, cannot be uniquely recovered from a single illumination and a single-wavelength PA data without additional prior information.27 Two more studies have proposed to decouple the Grüneisen parameter from the reconstructed optical parameters using prior information such as known absorption or scattering distributions,27,28 which, however, are not always available in clinical settings.

Increasing the amount of measured data is a common strategy to eliminate the need for prior information and reduce the degree of underdetermination. Multi-illumination in sequence, which can introduce more independent measurements, is widely used in diffuse optical tomography,30 but few similar strategies have been used in PAT. Zemp et al.31 and Ranasinghesagara et al.32 proposed a method on reconstructing optical parameters using multiple surface illuminations but with the assumption of a known Grüneisen parameter. Bal and Ren presented a noniterative procedure to reconstruct the Grüneisen parameter but with the assumption of known absorption or scattering coefficients.27 Shao et al. realized the simultaneous reconstruction of the Grüneisen parameter and the optical parameters, but the Born approximation was assumed, which limited the applicability of this method to the weak optical heterogeneity.33

Here, we propose an absolute temperature mapping method enhanced by an optical model based on multi-illumination PAT. Instead of the traditional calibration-based temperature measurements over time, this method focuses on recovering the optically relevant variables to quantify the absolute temperature directly. An optical diffusion model is introduced to estimate the optical fluence. A Jacobian matrix-based method without the Born approximation is used to reconstruct other optical parameters. With the known initial pressure rise reconstructed from the PA signals, the absolute temperature distribution inside the tissue can be recovered. The assumption that the tissue’s optical properties are temperature-independent is not required in our method. Moreover, this new strategy can correlate the deep-tissue temperature with the surface temperature, which can be readily measured.




Photoacoustic Temperature Mapping

In PAT, for a nonfluorescent target, the local heating due to a laser pulse Si incident at surface location rsi leads to a thermoelastic expansion and initial acoustic pressure rise, which can be expressed as

Eq. (1)

where Hi(r) is the absorbed energy density at location r. Hi(r) can be expressed as

Eq. (2)

where μa(r) and ϕ(r,rsi) are the optical absorption coefficient and the optical fluence, respectively. The Grüneisen parameter Γ(r), which is also position relevant, is defined as2

Eq. (3)

where β is the thermal coefficient of volume expansion, κ is the isothermal compressibility, ρ is the mass density, CV is the specific heat capacity at constant volume, vs is the speed of sound, CP is the specific heat capacity at a constant pressure, and T(r) is the temperature distribution in the tissue. The relationship between Grüneisen parameter and temperature can be modeled by curve fitting the measured data. In this study, the linear relationship for soft tissue from a previous study by Alaeian is used,34 which is given as

Eq. (4)

where T is the temperature in Celsius. Therefore, the initial acoustic pressure can be expressed as

Eq. (5)

From Eq. (5), to recover the temperature distribution, the initial acoustic pressure rise distribution reconstructed by the measured acoustic pressure signals, the optical absorption coefficient μa(r), and the optical fluence ϕ(r,rsi) are all needed.


Forward Problem: Optical Model-Enhanced Photoacoustic Tomography

For an optically heterogeneous media, the task of temperature mapping is converted to the recovery of optical parameters including absorption coefficient μa(r), reduced scattering coefficient μs(r), and optical fluence ϕ(r,rsi) from the reconstructed initial acoustic pressure rise. In the forward problem, an optical model of the photon transportation in the tissue is introduced to quantify the distribution of optical fluence ϕ(r,rsi). For the highly scattering biological tissue, the diffusion approximation model can be used, which is expressed as4

Eq. (6)

where D(r)={3[μa(r)+μs(r)]}1 is the optical diffusion coefficient and qi(r) is the source term of Si. If the optical parameters μa(r) and D(r) are obtained, the optical fluence can be solved from Eq. (6) using finite element method (FEM).35,36

Therefore, the optical forward problem can be formed by a nonlinear mapping F, linking the optical parameters to the optical fluence distribution

Eq. (7)

Substituting Eq. (7) into Eq. (5), the forward problem of optical model enhanced PAT can be expressed as

Eq. (8)

where P denotes a nonlinear mapping that relates the temperature and the optical parameters to the initial acoustic pressure rise.


Inverse Problem: Reconstruction of Optical Parameter and Temperature Distribution by Using Multiple Illuminations

In this work, the initial acoustic pressure rise is reconstructed from the measured PA signals and is considered as accurate measured data.25,37,38 The inverse problem is to reconstruct the optical parameters [μa(r) and D(r)], Grüneisen parameter Γ(r), and temperature T(r), from the measured data p0(r,rSi) and optical fluence ϕ(r,rsi), which can be calculated by solving the optical forward problem. From Eq. (8), the generation of initial acoustic pressure rise can be divided into two independent components: the temperature-dependent Grüneisen parameter and the temperature-irrelevant absorbed energy density. It is an underdetermined problem to simultaneously recover the two independent components, which are both heterogeneous in the tissue, from a single measurement p0(r,rSi). It is worth noting that the tissue’s optical properties can change with temperatures as a result of physiological responses (e.g., elevated blood perfusion), but as physical parameters, they are not temperature-dependent quantities (e.g., the optical absorption cross section of hemoglobin is not temperature dependent). Again, unlike previous methods taking ratiometric measurement at different temperatures, our method does not assume temperature-independent optical properties of the tissues.


Multiple-illumination PAT

First, the optical fluence and tissue’s optical properties need to be reconstructed. The multiple illuminations in sequence are used to generate independent measurements and eliminate the unknown heterogeneous Grüneisen parameter, which will be described in the following sections. As shown in Fig. 1, multiple illuminations with point optical sources at different locations can be achieved by rotating the imaged object or moving the optical sources. The multiple illuminations are performed sequentially in the PA imaging process, meaning that one independent PA image is acquired at each illumination location. Under each illumination source (location), the initial acoustic pressure rise can be reconstructed as one measurement. Then, the measurements are equally grouped into two sets, k and l. Each set has the same number of optical illuminations N.

Fig. 1

Schematic of multiple illuminations in sequence. The red circles denote the different illumination locations. Sk and Sl denote the two equally separated sets of illuminations with the same number N.


After discretization, the measured data can be assembled into two matrixes

Eq. (9)


The ratio of two sets of measured data is formed by the ratio of each corresponding element in the matrix and can be derived into the ratio of optical fluence of two sets as follows:

Eq. (10)

where the optical fluence ϕk(r) and ϕl(r) have the same form of the measured data.

In Eq. (10), the heterogeneous temperature-relevant terms are canceled out and only a ratiometric measurement of optical fluence corresponding to different illumination locations is left. Therefore, we can divide the modeling into two steps: first, the optical parameters and optical fluence can be reconstructed by the ratio of multiple-illumination data, and second, the temperature distribution will be recovered.


Nonlinear iterative reconstruction of optical parameters

In general, the optical absorption coefficient and diffusion coefficient are both heterogeneous and unknown in this step. If the perturbation of the optical fluence caused by the heterogeneity of optical parameters is not negligible, as in the biological tissues, the reconstruction of optical parameters is a nonlinear process. Therefore, the optical parameters and optical fluence need to be reconstructed iteratively, with the guidance of the measured data. Briefly, an iteration reconstruction architecture can be designed as follows:

  • 1. Input the initial homogeneous optical parameters as published in the literatures.39

  • 2. Solve the optical forward problem to update the optical fluence at each illumination location.

  • 3. Calculate the error between the ratio of two sets of optical fluence and the ratio of the measured data, and estimate the convergence. If the error satisfies the termination condition, terminate the iteration loop; otherwise, continue to the next step.

  • 4. Construct and solve the inverse problem to update the optical parameters.

  • 5. Return to step 2.

In step (2), at each illumination location, the corresponding optical fluence can be calculated by the diffusion model described in Sec. 2.2, with the optical parameters obtained from step (4) in the previous iteration.

In step (3), the error can be calculated as

Eq. (11)

where fc=ϕk(r)ϕl(r) and fm=p0  k(r)p0l(r) denote the calculated and measured values, respectively. A predefined small quantity ϵ* can be used as the termination condition. If ϵϵ*, the iteration will be terminated.

In step (4), the inverse problem is set up to calculate the updated value of optical parameters. The updating direction is determined by approaching fm from fc. Therefore, fm is expressed by the first-order Taylor expansion with the matrix form

Eq. (12)


Eq. (13)

where the Jacobian matrix is

Eq. (14)

The partial derivations can be calculated by the diffusion approximation model. Then, we define b=[fmfc] and x=[δμaδD], Eq. (13) can be rewritten as the simple matrix form

Eq. (15)

The inverse problem is solving Eq. (15) for the updated optical parameters in the current iteration, which is a linear problem. Tikhonov regularization is used to reduce the ill-condition of the inverse problem, and the solution is given as40

Eq. (16)

where I is the identity matrix and λ is the regularization parameter.


Recovery of the temperature distribution

The temperature distribution can be recovered in the second step with the optical fluence and optical parameters reconstructed from the first step. According to Eq. (5), the heterogeneous Grüneisen parameter Γ(r) can be estimated. To reduce the effect of noise, the estimated Grüneisen parameter Γ^(r) can be derived using all the illuminations as follows:

Eq. (17)


Finally, the temperature distribution can be calculated by

Eq. (18)

Here, Eqs. (17) and (18) show that the Grüneisen parameter and the absolute temperature can be reconstructed by combining multi-illumination PAT with a nonlinear forward model and inverse model. The most practical aspect of the proposed method is the capability of quantifying the temperature inside the tissue without the need of modulating the tissue’s temperature. Moreover, since the temperature at the tissue surface can be accurately measured by a thermometer, this method provides a practical solution to calibrate the deep-tissue temperature using the surface temperature, which, to our best knowledge, is the first in PA temperature mapping.


Simulation Setups

Two-dimensional (2-D) simulations are performed to test the feasibility of the proposed method. The imaging field is discretized into 3575 nodes by FEM within a 4-cm-radius-circle.35,36 As shown in Figs. 2(a) and 2(b), true absorption and reduced scattering coefficients are mapped into the imaging field, respectively. The absorption coefficient of the background is 0.01  mm1 while the reduced scattering coefficient of the background is 1  mm1.30 Three absorption anomalies with an absorption coefficient of 0.03  mm1 and three scattering anomalies with a reduced scattering coefficient of 3  mm1 are added to create the optical heterogeneity. The isotropic point optical source is placed at 1mm inside the boundary of the imaging field, according to the approximation of the optical diffusion model.30 For the heterogeneity of temperature, three temperature anomalies at 40°C, 45°C, and 55°C are added to the homogeneous background temperature of 37°C. The initial acoustic pressure rise distributions are modeled by multiplying the distribution of simulated optical fluence,30 absorption coefficient, and Grüneisen parameter. As an example, Figs. 2(c) and 2(f) show the distribution of optical fluence and initial acoustic pressure rise for the optical source S1. For all the simulations, the iteration terminated condition is set to be ϵ*=0.001 in the reconstruction process. A control group is also set for comparison in each simulation. The control group uses the homogeneous empirical optical parameters in the diffusion model.

Fig. 2

Basic simulation setups. (a) 2-D absorption coefficient distribution. (b) 2-D reduced scattering coefficient distribution. (c) Optical fluence distribution corresponding to a representative illumination location S1. (d) 2-D temperature distribution. (e) 2-D Grüneisen parameter distribution. (f) Initial acoustic pressure rise distribution corresponding to S1.



Influence of number of illuminations

The number of illuminations can impact the degree of underdetermination of the inverse problem. To test the influence of the number of illuminations, four simulations are configured. As shown in Fig. 3, 2, 4, 8, and 16 illuminations are placed around the imaging field, respectively. In each simulation setup, the source locations are distributed with equiangular intervals. In this section, all the simulations are noise-free.

Fig. 3

Simulation setups to test the influence of the number of illuminations. (a)–(d) 2, 4, 8, and 16 illuminations around the imaging field. The red circles denote the locations of each illumination, distributed with equiangular intervals.



Influence of noise

Noise level is also studied to test the noise-resistance performance of the proposed method. In this section, three simulations with 1%, 5%, and 10% Gaussian noise are performed. In each simulation, the number of illumination is set to eight, an optimal value obtained from the previous simulation.



Figure 4 shows the reconstructed results of the first four simulations in Sec. 2.4.1, corresponding to different numbers of illumination. The reconstructed temperature is presented with truncated range from 37°C to 55°C [Fig. 4(d)] and the absolute values [Fig. 4(e)]. Due to the assumption of optical homogeneity, the results with a different number of illuminations are the same in the control group (Fig. 4, the second row). Clearly, the temperature distributions are recovered incorrectly in the control group due to a failure to account for local optical fluence. With the false assumption of optical properties, the large artifacts at the absorption anomalies may mislead the temperature monitoring. By contrast, using the proposed method, the reconstruction quality improves with the increase in the number of illuminations. The temperature anomalies can be distinguished from the background in each simulation, and no crosstalk among them has been observed in the reconstructed temperature maps. However, the reconstructed temperature distributions still suffer from artifacts induced by the anomalies of the optical properties, especially the absorption coefficients. With reduced number of illuminations, the artifacts at the absorption anomalies in the simulation become more significant. In the proposed method, the temperature is a dependent variable and the optical parameters are independent variables. Thus, the artifacts in the reconstructed temperature maps are largely due to the inaccurate optical parameters reconstructed from the forward and inverse models. The variations of the reconstructed optical parameters also provide similar results, showing that more illuminations lead to less crosstalk and artifacts in the reconstructed temperature images.

Fig. 4

Simulation results of the optical parameters and temperature with different number of illuminations. From left to right: (a) absorption coefficient, (b) reduced scattering coefficient, (c) optical fluence, and (d–e) temperature. From the top row to the bottom row: ground truth, control group, 2, 4, 8, and 16 illuminations.


The relative errors are calculated for quantitative assessments of the reconstructed results. Table 1 presents the relative errors of the reconstructed variables in each simulation, with different number of illuminations. The true distributions are used as the gold standards. The relative error of optical fluence is defined as the average relative error of optical fluence within the imaging field. The reconstructed temperature distribution has an accuracy of 90% in the case of 8 and 16 illuminations. The relative errors of all the reconstructed variables present similar trends with the number of illuminations. When the illumination number is larger than four, the accuracy is significantly improved. We conclude that the improvement in accuracy is attributed to the reduction of the degree of underdetermination with the increase in the number of illuminations.

Table 1

Relative errors of reconstructed results and time cost with different number of illuminations.

Number of illuminationsTime cost (s)Relative error (%)
TemperatureAbsorption coefficientReduced scattering coefficientAverage optical fluence

The time cost of computation in each simulation is also assessed. More illuminations lead to larger size of the Jacobian matrix, which results in longer computational time mainly in the inversion process. Of course, more illuminations also lead to a longer imaging time when either the sample or the light source is rotated. Therefore, there is a tradeoff between the number of illuminations and time cost when both accuracy and efficiency are taken into consideration (Fig. 5). The time cost shown in Table 1 mostly accounts for the data processing time, not the actual acquisition time of the data. Ten image frames can be acquired within 1 s for all illuminators, if a typical photoacoustic computed tomography system with a 10-Hz laser is used. The temperature change in typical thermotherapy within 1 s is thought to be minimal,41,42 so it is safe to assume that the temperature measurement by our method at a single temperature is much less affected by the hemodynamic changes than previous methods using ratiometric measurement at different temperatures.

Fig. 5

Relative errors of the reconstructed temperature and time cost with different number of illuminations.


Figure 6 shows the reconstructed results with different noise levels. For the control group, temperature reconstructions are more blurred with increased noise levels. Although the quality of reconstructed temperature maps is sensitive to the noise level, the temperature anomalies can be resolved even with 10% noise, reflecting the noise resistance of the proposed method. The reconstructed temperature anomalies can be distinguished from the background, and the temperature distribution inside the targets is relatively smooth. Similar to the simulations above, the reconstructed temperature distributions are affected by the accuracy of reconstructed optical parameters. With the increased noise level, more artifacts induced by the reconstructed optical parameters, especially from the reconstructed absorption coefficients, are shown in the background and blur the boundaries of the temperature anomalies.

Fig. 6

Simulation results of the optical parameters and temperature with different noise levels. From left to right: (a) absorption coefficient, (b) reduced scattering coefficient, (c) optical fluence, (d, e) temperature with eight illuminations, and (f, g) temperature of the control group. From the top row to the bottom row: Ground truth, 1% noise, 5% noise, and 10% noise.


Detailed quantitative results can be found in Table 2, which shows the relative errors of each reconstructed variable with different added noise levels. Optical fluence (shown by the average value) has higher reconstruction accuracy but larger variation range than optical parameters, whereas the reconstruction accuracy of temperature is always worse than that of optical parameters with different noise levels. These results again support the conclusion that in PA mapping of the absolute temperatures, the reconstruction accuracy of optical parameters mainly contributes to the accuracy of the temperature estimation.

Table 2

Relative reconstruction errors with different noise levels.

Noise levelTemperature (%)Absorption coefficient (%)Reduced scattering coefficient (%)Average optical fluence (%)



We have reported an absolute temperature mapping method using PAT, which can assess the temperature by accounting for the optical fluence and the optical parameters inside the tissue. The major improvement of the proposed method over previous PA temperature mapping is that we do not have to modulate the tissue’s temperature due to the assumption that the tissue’s optical properties or the optical fluence inside the tissue do not change with the temperature. Our method has two steps. In the first step, the ratiometric measurement of initial acoustic pressure rise at different illumination locations in sequence is calculated to eliminate the influence of heterogeneous Grüneisen parameters, which are temperature-dependent. An optical forward model is then introduced to recover the distribution of optical fluence with different illumination locations. In the second step, the optical parameters and the optical fluence are updated iteratively. After the iteration converges, the absolute temperature distribution can be recovered by dividing the initial acoustic pressure rise by the product of reconstructed optical fluence and the absorption coefficient. Results of 2-D simulations have shown the feasibility of the proposed method. Relative errors of 10.6% to 15.7% can be achieved at different noise levels. The absolute temperature mapping ability of this method may lead to potential applications of temperature monitoring during various cancer thermotherapies.

There are two major reasons to adopt the multi-illumination strategy. First, it is necessary to eliminate the unknown heterogeneous Grüneisen parameter in the forward optical model. Second, it is important to reduce the degree of underdetermination of the inverse problem. Different number of illuminations is used to test the performance of the proposed method. It is worth reiterating that the multiple illuminations are performed sequentially in the PA imaging process, meaning that one independent PA image is acquired at each illumination location. Since the measurements are based on the ratio of different illumination locations, at least two illuminations are required. The results shown in Fig. 4 and Table 1 demonstrate that a higher temperature mapping accuracy can be achieved with more illuminations. This accuracy improvement can be attributed to the reduction of underdetermination of the inverse problem. More independent information can enlarge the dimensions of the Jacobian matrix and reduce the degree of underdetermination.

The optical parameters need to be reconstructed iteratively. In each iteration, the values are updated by solving the inverse problem. Due to the strong irregularity of photon transport in the tissue, the inverse problem becomes highly ill-posed, which is the primary factor that impacts the accuracy of the temperature mapping in our method. Both Tables 1 and 2 show that the relative errors of temperature distributions are larger than those of optical parameters. Since the temperature is derived from initial acoustic pressure rise, optical fluence, and optical parameters, the accuracy of the temperature assessment will not be better than that of the optical variables in theory, especially with the noise. In Table 2, the relative errors are increased by 42.27% (for the absorption) and 61.18% (for the scattering) when the noise level changes from 1% to 10%. The regularization strategy used in this method relieves the ill-posed nature of the reconstruction problem and improves the accuracy. For future studies, other regularization methods such as L-p norm43 and total variation norm44 can be used to optimize the inverse problem.

In this study, we mainly focus on the optical perspectives of the physical model. Therefore, we assume an ideal ultrasound detection system, a simplification strategy used by other studies as well.1820,23,29 The reconstruction error of the initial acoustic pressure rise is not considered and is assumed to be accurate to monitor the direction of iteration in the temperature reconstruction. Nevertheless, the accuracy of the ultrasound reconstruction will eventually impact the accuracy of the temperature mapping in our method. To improve the ultrasound reconstruction, several important factors should be taken into consideration, including the speed of sound, the acoustic attenuation, and the frequency responses of the ultrasound transducer (array), as detailed by numerous previous studies.25,37,38,4554

Other factors of the PA imaging system are also not included in the mathematical model, such as the detection sensitivity of the ultrasound transducer (array). Based on the calculated temperature distribution given by this method, the deep-tissue temperature can be spatially calibrated by the absolute temperature directly measured at the surface to reduce the influence of the imaging system. So even though the calibration is still needed when the unknown factors of the imaging system are taken into consideration, this method provides a spatial calibration instead of changing the tissue’s temperature over time, which is more convenient for the clinical applications. Another factor contributing to the temperature reconstruction error is the forward model, which directly determines the accuracy of the optical fluence. The diffusion approximation model is used in this method in the case of optical heterogeneous media, which has modified the assumption of optical homogeneity in the traditional PAT studies. However, an important approximation of this forward model is that the optical fluence has to be dominated by diffuse photons, which holds only for the photons traveling deeper than one transport mean-free path (TMFP). This approximation leads to the inability of the diffusion approximation model to accurately describe the photon transportation in the region within one TMFP. From this perspective, both the optical fluence and the optical parameters are inaccurate even if the ill-posed nature of the inverse problem can be reduced. To achieve more accurate quantification of photon transportation, we will consider using the analytical models such as the RTE55 and delta-Eddington approximation,56 and a numerical method such as the Monte Carlo method.

Our method eventually utilizes a linear relationship between the Grüneisen parameter and the temperature, which, however, is not a necessary condition. Based on the literature, the linear temperature-dependence range is approximately from 4°C to 55°C,1 which is actually achievable in most thermotherapy. Several studies have reported on the nonlinear temperature dependence of the Grüneisen parameter when the temperature exceeds the linear range.5759 In our method, the temperature is recovered from the estimated Grüneisen parameter only at the last step of the process, so a nonlinear temperature dependence may also work in our model.

Computational time cost is also an important performance index of the method, especially for applications that need fast temperature mapping. The results in Sec. 2.4.1 show the tradeoff between the reconstruction accuracy and the time cost. Moreover, a longer computational time will be taken if the reconstruction is extended to three-dimensional case. The number of illuminators and the number of the discretized elements of the solution domain, which contribute to the scale of the Jacobian matrix, primarily determine the time cost of solving the forward and inverse problems. The iteration algorithm for solving the inverse problem also impacts the time cost of the reconstruction process. Therefore, for time-sensitive applications, the time cost of our method needs to be reduced by optimizing the reconstruction algorithm in the inverse problem60,61 or using parallel computing techniques.62


All authors have no financial interest.


This work was supported by Duke MEDx fund (to J. Yao), National Natural Science Foundation of China (81227901, 81271617, 61322101, and 61361160418), and National Major Scientific Instrument and Equipment Development Project (2011YQ030114). The authors thank Chuangjian Cai, Mucong Li, and Yuqi Tang for their helpful discussion.



I. V. Larina, K. V. Larin and R. O. Esenaliev, “Real-time optoacoustic monitoring of temperature in tissues,” J. Phys. D: Appl. Phys., 38 (15), 2633 –2639 (2005). JPAPBE 0022-3727 Google Scholar


M. Pramanik and L. V. Wang, “Thermoacoustic and photoacoustic sensing of temperature,” J. Biomed. Opt., 14 (5), 054024 (2009). JBOPFO 1083-3668 Google Scholar


S.-H. Wang et al., “Photoacoustic temperature measurements for monitoring of thermal therapy,” Proc. SPIE, 7177 71771S (2009). Google Scholar


B. Cox et al., “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt., 17 (6), 061202 (2012). JBOPFO 1083-3668 Google Scholar


N. Bilaniuk and G. S. Wong, “Speed of sound in pure water as a function of temperature,” J. Acoust. Soc. Am., 93 (3), 1609 –1612 (1993). Google Scholar


H. Watanabe, N. Yamada and M. Okaji, “Linear thermal expansion coefficient of silicon from 293 to 1000 K,” Int. J. Thermophys., 25 (1), 221 –236 (2004). IJTHDY 0195-928X Google Scholar


Y. S. Chen et al., “Sensitivity enhanced nanothermal sensors for photoacoustic temperature mapping,” J. Biophotonics, 6 (6–7), 534 –542 (2013). Google Scholar


H. Ke et al., “Temperature mapping using photoacoustic and thermoacoustic tomography,” Proc. SPIE, 8223 82230T (2012). Google Scholar


M. Pramanik et al., “Tissue temperature monitoring using thermoacoustic and photoacoustic techniques,” Proc. SPIE, 7564 75641Y (2010). Google Scholar


X. Wu et al., “Photoacoustic-imaging-based temperature monitoring for high-intensity focused ultrasound therapy,” in IEEE 38th Annual Int. Conf. of the Engineering in Medicine and Biology Society (EMBC ’16), 3235 –3238 (2016). Google Scholar


J. Yao et al., “Absolute photoacoustic thermometry in deep tissue,” Opt. Lett., 38 (24), 5228 –5231 (2013). OPLEDP 0146-9592 Google Scholar


H. Ke, S. Tai and L. V. Wang, “Photoacoustic thermography of tissue,” J. Biomed. Opt., 19 (2), 026003 (2014). JBOPFO 1083-3668 Google Scholar


J. Shah et al., “Photoacoustic imaging and temperature measurement for photothermal cancer therapy,” J. Biomed. Opt., 13 (3), 034024 (2008). JBOPFO 1083-3668 Google Scholar


E. V. Petrova et al., “In vivo cryoablation of prostate tissue with temperature monitoring by optoacoustic imaging,” Proc. SPIE, 9708 97080G (2016). Google Scholar


D. K. Kelleher and P. Vaupel, “No sustained improvement in tumor oxygenation after localized mild hyperthermia,” Adv. Exp. Med. Biol., 662 393 –398 (2010). Google Scholar


C. W. Song, H. Park and R. J. Griffin, “Improvement of tumor oxygenation by mild hyperthermia,” Radiat. Res., 155 (4), 515 –528 (2001).[0515:IOTOBM]2.0.CO;2 Google Scholar


C. W. Song et al., “Improvement of tumor oxygenation status by mild temperature hyperthermia alone or in combination with carbogen,” Semin. Oncol., 24 (6), 626 –632 (1997). Google Scholar


B. T. Cox et al., “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt., 45 (8), 1866 –1875 (2006). APOPAI 0003-6935 Google Scholar


L. Yao, Y. Sun and H. Jiang, “Quantitative photoacoustic tomography based on the radiative transfer equation,” Opt. Lett., 34 (12), 1765 –1767 (2009). OPLEDP 0146-9592 Google Scholar


L. Yao, Y. Sun and H. Jiang, “Transport-based quantitative photoacoustic tomography: simulations and experiments,” Phys. Med. Biol., 55 (7), 1917 –1934 (2010). PHMBA7 0031-9155 Google Scholar


M. Haltmeier, L. Neumann and S. Rabanser, “Single-stage reconstruction algorithm for quantitative photoacoustic tomography,” Inverse Prob., 31 (6), 065005 (2015). INPEEY 0266-5611 Google Scholar


Z. Yuan and H. Jiang, “Quantitative photoacoustic tomography: recovery of optical absorption coefficient maps of heterogeneous media,” Appl. Phys. Lett., 88 (23), 231101 (2006). APPLAB 0003-6951 Google Scholar


A. Pulkkinen et al., “Direct estimation of optical parameters from photoacoustic time series in quantitative photoacoustic tomography,” IEEE Trans. Med. Imaging, 35 (11), 2497 –2508 (2016). ITMID4 0278-0062 Google Scholar


B. Cox, S. Arridge and P. Beard, “Gradient-based quantitative photoacoustic image reconstruction for molecular imaging,” Proc. SPIE, 6437 64371T (2007). PSISDG 0277-786X Google Scholar


M. Fink et al., “Time-reversed acoustics,” Rep. Prog. Phys., 63 (12), 1933 –1995 (2000). RPPHAG 0034-4885 Google Scholar


Y. Hristova, P. Kuchment and L. Nguyen, “Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,” Inverse Prob., 24 (5), 055006 (2008). INPEEY 0266-5611 Google Scholar


G. Bal and K. Ren, “Multi-source quantitative photoacoustic tomography in a diffusive regime,” Inverse Prob., 27 (7), 075003 (2011). INPEEY 0266-5611 Google Scholar


G. Bal and K. Ren, “On multi-spectral quantitative photoacoustic tomography in diffusive regime,” Inverse Prob., 28 (2), 025010 (2012). INPEEY 0266-5611 Google Scholar


B. T. Cox, S. R. Arridge and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A, 26 (2), 443 –455 (2009). Google Scholar


L. V. Wang and H.-I. Wu, Biomedical Optics: Principles and Imaging, John Wiley & Sons, Hoboken, New Jersey (2012). Google Scholar


R. J. Zemp et al., “A photoacoustic method for optical scattering measurements in turbid media,” Proc. SPIE, 7177 71770Q (2009). PSISDG 0277-786X Google Scholar


J. C. Ranasinghesagara et al., “Photoacoustic technique for assessing optical scattering properties of turbid media,” J. Biomed. Opt., 14 (4), 040504 (2009). JBOPFO 1083-3668 Google Scholar


P. Shao, B. Cox and R. J. Zemp, “Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography,” Appl. Opt., 50 (19), 3145 –3154 (2011). APOPAI 0003-6935 Google Scholar


M. Alaeian, “Parameter and temperature estimation using photoacoustic technique based on computational simulation,” (2017). Google Scholar


S. R. Arridge, “Optical tomography in medical imaging,” Inverse Prob., 15 (2), R41 (1999). INPEEY 0266-5611 Google Scholar


S. Arridge et al., “A finite element approach for modeling photon transport in tissue,” Med. Phys., 20 (2), 299 –309 (1993). MPHYA6 0094-2405 Google Scholar


M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, 71 (1), 016706 (2005). PLEEE8 1539-3755 Google Scholar


M. Fink, “Time reversal of ultrasonic fields. I. Basic principles,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 39 (5), 555 –566 (1992). Google Scholar


D. Hyde et al., “Performance dependence of hybrid x-ray computed tomography/fluorescence molecular tomography on the optical forward problem,” J. Opt. Soc. Am. A, 26 (4), 919 –923 (2009). Google Scholar


A. Neumaier, “Solving ill-conditioned and singular linear systems: a tutorial on regularization,” SIAM Rev., 40 (3), 636 –666 (1998). SIREAD 0036-1445 Google Scholar


A. Jordan et al., “The effect of thermotherapy using magnetic nanoparticles on rat malignant glioma,” J. Neuro-Oncol., 78 (1), 7 –14 (2006). JNODD2 0167-594X Google Scholar


S. Maenosono and S. Saita, “Theoretical assessment of FePt nanoparticles as heating elements for magnetic hyperthermia,” IEEE Trans. Magn., 42 (6), 1638 –1642 (2006). IEMGAQ 0018-9464 Google Scholar


S. Okawa, Y. Hoshi and Y. Yamada, “Improvement of image quality of time-domain diffuse optical tomography with lp sparsity regularization,” Biomed. Opt. Express, 2 (12), 3334 –3348 (2011). BOEICL 2156-7085 Google Scholar


J. Dutta et al., “Joint L1 and total variation regularization for fluorescence molecular tomography,” Phys. Med. Biol., 57 (6), 1459 –1476 (2012). PHMBA7 0031-9155 Google Scholar


C. Hoelen et al., “Three-dimensional photoacoustic imaging of blood vessels in tissue,” Opt. Lett., 23 (8), 648 –650 (1998). OPLEDP 0146-9592 Google Scholar


M. Pramanik, “Improving tangential resolution with a modified delay-and-sum reconstruction algorithm in photoacoustic and thermoacoustic tomography,” J. Opt. Soc. Am. A, 31 (3), 621 –627 (2014). Google Scholar


S. Park et al., “Adaptive beamforming for photoacoustic imaging,” Opt. Lett., 33 (12), 1291 –1293 (2008). OPLEDP 0146-9592 Google Scholar


D. Yang et al., “Fast multielement phase-controlled photoacoustic imaging based on limited-field-filtered back-projection algorithm,” Appl. Phys. Lett., 87 (19), 194101 (2005). APPLAB 0003-6951 Google Scholar


P. Burgholzer et al., “Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors,” Inverse Prob., 23 (6), S65 (2007). INPEEY 0266-5611 Google Scholar


Y. Wang et al., “Photoacoustic imaging with deconvolution algorithm,” Phys. Med. Biol., 49 (14), 3117 –3124 (2004). PHMBA7 0031-9155 Google Scholar


L. Zeng et al., “High antinoise photoacoustic tomography based on a modified filtered backprojection algorithm with combination wavelet,” Med. Phys., 34 (2), 556 –563 (2007). MPHYA6 0094-2405 Google Scholar


J. Laufer et al., “Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt., 49 (8), 1219 –1233 (2010). APOPAI 0003-6935 Google Scholar


H. Grün et al., “Photoacoustic tomography using a fiber based Fabry-Perot interferometer as an integrating line detector and image reconstruction by model-based time reversal method,” in European Conf. on Biomedical Optics, (2007). Google Scholar


X. L. Dean-Ben et al., “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging, 31 (10), 1922 –1928 (2012). ITMID4 0278-0062 Google Scholar


A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, New York (1978). Google Scholar


W. Cong et al., “Modeling photon propagation in biological tissues using a generalized Delta-Eddington phase function,” Phys. Rev. E, 76 (5), 051913 (2007). PLEEE8 1539-3755 Google Scholar


S. M. Nikitin, T. D. Khokhlova and I. M. Pelivanov, “Temperature dependence of the optoacoustic transformation efficiency in ex vivo tissues for application in monitoring thermal therapies,” J. Biomed. Opt., 17 (6), 061214 (2012). JBOPFO 1083-3668 Google Scholar


R. Brinkmann et al., “Real-time temperature determination during retinal photocoagulation on patients,” J. Biomed. Opt., 17 (6), 061219 (2012). JBOPFO 1083-3668 Google Scholar


E. Petrova et al., “Imaging technique for real-time temperature monitoring during cryotherapy of lesions,” J. Biomed. Opt., 21 (11), 116007 (2016). JBOPFO 1083-3668 Google Scholar


D. Kingma and J. Ba, “Adam: a method for stochastic optimization,” (2014). Google Scholar


X. Cao et al., “Accelerated image reconstruction in fluorescence molecular tomography using dimension reduction,” Biomed. Opt. Express, 4 (1), 1 –14 (2013). BOEICL 2156-7085 Google Scholar


M. Chen et al., “Fast direct reconstruction strategy of dynamic fluorescence molecular tomography using graphics processing units,” J. Biomed. Opt., 21 (6), 066010 (2016). JBOPFO 1083-3668 Google Scholar


Yuan Zhou is a PhD student at the Department of Biomedical Engineering, Tsinghua University, Beijing, China. He is currently a visiting scholar at the Department of Biomedical Engineering, Duke University, North Carolina, USA. His research interests are modeling and algorithms in optical molecular imaging and photoacoustic imaging.

Eric Tang is an undergraduate student at the Department of Biomedical Engineering and Computer Science, Duke University, North Carolina, USA. He currently works at Duke Photoacoustic Imaging Lab.

Jianwen Luo received his BS and PhD degrees in biomedical engineering from Tsinghua University, Beijing, China, in 2000 and 2005, respectively. He became a Principal Investigator (Tenure Track Associate Professor) in the Department of Biomedical Engineering and the Center for Biomedical Imaging Research at Tsinghua University in 2011, and obtained his tenure position (associate professor) in 2017. His research interest includes ultrasound imaging, fluorescence imaging, and signal processing. He is a senior member of IEEE and serves as an advisory editorial board member of Journal of Ultrasound in Medicine and associate editor of Medical Physics.

Junjie Yao is an assistant professor of biomedical engineering at Duke University, North Carolina, USA, and a faculty member of Duke Center for In-Vivo Microscopy, Fitzpatrick Institute for Photonics, and Duke Cancer Institute. He received his BE and ME degrees at Tsinghua University, and his PhD in biomedical engineering at Washington University, St. Louis, Missouri, USA. His research interest is in photoacoustic tomography technologies in life sciences, especially in functional brain imaging and early cancer detection.

© 2018 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2018/$25.00 © 2018 SPIE
Yuan Zhou, Eric Tang, Jianwen Luo, and Junjie Yao "Deep-tissue temperature mapping by multi-illumination photoacoustic tomography aided by a diffusion optical model: a numerical study," Journal of Biomedical Optics 23(1), 016014 (29 January 2018).
Received: 26 August 2017; Accepted: 11 January 2018; Published: 29 January 2018

Back to Top