Pixel-wise reconstruction of tissue absorption coefficients in photoacoustic tomography using a non-segmentation iterative method.

In Photoacoustic Tomography (PAT), the acquired image represents a light energy deposition map of the imaging object. For quantitative imaging, the PAT image is converted into an absorption coefficient ( μ a ) map by dividing the light fluence (LF). Previous methods usually assume a uniform tissue μ a distribution, and consequently degrade the LF correction results. Here, we propose a simple method to reconstruct the pixel-wise μ a map. Our method is based on a non-segmentation-based iterative algorithm, which alternately optimizes the LF distribution and the μ a map. Using simulation data, as well as phantom and animal data, we implemented our algorithm and compared it to segmentation-based correction methods. The results show that our method can obtain accurate estimation of the LF distribution and therefore improve the image quality and feature visibility of the μ a map. Our method may facilitate efficient calculation of the concentration distributions of endogenous and exogenous agents in vivo.

Introduction

Photoacoustic tomography (PAT) or optoacoustic tomography (OAT), as a hybrid biomedical imaging modality, is able to acquire the instant light energy deposition inside an imaged object by detecting and processing the ultrasonic signal generated by laser illumination [1], [2], [3], [4]. Utilizing multiple wavelength excitation, multispectral PAT can distinguish the distribution of endogenous tissue absorbers, such as oxyhemoglobin (HbO2) and de-oxyhemoglobin (Hb), and exogenous optical probes such as the FDA-approved Indocyanine Green (ICG) based on their specific absorption spectrum [5]. Benefiting from high imaging speed, good detection sensitivity and centimeter-scale imaging depth, PAT has been used in various preclinical researches and clinical trials [6], [7], [8], [9], [10].

Within a given voxel, the PAT image intensity is proportional to the absorbed light energy density, which is the product of the light fluence (LF) arriving at that voxel and the absorption coefficient () of the enclosed absorbers [11]. Deriving the distribution of the estimated by removing the LF distribution can help obtain a quantitative concentration map of the absorbers [12]. Although there are currently some methods that can achieve quantitative blood oxygen saturation (sO2) measurement without LF correction (e.g., acoustic-spectrum-based methods [13], [14] and the absorption-saturation-based method [15]), estimating and correcting for the LF can further improve the accuracy of the measurement. Indeed, the true LF distribution of the imaged object is very difficult to measure, and thus is usually approximated by using light transport models. These models include the radiative transfer equation (RTE) [16], [17], diffusion equation (DE, a simplified form of the RTE) [18], [19], and numerical simulations based on the Monte Carlo method [20], [21]. Despite this, the estimation of the LF distribution relies on the prior information of the absorption distribution, which poses challenges for quantitative PAT imaging.

Previous studies for quantitative PAT imaging can be divided into experimental methods and simulation-based methods [11]. Experimental methods aim to measure the absorption and scattering of light within samples using external gauging devices. For example, Bauer et al. proposed the use of diffuse optical tomography to measure the surface scattered light and then used it to infer the LF distribution within the sample [22]. Steenbergen et al. used acousto-optic tomography to measure the LF distribution [23]. These methods required additional instruments and thus their application in PAT was limited. Simulation-based methods estimate the LF distribution based on the acquired PAT image or combined with co-registered images of other modalities, such as ultrasound (US) [24] and magnetic resonance imaging (MRI) [18]. Representative simulation-based techniques include the direct correction method [25], fixed-point iteration method [12], and model-based methods [26], [27], [28]. Banerjee et al. proposed the direct correction method to obtain the LF distribution by assuming that the reduced scattering coefficient () is known a priori and thus the diffusion coefficient [25]. Cox et al. proposed the fixed-point iteration method to obtain map without making linear assumptions [12], but the method diverges quickly if the prior knowledge of the scattering distribution is incorrect [11], [29]. The model-based minimization methods aim to minimize a functional quantifying the difference between the model output and the measured data by adjusting the distribution. For example, Jiang et al. applied the Gauss-Newton method to minimize the objective function and tackle the inversion of [26]. Cox et al. achieved greater computational efficiency by accelerating the calculation of the functional gradient vector with an adjoint model, but at the cost of requiring more iterations to converge [27], [28]. Compared to other techniques, the model-based methods usually yield better results.

