1.

Introduction

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.^1–3 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.⁴ The Grüneisen parameter $\mathrm{\Gamma }$ is related to the specific heat capacity $c_p$, thermal expansion coefficient $\beta$, and speed of sound $v_s$. $\mathrm{\Gamma }$ is a temperature-dependent, mainly due to the temperature dependence of $\beta$, which leads to the temperature mapping capability of PAT. Although $v_s$ in tissues is also temperature-dependent ($\sim 0.1\%$ to 0.3% increase per degree temperature rise), its effect on Γ is usually negligible compared with that of $\beta$ ($\sim 5\%$ increase per degree temperature rise).^5,6

Several PA methods have been developed to measure temperatures.^1,3,7–11 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 $\sim 0.1°\mathrm{C}$ in relative temperature mapping and have been used in photothermal therapy,^10,13 cryoablation of prostate tissue,¹⁴ and development of gold nanorods.⁷ 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.¹¹ 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.^15–17 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,¹⁸ radiative transfer equation (RTE),^19–21 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,^{18–20,22–24} 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.^25–29 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.²⁷ 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,³⁰ but few similar strategies have been used in PAT. Zemp et al.³¹ and Ranasinghesagara et al.³² 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.²⁷ 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.³³

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.

2.

Methods

2.3.

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 [${\mu }_a(r)$ and $D(r)$], Grüneisen parameter $\mathrm{\Gamma }(r)$, and temperature $T(r)$, from the measured data $p_0(r,r_{S_i})$ and optical fluence $\mathrm{\varphi }(r,r_{s_i})$, 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 $p_0(r,r_{S_i})$. 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.

2.3.1.

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. $S_k$ and $S_l$ 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)

$$\begin{gathered} p_{0k}(r)={[p_0(r,r_{S_1}),p_0(r,r_{S_2}),\dots ,p_0(r,r_{S_I})]}^T,\{S_1,S_2,\dots ,S_N\}\in k, \\ {p}_{0l}(r)={[p_0(r,r_{S_{I+1}}),p_0(r,r_{S_{I+2}}),\dots ,p_0(r,r_{S_{2I}})]}^T,\{S_{N+1},S_{N+2},\dots ,S_{2N}\}\in l. \end{gathered}$$

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)

$$\begin{gathered} \frac{p_{0k}(r)}{p_{0l}(r)}=\frac{P_k[T(r),{\mu }_a(r),D(r)]}{P_l[T(r),{\mu }_a(r),D(r)]} \\ =\frac{{\mathrm{\varphi }}_k(r)}{{\mathrm{\varphi }}_l(r)}, \end{gathered}$$

where the optical fluence ${\mathrm{\varphi }}_k(r)$ and ${\mathrm{\varphi }}_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.

2.3.2.

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.³⁹

• 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)

$$ϵ=\frac{{\Vert f_c-f_m\Vert }_2^2}{{\Vert f_m\Vert }_2^2},$$

where $f_c=\frac{{\mathrm{\varphi }}_k(r)}{{\mathrm{\varphi }}_l(r)}$ and $f_m=\frac{p_{0k}(r)}{p_{0l}(r)}$ denote the calculated and measured values, respectively. A predefined small quantity ${ϵ}^{*}$ can be used as the termination condition. If $ϵ\le {ϵ}^{*}$, 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 $f_m$ from $f_c$. Therefore, $f_m$ is expressed by the first-order Taylor expansion with the matrix form

Eq. (12)

$$[f_m]=[f_c]+\left[\frac{\partial f_c}{\partial {\mu }_a}\right][\delta {\mu }_a]+\left[\frac{\partial f_c}{\partial D}\right][\delta D],$$

Eq. (13)

$$J·\left[\begin{array}{l}\delta {\mu }_a\\ \delta D\end{array}\right]=[f_m-f_c],$$

where the Jacobian matrix is

Eq. (14)

$$[J]=\left[\begin{array}{l}\frac{\partial f_c}{\partial {\mu }_a}\\ \frac{\partial f_c}{\partial D}\end{array}\right]=\left[\begin{array}{l}\frac{\partial {\mathrm{\varphi }}_k}{\partial {\mu }_a}·\frac{1}{{\mathrm{\varphi }}_l}-\frac{{\mathrm{\varphi }}_k}{{\mathrm{\varphi }}_l^2}·\frac{\partial {\mathrm{\varphi }}_l}{\partial {\mu }_a}\\ \frac{\partial {\mathrm{\varphi }}_k}{\partial D}·\frac{1}{{\mathrm{\varphi }}_l}-\frac{{\mathrm{\varphi }}_k}{{\mathrm{\varphi }}_l^2}·\frac{\partial {\mathrm{\varphi }}_l}{\partial D}\end{array}\right].$$

The partial derivations can be calculated by the diffusion approximation model. Then, we define $b=[f_m-f_c]$ and $x=\left[\begin{array}{l}\delta {\mu }_a\\ \delta D\end{array}\right]$, Eq. (13) can be rewritten as the simple matrix form

Eq. (15)

$$Jx=b.$$

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 as⁴⁰

Eq. (16)

$$x={(J^{T}J+\lambda I)}^{-1}{J}^{T}b,$$

where $I$ is the identity matrix and $\lambda$ is the regularization parameter.

2.3.3.

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 $\mathrm{\Gamma }(r)$ can be estimated. To reduce the effect of noise, the estimated Grüneisen parameter $\widehat{\mathrm{\Gamma }}(r)$ can be derived using all the illuminations as follows:

Eq. (17)

$$\widehat{\mathrm{\Gamma }}(r)=\frac{\sum _i{p}_{0i}(r)+\sum _j{p}_{0j}(r)}{{\mu }_a(r)[\sum _i\mathrm{\varphi }(r,r_{s_i})+\sum _j\mathrm{\varphi }(r,r_{s_j})]}.$$

Finally, the temperature distribution can be calculated by

Eq. (18)

$$T(r)=f^{-1}[\mathrm{\Gamma }(r)].$$

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.

2.4.

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{\mathrm{mm}}^{-1}$ while the reduced scattering coefficient of the background is $1{\mathrm{mm}}^{-1}$.³⁰ Three absorption anomalies with an absorption coefficient of $0.03{\mathrm{mm}}^{-1}$ and three scattering anomalies with a reduced scattering coefficient of $3{\mathrm{mm}}^{-1}$ 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.³⁰ 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,³⁰ 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 $S_1$. 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$.

2.4.1.

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.

2.4.2.

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.

3.

Results

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 $\sim 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.

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.

4.

Discussion

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 norm⁴³ and total variation norm⁴⁴ 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.^{18–20,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,45–54}

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 RTE⁵⁵ and delta-Eddington approximation,⁵⁶ 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,¹ 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.^57–59 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 problem^60,61 or using parallel computing techniques.⁶²