To alleviate the computational complexity, segmentation-based methods have been proposed. Such methods generally assume that the imaging object has the same optical parameters, and require segmentation of the object contours. After assigning initial optical parameters to the segmented image, LF estimation can be performed using light transport models. Based on their solution algorithms, the simulation-based methods can be divided into two categories: segmentation-based direct correction (SBDC) methods, and segmentation-based iterative correction (SBIC) methods. The SBDC methods, such as Ref [30], directly obtained the map by dividing the PAT image by the obtained LF map [Fig. 1(c)]. However, these methods require manual selection of ideal initial optical parameters for each segmented region. The SBIC methods first segment the object contours in the PAT image and assume a uniform value [Fig. 1(d)], and then iteratively estimates the LF distribution by minimizing the error between the acquired PAT image and an estimated PAT image. To ease the manual intervention requirement of the simulation-based methods, Liang et al. proposed an automatic segmentation algorithm to obtain the animal contours for LF estimation [31]. Brochu et al. proposed organ-level segmentation to improve the accuracy of LF estimation [32]. However, the poor soft-tissue contrast of PAT images makes precise organ segmentation challenging. Recently, Pattyn et al. proposed to use co-registered US images to segment phantom images [24]. Zhang et al. proposed to use co-registered MR images to obtain a more refined segmentation of animal organs to improve LF estimation [33]. These attempts focused on improving segmentation efficiency and accuracy. However, they all started from the assumption of a uniform distribution within a region (or organ), whereas tissue values are actually spatially inhomogeneous. This inaccurate assumption leads to errors in the estimation of distribution. At present, there is still a lack of simple and convenient method that can accurately reconstruct a pixel-wise map for PAT.

Fig. 1

Schematic diagram of quantitative reconstruction in PAT. (a) Data acquisition: multispectral PAT images are acquired. (b) Un-corrected: un-corrected PAT images. (c) SBDC: segmentation-based direct correction. (d) SBIC: segmentation-based iterative correction. (e) Proposed: our proposed non-segmentation iterative algorithm. Segmentation: prior images for SBDC and SBIC methods. Initialization: initialization of values. reconstruction: Schematic diagram of different methods to obtain images. The fluence forward model is employed to estimate LF distribution. Spectral un-mixing: identify the distribution of endogenous absorbers (HbO2, Hb) from the background. Scale bar, 3 mm.

To meet the above challenges, we propose a simple, non-segmentation method to directly reconstruct the spatially varying distribution. We treat the tissue distribution as a pixel-wise spatially varying map and thus avoid the need for segmentation. Unlike previous optimization schemes for numerical models [16], [34], our method does not need to calculate Hessian matrices or gradient vectors of the objective function, but instead employs a two-step iterative algorithm to alternately optimize the LF distribution and the map [Fig. 1(e)]. We tested our method in numerical simulation experiments, tissue-mimicking phantom experiments, and live animal experiments. The results show that, compared with the SBDC and SBIC methods, our method reconstructs an accurate map that removes the signal inhomogeneity induced by light fluence attenuation. Based on the obtained map, spectral unmixing experiments on in vivo animals further show enhanced visualization of endogenous and exogenous contrast agents in deep tissues. This simple and convenient reconstruction method has the potential for the efficient quantification of absorber concentration in future PAT applications.

Method

The imaging model of PAT

The PAT images are formed by reconstructing the original point source of the ultrasonic waves generated by absorbing the laser pulse, and the pixel value in the image can be expressed as:where Γ denotes the thermo-elastic Grüneisen parameter, and denote the absorption and reduced scattering coefficients, Φ denotes the LF within a voxel at the position . In the beginning, we make an assumption that has been reconstructed from the acoustic measurements accurately and with negligible structural distortion [35], [36], [37], [38], [39]. In biological soft tissue, the thermo-elastic Grüneisen parameter has a small variation so that it is assumed to be constant [32]. Then, the distribution of absorbed energy is obtained by eliminating the constant terms from . Finally, the PAT image can be expressed as the product of the LF distribution Φ and the map:

However, solution of from Eq. (2) is considered difficult because the LF matrix Φ depends on both and . is the reduced scattering coefficient, which can be calculated by , where represents the scattering coefficient, represents the scattering anisotropy. Thus, it is often assumed that is known [12], [32], [40], [41]. Accordingly, the imaging model can be expressed as:where represents the mapping function that relates the map to the PAT image . The objective of quantitative PAT imaging is then to formulate an inverse mapping to obtain , which can be expressed as:

The segmentation-based LF correction methods

Since organ-level segmentation is difficult to achieve in animal experiments, all the segmentation-based methods involved in this paper segment the object contours and assume a uniform optical parameter within the whole imaging object. The SBDC method assumes a uniform distribution of within the imaged object [Fig. 1(c)], and the LF map can be estimated by manually selecting the initial value and then obtain the estimated () image [30], [31]. It can be expressed as [Fig. 1(c)]:

The SBIC method first segments the object contours in the PAT image and assumes a uniform value [Fig. 1(d)], and then iteratively calculate the value to minimize the error between the product and the un-corrected PAT image [32], [33]. It can be expressed as [Fig. 1(d)]:

The LF map can be estimated using the final value . Then, by dividing the un-corrected PAT image by the LF map, we can obtain the map:

The proposed iterative algorithm

With Eq. (3), the imaging model can be represented as a nonlinear model:where and are in vector representation, Φ is in sparse matrix representation. We can then formulate the solution of as a least square optimization problem:

The direct solution of the above objective function is difficult because the LF matrix Φ depends on μa. However, if we can solve Φ and μa one by one and then alternately update each of them, we might be able to reach their optimized solutions. Based on this idea, a two-step iterative optimization method based on gradient descent was developed [Fig. 1(e)] to solve the nonlinear model in Eq. (9). To start with, the μa map is initialized as 0. During each iteration, we first use the updated μa map to obtain the LF matrix Φ. Then, the objective function is converted into a linear model, from which the μa map can be updated by gradient descent algorithm. These two steps are repeated until convergence. The detailed steps of our method are summarized in Algorithm 1.

In the above two-step iterative algorithm, the LF distribution was modeled in the 2D plane by the DE, which was performed in MATLAB (The MathWorks, Inc., USA) using the open-source finite-element-based LF simulation software Toast+ + toolbox [42]. During the optical transport modeling process, the background medium is set as water and the map is initialized as 0 [Fig. 1(e)]. For all experiments, the Iter max is set to 30 and ε is set to 10−12. The average computation time for a single-slice image is 92.49 s

Compared to previous SBDC and SBIC methods, our two-step iterative algorithm solves for a pixel-wise map rather than a region-wise map. It does not require manual segmentation of the sample as well as initialization with a coarse estimation, and therefore avoids subjective variation of the observers. Also, allowing the map to vary for each image pixel, the objective function of Eq. (9) should converge to a more ideal solution.

Spectral un-mixing

The reconstructed image is a combination of various endogenous and exogenous absorbers. There is a linear relationship between the image and the value of the identified chromophore i:where, is the concentration of chromophore i and is the molar extinction coefficient at λ wavelength. In multispectral PAT, the core idea of spectral separation is to decompose the distribution of each absorber by identifying the specific absorption curve of each absorber. In this paper, we choose the most commonly used linear un-mixing method to calculate the distribution of HbO2, Hb and ICG [43]. After obtaining the concentrations of HbO2 () and Hb (), sO2 can be computed as:

Experimental setups

PAT imaging

A commercial small animal multispectral optoacoustic tomography system (MSOT inVision128, iThera Medical, Germany) was employed for imaging. Fig. 2(a) shows a diagram of the imaging chamber of the MSOT system. Pulsed laser (670–960 nm tunable) with pulse width < 10 ns, repetition rate of 10 Hz, and a peak pulse energy of 60 mJ at 760 nm excite the sample through a ten-arm fiber bundle, which provides homogeneous, 360-degree illumination of approximately 8 mm width over the surface of the imaged object. The generated ultrasonic waves are detected by ring-array ultrasound transducer with a center frequency of 5 MHz (60 % bandwidth). The ring-array transducer has 128 elements in the 270-degrees range [44]. The array radius is 40.5 mm. Fig. 2(b) shows the schematic of the ring-array ultrasound transducer setting. PAT images are averaged with signals from 10 frames per wavelength. The ultrasound time-series signals are then reconstructed into 2D PAT images using a model-based iterative reconstruction algorithm [37] with a field of view 30 mm × 30 mm and matrix 300 × 300. Negative values of the reconstructed image are set to zero on the basis that the absorbed energy cannot be negative.

Fig. 2

(a) Diagram of the PAT system. (b) Schematic of the ring-array ultrasound transducer setting. (c) The process of obtaining the un-corrected PAT images used in the simulation experiments. : ideal image. LF map (Φ): estimated light fluence distribution using image. Noise: noise with a mean of 0 and a standard deviation of 2 × 10−5 was added to the un-corrected PAT image. PAT: un-corrected PAT image.

Simulation experiments

In the simulation experiments, we first designed a mouse example to compare the performance between our proposed method and two other comparing methods. Similar to Ref [12], [36], our simulation experiments only incorporate the optical inversion problem that needs to be addressed in this paper, but not the acoustic inversion. We employed a mouse organ model to generate the simulation images, as shown in Fig. 2(c). The value of each organ can be calculated by referring to Ref [45]:where x was sO2, S was a scaling factor relative to whole blood. The values of S and x for different organs are listed in Supplementary Table 1. The values of HbO2 and Hb at different wavelengths were determined by referring to Ref [36]. These values are listed in Supplementary Table 2. The value of each organ was calculated by referring to Ref [45]. These values are listed in Supplementary Table 3. Next, we designed a simulation example that includes vascular network with varying diameters (minimum: 200 µm). The vascular network contains arteries and veins, and the sO2 values of the arteries and veins are set to 100 % and 75 %, respectively. The scattering coefficient and scattering anisotropy values of HbO2 and Hb at different wavelengths were also derived from Ref [36]. The calculated values are listed in Supplementary Table 2.

Fig. 2(c) presents the process of obtaining the un-corrected PAT images. We used Toast+ + software to obtain the LF distribution Φ. As a result, the absorbed energy is proportional to the product of the and the LF intensity within a given voxel. We multiplied the ideal () images by the LF images Φ to obtain the un-corrected PAT images. To further approximate the un-corrected PAT image, we added noise with a mean of 0 and a standard deviation of 2 × 10−5 to the un-corrected PAT images. According to the above procedure, multispectral PAT images were acquired at five different illumination wavelengths: 700, 730, 760, 800 and 850 nm. In the later reconstruction procedure, the map was initialized to 0 and the values were fixed to the ideal value. Then, a linear spectral un-mixing algorithm was performed to decompose the distribution of the absorbers (HbO2 and Hb) from the , PAT and images. Lastly, a high-pass filter was used to filter out the HbO2 and Hb signals in the background.

Phantom experiments

We constructed a tissue-mimicking phantom containing rod-shaped inclusions to further evaluate our proposed method. The standard and method for making the phantom are based on Ref [33], [46]. The background of the phantom was an agar solution mixed with 0.5% Intralipid designed to enhance scattering. Before its solidification, three 3D printed rod molds were inserted into the background solution. After the solution cooled and solidified, the molds were pulled out and three rod-shaped absorbers containing agar solution mixed with 0.5 % Intralipid and 0.02 % Chinese ink were inserted. Referring to Ref [46], the and values were determined to be 0.12 mm−1, 0.5 mm−1 for the absorbers and 0.006 mm−1, 0.5 mm−1 for the background solution at 700 nm, respectively. The diameter of the phantom is 20 mm and the diameter of each inclusion is 3 mm. In addition, we made a phantom containing three absorbers with different values. The values of the three absorbers are 0.05 mm−1, 0.1 mm−1, and 0.15 mm−1, respectively. PAT images at 700 nm were acquired and the speed of sound was set to 1500 m/s. In the later reconstruction procedure, the map was initialized to 0 and the value was fixed to 0.5 mm−1.

Animal experiments

In in vivo animal imaging experiments, four healthy nude mice (female, 8 weeks, Southern Medical University, Guangzhou, China) were used. Animals were kept in ventilated cages inside a temperature-controlled room, under a 12-h dark/light cycle. To reduce abdominal peristalsis artifacts caused by food digestion and to prevent the mice from excreting and polluting the imaging environment during PAT imaging, the nude mice have fasted for 8 h before imaging. All animal experiments were approved by the local Animal Ethics Committee of Southern Medical University and were performed under current guidelines.

To reduce the image artifacts caused by respiratory movements, medical oxygen mixed with 1 % isoflurane (RWD Life Science Co., Ltd, China) was transmitted to the breathing mask, so that the respiratory rate of nude mice was maintained at 15–20 breath/min. For PAT imaging without contrast enhancement, multispectral PAT images of kidney, liver, head, neck and other positions were acquired at five different illumination wavelengths: 700, 730, 760, 800 and 850 nm. For contrast-enhanced PAT imaging, an insulin injection needle was embedded into the tail vein in advance, and was connected to a long Polyethylene Tubing 10 (PE 10) that enabled probe injection (ICG) from outside the imaging chamber. Multispectral PAT images of kidney, liver and other positions were acquired at seven different illumination wavelengths: 740, 760, 780, 800, 820, 840 and 860 nm. For the reconstruction, the speed of sound was set to 1536 m/s. In the later reconstruction procedure, the map was initialized to 0 and the value was fixed by referring to Supplementary Table 4. Finally, the linear un-mixing algorithm was used to calculate the distribution of HbO2, Hb.

Image quality metrics

The quantitative performance metrics used to evaluate the performance of the correction algorithm include L2-norm PAT image error (Err) and the sum of squares error of the images (SSE). Here, the product of the map and the LF distribution Φ is defined as the estimated PAT image (ePAT) to approximate the un-corrected PAT image.

The L2-norm PAT image error is a metric to evaluate the convergence performance of the correction algorithm by calculating the error between the un-corrected PAT image and the ePAT image as:where, PAT is the un-corrected PAT image and ePAT is the estimated PAT image. N is the total pixel number of the image.

The SSE is used to evaluate the error between the image and the reconstructed images in the simulation experiments. It is defined as:where, and are the ideal image and the reconstructed image, respectively.

Results

Simulation Results

The results of the simulation experiment are shown in Fig. 3. In the images, the same concentration of absorber has a uniform signal intensity at each organ, as shown in Fig. 3(a). However, due to the attenuation of laser energy in the transmission process, the PAT signal weakens with the increase of imaging depth, as shown in Fig. 3(b). The prior images of SBDC and SBIC methods are shown in Fig. 3(c), which is generated by manually segmenting the object contours from the PAT image. The LF distribution maps estimated using different methods are shown in Fig. 3(d). The ePAT images are shown in Fig. 3(e) and the reconstructed images are shown in Fig. 3(f). In the SBDC method, we choose 0.05832 mm−1 (average of all organs) as the value. In the SBIC method, the optimal value calculated by iteration is 0.01911 mm−1.

Fig. 3

Simulation results: (a) : ideal distribution image at 850 nm. (b) PAT: un-corrected PAT image at 850 nm obtained by multiplying the image with the LF map and adding ~40 dB of noise. (c) Segmentation prior: segmentation results for SBDC and SBIC methods. (d) LF map (Φ): light fluence distribution map estimated using different methods. (e) ePAT: estimated PAT image derived by using different methods. (f) : image solved by different methods. (g) Difference: the difference images between PAT and ePAT images. (h) Difference: the difference images between and images. (i) Profiles of and images drawn along the white solid line in (a). (j) Profiles of PAT and ePAT images drawn along the white solid line in (a). (k) The Err values between PAT and ePAT images for all positions. (l) The SSE values between and images for all positions. Description of markers: Sk: skin; Int: intestines; K: kidney; Sp: spine; M: muscle; S: spleen. Scale bar, 3 mm.

As can be seen from the result of the SBDC method, the image shows excessive enhancement of central features (such as muscles and kidneys) [Fig. 3(f)]. In the image obtained from the SBIC method, the degree of enhancement is appropriate [Fig. 3(f)], but the signal within the same organ of the image is still inhomogeneous (e.g., spleen). The image obtained by our proposed method eliminates the phenomenon of a lower signal at larger depths, and the signal in each organ has reached an even distribution [Fig. 3(f)]. The reason for this improvement can be found in the LF map and the estimated PAT images shown in Fig. 3(d) and Fig. 3(e) respectively. The LF map estimated by our method has a rapid decay in the spleen, which is caused by the higher value. However, the LF map obtained by SBDC and SBIC methods can not reflect this phenomenon. The ePAT image obtained by the SBIC method does not approximate the un-corrected PAT image because of the assumption of uniform . In contrast, using our method, the ePAT image is very close to the un-corrected PAT image. This proves the accuracy of our proposed iterative algorithm for solving the nonlinear model.

Next, we have calculated the difference between PAT and ePAT, and images. The results are shown in Fig. 3(g, h). It can be observed that the difference of the SBDC method is the largest, whereas the SBIC method only has a greater difference at organs with higher values (e.g., spleen), and the difference gradually increases with increasing imaging depth. In contrast, the difference in our proposed method is so small that it cannot be observed.

For further comparison, Fig. 3(i, j) shows the image profiles of the , , PAT and ePAT. The position of the profiles is along the white solid line in Fig. 3(a). The profiles of the image show that there are errors in the SBDC and SBIC methods compared to the image. These problems have been addressed by the image obtained by our proposed method, as evidenced by the highly overlapping profiles between the and reconstructed image [Fig. 3(i)]. Fig. 3(j) also shows a high degree of overlap between the un-corrected PAT and ePAT images obtained by our proposed method, which demonstrates the high reliability of our model solutions. Moreover, the Err between the entire PAT and ePAT images, the SSE between the entire and images for all positions are shown in Fig. 3(k, l). For both Err and SSE, our proposed method possesses the smallest value (two orders of magnitude smaller than the SBIC method).

Next, we acquired HbO2 and Hb distribution images and further calculated the sO2, as shown in Fig. 4(a-e). As can be seen from the un-corrected results, the light energy attenuation process affects the spectral un-mixing, resulting in the concentration decrease with the increase of depth in the HbO2 and Hb images [Fig. 4(b)]. This further leads to varying sO2 values within the same organ (contrary to the ideal value). Greater error is observed at the organs with a higher absorption coefficient (e.g., spleen) and deeper spatial distribution (e.g. kidney) [Fig. 4(b)]. The SBDC and SBIC methods enhance the concentrations of HbO2 and Hb in the depth of the image [Fig. 4(c, d)], but failed to accurately compensate for the LF in the kidneys, muscle, spine and spleen. There are still quantitative errors in the sO2 values of these regions. Our proposed method brings the concentrations of HbO2 and Hb in each organ to a uniform distribution, and the sO2 values are consistent with the ideal values [Fig. 4(e)].

Fig. 4

Spectral un-mixing results of simulation experiment: (a-e) The concentration distribution (HbO2, Hb) and sO2 images obtained by spectral un-mixing from , un-corrected PAT and images. (f) Profiles of HbO2, Hb and sO2 drawn along the white solid line in (a). Description of markers: Sk: skin; Int: intestines; K: kidney; Sp: spine; M: muscle; S: spleen. Scale bar, 3 mm.

The image profiles of HbO2, Hb and sO2 are illustrated in Fig. 4(f). The position of the profiles is along the white solid line in Fig. 4(a). The HbO2, Hb and sO2 of un-corrected are lower than ideal values, and they are unevenly distributed in the skin, spleen, and kidneys. The results obtained by the SBDC method are higher than the ideal values. The SBIC method restores the HbO2 and Hb concentration estimates of the left skin, intestine and kidney, but the quantitative estimation error of the right organs is still large. In all profiles of HbO2, Hb and sO2, the reconstructed by our proposed method coincides with the . The SSE values between the ideal and estimated HbO2, Hb and sO2 images for all positions are listed in Table 1. For all images, our proposed method possesses the smallest value. The above results show that the HbO2 and Hb concentrations of all organs have been restored to ideal values with our proposed method, as well as sO2.

Table 1

The SSE values between the ideal and estimated HbO2, Hb and sO2 images for all positions.

SSE	Methods
	SBDC	SBIC	Proposed
HbO2	180.351 ± 136.0068	17.2596 ± 35.3555	0.0028 ± 0.001
Hb	6.8366 ± 5.1225	1.1563 ± 2.4457	0.001 ± 0.0003
sO2	10.6515 ± 8.7329	1.0628 ± 0.9624	0.0316 ± 0.0228

Phantom experiment results

We compared our proposed method and two other comparing methods in the tissue-mimicking phantom imaging experiment, and the results at 700 nm are shown in Fig. 5. The three rod-shaped inclusions of the phantom contained the same absorption material in the same concentration [Fig. 5(a)]. In the un-corrected PAT image, the signal of the inner rod is lower than that of the two outer rods due to laser attenuation [Fig. 5(a)]. The prior images of SBDC and SBIC methods are shown in Fig. 5(b). It is generated by manually segmenting the object contours from the PAT image. The LF distribution maps estimated using different methods are shown in Fig. 5(c). The ePAT images are shown in Fig. 5(d). The difference images between PAT and ePAT are shown in Fig. 5(e) and the reconstructed images are shown in Fig. 5(f). In the SBDC method, we choose 0.04 mm−1 as the value. In the SBIC method, the optimal value calculated by iteration is 0.01715 mm−1.

Fig. 5

Tissue-mimicking phantom experiment: (a) PAT: un-corrected PAT image. (b) Segmentation prior: segmentation results for SBDC and SBIC methods. (c) LF map (Φ): light fluence distribution maps estimated using different methods. (d) ePAT: estimated PAT images derived by using different methods. (e) Difference: the difference images between PAT and ePAT images. (f) : images solved by different methods. (g) Profiles of PAT and ePAT images drawn along the white line in (a). (h) Profiles of images drawn along the white line in (a). Scale bar, 3 mm.

At 700 nm, the value is relatively high in Chinese ink (0.12 mm−1), moderate in Intralipid (0.006 mm−1), and very small in water [Fig. 5(a)]. As expected, the LF map derived from our method shows that the LF does not change obviously in water, gradually decreases in the background containing Intralipid, and rapidly decays in the rod-shaped absorbers containing Chinese ink [Fig. 5(c)]. However, this feature has not been reflected in SBDC and SBIC methods. Similar to the un-corrected PAT image, the ePAT image derived from all three methods has an attenuated signal with increasing imaging depth [Fig. 5(d)]. However, the ePAT image derived from our method has a visually more similar attenuation variation to PAT: the signal intensity of the inner rod is lower than that of the two outer rods [Fig. 5(d)]. Further, it can be observed that the difference in SBDC and SBIC methods is larger than our proposed method, especially in regions with higher values (e.g., rod-shaped inclusions) [Fig. 5(e)]. In the reconstructed image obtained by our proposed method, the three rod-shaped absorbers have very similar pixel values, and for the same rod-shaped absorber, the pixel values become more homogeneous [Fig. 5(f)].

For further comparison, the image profiles of PAT, ePAT and are illustrated in Fig. 5(g, h), with the position along the white solid line in Fig. 5(a). The profiles show a high degree of overlap between the PAT and ePAT obtained by our proposed method, and the PAT signal of the two outer rod-shaped absorbers decreases with the increase of depth [Fig. 5(g)]. By contrast, there are large errors between the profiles of PAT and ePAT obtained by the SBDC and SBIC methods. The profile of the image obtained by our method shows that the values of these three rod-shaped inclusions have reached the same level, and the values of the two outer absorbers have become more homogeneous [Fig. 5(h)]. These results suggested that our method can reconstruct the real distribution in the tissue-mimicking phantom experiment.

Animal results

The results of the in vivo imaging of nude mice at 700 nm are shown in Fig. 6. Fig. 6(a) illustrates the PAT signal attenuation in the central features [artery, portal veins (PVs), inferior vena cava (IVC), the interior of the kidney and spine] due to the attenuation of the laser energy. The prior images of SBDC and SBIC methods generated by manually segmenting the object contours from the PAT image are shown in Fig. 6(b). The LF distribution maps estimated using different methods are shown in Fig. 6(c). The ePAT images are shown in Fig. 6(d) and the difference images between PAT and ePAT are shown in Fig. 6(e). The reconstructed images are shown in Fig. 6(f).

Fig. 6

In vivo animal experiment: (a) Un-corrected PAT image at the kidney position. (b) Segmentation prior: segmentation results for SBDC and SBIC methods. (c) LF map (Φ): light fluence distribution maps estimated using different methods. (d) ePAT: estimated PAT images derived by using different methods. (e) Difference: the difference images between PAT and ePAT images. (f) : images solved by different methods. (g) Profiles of PAT and ePAT images drawn along the white line in (a). (h) Profiles of images drawn along the white line in (a). (i) The Err values between PAT and ePAT images for all nude mice at different wavelengths. Description of markers: A: artery; PV: portal vein; IVC: inferior vena cava; K: kidney; Sp: spine. Scale bar, 3 mm.

As can be seen from the results of SBDC and SBIC methods, the LF map and the ePAT image also have an attenuated signal with the increasing imaging depth [Fig. 6(c, d)]. However, this isotropic attenuation resulting from the uniform value makes the difference between PAT and ePAT very large [Fig. 6(e)]. The LF map estimated by our proposed method has a rapid decay in regions with higher values (kidneys, PV and spine), whereas the LF maps obtained by SBDC and SBIC methods attenuate isotropically [Fig. 6(c)]. The difference between PAT and ePAT obtained by our method is the smallest [Fig. 6(e)], which proves that one of the advantages of our method over SBDC and SBIC methods is that it can obtain a more accurate solution.

In the SBDC method, the accuracy of the results is depended on the manual choice of the value. For example, the image obtained by will cause excessive enhancement of central features (such as artery, IVC, and PVs) [Fig. 6(f)]. The image obtained by the SBIC method appears to be appropriate for the enhancement of the signal from deeper tissues (), but the signal within the same organ of the image is still inhomogeneous (e.g. kidney) [Fig. 6(f)]. The image obtained by our method effectively enhances the visibility of deep tissues, e.g., image intensity in the kidney has reached an even distribution [Fig. 6(f)].

To further assess the performance of these three methods, the image profiles of PAT, ePAT and at the position along the white solid line in Fig. 6(a) are illustrated in Fig. 6(g, h). The profiles show a high degree of overlap between the PAT and the ePAT image obtained by our proposed method [Fig. 6(g)]. By contrast, there are large errors between the profiles of PAT and ePAT obtained by the SBDC and SBIC methods. In the profile of the PAT image, the signals in the superficial layers of both the left and right kidneys are higher than those in deeper layers [Fig. 6(g), K arrow]. In the SBDC method, the signals of the deeper layers of the kidney are enhanced indiscriminately, resulting in much higher internal signals than external ones [Fig. 6(h), K arrow]. In the SBIC method, there is a small difference in the profile amplitudes within the left and right kidneys [Fig. 6(h), K arrow], which cannot be addressed by the setting of uniform value. In our proposed method, there is no obvious difference in the profile amplitudes within the left and right kidneys [Fig. 6(h)], and the signal amplitude of the central IVC is consistent with that of the SBIC method.

Furthermore, we measured the Err between the un-corrected PAT and ePAT images in the in vivo animal experiment of all four mice at different wavelengths, and the results are shown in Fig. 6(i). For the Err results, our proposed method is on average two orders of magnitude smaller than SBDC and SBIC methods at all wavelengths. Specifically, at 760 nm, the mean Err value of our proposed method for all animals is 0.091 ± 0.13, compared to 2.8906 ± 0.2711 for the SBDC method and 1.6342 ± 0.4762 for the SBIC method. This indicates that the ePAT image estimated by our method has reached a better approximation of the un-corrected PAT image.

Fig. 7(a-d) shows the HbO2, Hb and sO2 distribution images obtained from linear spectral unmixing. The HbO2 of oxygen-rich (e.g., the artery) and Hb of hypoxic (the PV and IVC) blood vessels should be higher than in other tissues. However, these features are not reflected in the un-corrected images due to the laser energy attenuation [Fig. 7(a)]. This further leads to increased sO2 value with increasing depth. The SBDC and SBIC methods enhance the concentrations of HbO2 and Hb in the depth of the image [Fig. 7(b, c)]. However, due to the uniform value setting under different wavelengths in SBDC and SBIC methods, the sO2 value still tends to increase with depth. There are still quantitative errors in the sO2 values of the artery [Fig. 7(b, c)]. In the HbO2 and Hb images obtained by our proposed method, all blood vessels are visible [Fig. 7(d)].

Fig. 7

Spectral un-mixing results of in vivo animal experiment: (a-d) The concentration distribution (HbO2, Hb) and SO2 images obtained by spectral un-mixing from un-corrected PAT and images. (e-g) The difference images of HbO2, Hb and sO2 by using un-corrected images as reference. Description of markers: A: artery; PV: portal vein; IVC: inferior vena cava; K: kidney; Sp: spine. Scale bar, 3 mm.

We have also calculated the difference images of HbO2, Hb and sO2 by using the un-corrected images as reference, and the results are shown in Fig. 7(e-g). As can be seen, the HbO2 concentration in the artery and the Hb concentration in the veins (PV and IVC) have increased in the results. Nevertheless, the sO2 value at the artery obtained by our proposed method is still high, while the sO2 value at the kidney and PV is significantly reduced. This indicates that our proposed method can obtain a more accurate endogenous absorber concentration distribution.

Finally, we applied these three methods to multispectral PAT images of the kidney region acquired from living nude mice after intravenous administration of 200 μl ICG (0.05 μg/μl). The spectrally un-mixed ICG images are shown in Fig. 8(a-d), and the difference images of ICG by using the un-corrected image as reference are shown in Fig. 8(e-g). the ICG is metabolized from the veins into the kidneys and spleen, and thus these positions should have high ICG signals. However, in the un-corrected spectral un-mixing results, the ICG concentration in the kidneys is not homogeneous (inside is lower than the outside), and the ICG signal of the IVC is very weak, which seriously affected the quantitative study of ICG metabolism [Fig. 8(a)]. The center concentration of ICG obtained by the SBDC method is over-enhanced [Fig. 8(b, e)], while the SBIC method has no significant improvement [Fig. 8(c, f)]. Both these two methods cannot restore the higher ICG concentration distribution in the spleen due to metabolism. The result derived from our proposed method shows that the visualization of deep ICG concentration, including the IVC, kidney and spleen, is enhanced [Fig. 8(d, g)]. These results indicate that our method can help obtain a more accurate exogenous probe distribution.

Fig. 8

ICG experiment: (a-d) The concentration distribution images of ICG obtained by spectral un-mixing. (e-g) The difference images by using the un-corrected image as reference. Description of markers: IVC: inferior vena cava; K: kidney; S: spleen. Scale bar, 3 mm.

Discussions

The above results demonstrate that our method can directly reconstruct pixel-wise tissue maps and obtain the correct absorber distributions, revealing its attractive potential for the quantitative measurement of blood oxygenation and the concentration of various exogenous contrast agents. Our proposed method includes the following advantages:

Firstly, to ensure accuracy, our method converts the reconstruction problem to a nonlinear least square model, and proposes a two-step iterative optimization method to update the LF distribution and the map until convergence. This iterative algorithm allows optimizing both the LF map and the image simultaneously, and thus improves the convergence of the solution.

Secondly, previous methods such as SBDC and SBIC consider uniform distributions either on the whole body or within each organ. They do not match the real, spatially varying optical parameters in tissues. By contrast, our method directly solves this inhomogeneous map at pixel-level. By solving it in an iterative way, the obtained map is more accurate than region-wise methods, as proved by our experiments.

Thirdly, our proposed method is simple and does not require manual intervention. Both the SBDC and SBIC methods require segmentation of PAT images, and the SBDC method even requires manual assignment of optical parameters. Furthermore, the performance of the SBDC method is highly dependent on manually selected value, choosing a different value can lead to a deficient or excessive correction. Our method can automatically calculate the distribution without manual segmentation of the animal body or organs, thus largely simplifies the correction process and avoids the influence of subjective factors.

There are still some limitations in our study. Firstly, as in most previous studies [12], [33], [40], [41], we assume that the reduced scattering coefficient is known because it is difficult to solve both the absorption and reduced scattering coefficients at the same time [32]. Secondly, we use a two-dimensional light transport model to estimate the LF distribution as in Ref [32]. However, in practice, the LF estimation should be carried out in three-dimension because PAT systems typically use diffuse illumination. Future work will require the design of the special objective functions for the 3D optical transport model for LF estimation, which in turn enables 3D reconstruction. Thirdly, our method requires a longer computation time while improving computational accuracy. The average computation time of our proposed method for a single-slice image is 92.49 s, so real-time imaging cannot be achieved. Finally, during image reconstruction, we set a uniform speed of sound distribution, but biological tissues are acoustically heterogeneous, therefore in our future work we also need to correct for the heterogeneous speed of sound.

Conclusions

In this paper, we developed a non-segmentation iterative method to directly reconstruct the images for quantitative PAT. Demonstrated in the simulation experiments, phantom experiments and in vivo animal experiments, our method enables accurate reconstruction of spatially varying images at pixel-level. These findings provide a reliable foundation for the quantitative measurements of blood oxygenation and the various exogenous contrast agents, and therefore contribute to the advancement and future applications of the PAT technology.

CRediT authorship contribution statement

Shuangyang Zhang and Jiaming Liu: Methodology, Software, Investigation, Data Analysis, Writing - Original Draft Zhichao Liang and Jia Ge: Resources, Animal & Phantom Experiments Yanqiu Feng, Wufan Chen, and Li Qi: Conceptualization, Supervision, Funding acquisition, Writing - Review & Editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table

Algorithm 1 Two-Step Iterative Algorithm