Enhancement of photoacoustic tomography of acoustically inhomogeneous tissue by utilizing a memory effect

Jie Yin; Jie Yin; Junjun He; Chao Tao; Chao Tao; Xiaojun Liu; Xiaojun Liu

doi:10.1364/OE.388902

1. Introduction

Photoacoustic tomography (PAT) is an emerging modality that combines rich optical contrasts with good ultrasound resolution in deep tissue. Under the illumination of pulse laser, the absorbers in tissue will generate ultrasound signals, i.e., photoacoustic (PA) signals, due to thermoelastic expansion after absorbing laser energy. If the PA signals are recorded by ultrasound transducers arranged around the tissue, an image of the optical absorption distribution can be decoded from the detected PA signals by using a reconstruction algorithm [1–6]. In recent years, PAT has shown great potential for tumor detection [7,8], osteoarthritis assessment [9], drug delivery monitoring [10], small-animal whole body panoramic imaging [11], surgery guidance [12], identification of atherosclerosis plaque [13,14] and so on [15–18].

Reconstruction algorithms play key roles in PAT. Most PAT reconstruction algorithms are based on the principle of coherent beamforming, which usually assumes a homogeneous distribution of speed of sound (SOS) in tissue. However, in the most realistic tissues, the distribution of SOS is inhomogeneous due to acoustic properties discrepancy of various tissue. For soft tissue, the SOS changes from 1400 m/s (e.g., fat) to 1600 m/s (e.g., muscle), depending on the compositions of the tissue [19]. Inhomogeneous distribution of SOS would lead to poor image contrast and resolution [20–22].

A lot of work have been done to improve the performance of PAT for acoustically inhomogeneous tissue. The autofocus approach is utilized to automatically select SOS by maximizing the sharpness of the reconstructed images [23]. The coherent factor method is developed to estimate the average SOS [24]. The full-wave iteration method is proposed to mitigate the effects of SOS information incompleteness and noise [25]. The photoacoustic and ultrasound tomography measurements are combined to acquire the SOS distribution and compensate the photoacoustic imaging distortion induced by SOS deviation [26,27]. The time reversal method based on Green’s function’s retrieval is also established to improve the PA image quality in scattering medium [28]. Based on Eikonal equation, which uses ray refracted paths instead of straight ray paths, Jose et al. proposed an iterative approach for reconstruction of SOS map with which an accurate PA image can be achieved [29]. Several methods for automatic calibration of a uniform SOS were investigated experimentally by Mandal et al., and two of them were found superior to the others [30]. Cai et al. established a feature coupling method to reconstruct initial pressure and SOS jointly in vivo, which iteratively optimizes a SOS distribution by maximizing the similarity of images from the subsets of a full transducer array [31]. However, imaging in acoustically inhomogeneous tissue is still a challenge for PAT [32–34].

It was found that the incident light would give rise to a transmitted speckle when it propagates through a thin inhomogeneous medium. As the incident wave is rotated by a small degree, a new transmitted speckle pattern can be acquired. The new pattern can be deemed as a replica of the pattern against the original, but was shifted a few pixels along a specific direction. That is referred to as “memory effect” [35,36].

In this work, we prove that a similar effect can be observed as PA signals propagate through a medium with inhomogeneous distribution of SOS. The PA signal received by one ultrasound transducer has a memory of the waveform of PA signals received by adjacent transducers. Taking into account this effect, we develop a method to estimate an optimal SOS for each transducer in PAT. The optimized SOS will eventually improve the PAT image quality for inhomogeneous tissue. Both numerical simulations and phantom experiments have been employed to validate the proposed method.

2. Theory

Considering the scenario in Fig. 1, the yellow disk is the imaging target we set, four black dots are optical absorbers, and a ring-shape ultrasound transducer array is arranged at a closed curve surrounding the sample to detect the PA signals.

Fig. 1. (a) The schematic of the scenario used for theoretical derivation. Four black dots are optical absorbers, and the medium consists of two regions, the SOS in region I and II are 1460 m/s, and 1540 m/s, respectively. (b) PA signals received by transducers at π/2, π, 3π/2 and 2π (blue solid lines), and the adjacent transducers (red dash lines).

Download Full Size | PDF

A pulse laser beam illuminates the target, then the PA signals stimulated by absorbers are recorded by the transducer array. The PA signal recorded by a transducer at r_d can be expressed by Eq. (1),

(1)$$p({{{\textbf r}_d},t} )= \frac{1}{{4\pi {c^2}}}\int_V {\frac{{A({\textbf r}){{dH(t - |{{\textbf r} - {{\textbf r}_d}} |/c)} \mathord{\left/ {\vphantom {{dH(t - |{{\textbf r} - {{\textbf r}_d}} |/c)} {dt}}} \right.} {dt}}}}{{|{{\textbf r} - {{\textbf r}_d}} |}}dV} ,$$

where c is the SOS in soft tissues, which is generally assumed to a constant of 1500 m/s, A(r) is the spatial optical absorption function, H(t) is the temporal profile of the laser beam. Since the laser pulse is extremely narrow, H(t) can be approximately considered as the Dirac Delta function. Therefore p(r_d, t) can be obtained by integrating through the path S,

(2)$$p({{{\textbf r}_d},t} )= \frac{1}{{4\pi {c^2}}}\frac{\textrm{d}}{{\textrm{d}t}}\int_S {\frac{{A({\textbf r})}}{{|{{\textbf r} - {{\textbf r}_d}} |}}\textrm{d}s} ,$$

where S represents the path of integration |r – r_d| = ct, in 2-D applications, r = (x, y), and |r – r_d| = ct can be written as (x – x_d)² – (y – y_d)²= c²t². Then we have x = r_dcosθ + ctcosφ and y = r_dsinθ + ctsinφ, where r_d is the scanning radius, x_d = r_dcosθ and y_d = r_dsinθ. Then substitute them into Eq. (2), in order to simplify the issues, coefficient 1/4πc² which causes minor effect to the waveform is ignored. For a given t (R = ct), it has,

(3)$$\begin{array}{l} p({{{\textbf r}_d},t} )= \frac{d}{{dt}}\int_{{\varphi _1}}^{{\varphi _2}} {A({x,y} )d\varphi } \\ \textrm{ = }\frac{d}{{dt}}\int_{{\varphi _1}}^{{\varphi _2}} {A({{r_d}\cos \theta + R\cos \varphi ,{r_d}\sin \theta + R\sin \varphi } )d\varphi } . \end{array}$$

For a transducer which is located at r_d’, as shown in Fig. 1, the PA signal it received can be expressed as,

(4)$$p({{{\textbf r}_d}^{\prime},t} )= \frac{d}{{dt}}\int_{{\varphi _1}^{\prime}}^{{\varphi _2}^{\prime}} {A({{r_d}\cos \theta^{\prime} + R^{\prime}\cos \varphi^{\prime},{r_d}\sin \theta^{\prime} + R^{\prime}\sin \varphi^{\prime}} )d\varphi ^{\prime}} ,$$

where θ’ = θ+ Δθ, φ₁’ =φ₁ + Δθ, φ₂’ =φ₂ + Δθ and φ’ =φ + Δθ. As the detecting position changes, the propagating path of the PA signal also changes along with it, consequently the average SOS changes. For a given time t [t is the same as in Eq. (3)], R’ = c't = (c +Δc)t, where Δc is the increment of SOS induced by the changing of the propagating path. Then reorganized the above equations by using the trigonometric formula, Eq. (4) can be written as,

(5)$$p({{{\textbf r}_d}^{\prime},t} )= \frac{d}{{dt}}\int_{{\varphi _1}}^{{\varphi _2}} {A\left( \begin{array}{l} ({r_d}\cos \theta + R^{\prime}\cos \varphi )\cos \Delta \theta - ({r_d}\sin \theta + R^{\prime}\sin \varphi )\sin \Delta \theta ),\\ ({r_d}\sin \theta + R^{\prime}\sin \varphi )\cos \Delta \theta + ({r_d}\cos \theta + R^{\prime}\cos \varphi )\sin \Delta \theta ) \end{array} \right)d\varphi } .$$

For Δθ is very small, we can expand cosΔθ and sinΔθ by using the Taylor series, then ignore the higher order and keep the first order, we have,

(6)$$p({{{\textbf r}_d}^{\prime},t} )= \frac{d}{{dt}}\int_{{\varphi _1}}^{{\varphi _2}} {A({x - y\Delta \theta + \cos \varphi^{\prime}\Delta ct,y + x\Delta \theta + \sin \varphi^{\prime}\Delta ct} )d\varphi } ,$$

compared with Eq. (3), then Eq. (6) can be rewritten as,

(7)$$p({{{\textbf r}_d}^{\prime},t} )= \frac{d}{{dt}}\int_{{\varphi _1}}^{{\varphi _2}} {A({{\textbf r} + \Delta {{\textbf r}_\theta } + \Delta {{\textbf r}_c}} )d\varphi } ,$$

where Δr_θ = (–yΔθ, xΔθ) and Δr_c = (Δctcos(φ + Δθ), Δctsin(φ + Δθ)) are disturbing terms introduced by detecting position offset and SOS variance respectively.

Now two conclusions can be deduced from the aforementioned theoretically analysis. Firstly, the PA signals received at two adjacent positions are highly correlative, and the correlation between them is affected by two factors, one is the detecting position offset, the other is the SOS variance of the tissue. Secondly, the correlation of two PA signals detected at two adjacent positions contains the information of SOS variance. For |Δr_θ| = (x² + y²)^1/2Δθ, |Δr_c| = Δct = ct(Δc/c) = R(Δc/c), if r_d and r_d’ are close enough, i.e., Δθ is very small, and the radius r₀ of imaging region is much smaller than the array radius r_d, i.e., r₀ << r_d or r₀ << R, we can learn that, |Δr_θ| << |Δr_c|. In other words, the variance of PA signals induced by detecting position offset is much smaller than that induced by SOS variance. Therefore, by ignoring the small term Δr_θ in Eq. (7), we have,

(8)$$p({{{\textbf r}_d}^{\prime},t} )\approx \frac{d}{{dt}}\int_{{\varphi _1}}^{{\varphi _2}} {A({{\textbf r} + \Delta {{\textbf r}_c}} )d\varphi } ,$$

for the sake of contrast, Eq. (8) is rearranged as,

(9)$$p({{{\textbf r}_d}^{\prime},t} )\approx \frac{d}{{dt}}\int_{{\varphi _1}}^{{\varphi _2}} {A({{r_d}\cos \theta + c(t + \tau )\cos \varphi ,{r_d}\sin \theta + c(t + \tau )\sin \varphi } )d\varphi } ,$$

that is,

(10)$$p({{{\textbf r}_d}^{\prime},t} )\approx p({{{\textbf r}_d},t + \tau } ),$$

where τ = (Δc/c)t. Equation (10) shows that, if the scanning step Δθ of PAT is small enough, and the size of imaging target is much smaller than the scanning radius, the PA signal received by one transducer has a “memory” of the PA signals received by adjacent transducers, which contains the information of SOS variance.

According to the theoretical derivation, we can know that if a set of accurate time delays were applied to the received PA signals, the correlation coefficients between PA signals received by two adjacent transducers would reach the maximum value, and the optimal SOS of each transducer can be determined by τ = (Δc/c)t. An iterative algorithm for optimal SOS estimation is proposed. For biomedical applications, we set $c_0^1 = c_0^2 = \ldots c_0^n =$ 1500 m/s as the initial values, where the subscripts denote the times of iterations, while the superscripts denote the number of transducers, then we set Δc as a variation range of SOS, i.e., Δc = [ − 10 : 0.05 : 10] m/s, the correlation coefficient is utilized to find the optimal Δc, i.e., Δc_opt. For the n-th and (n − 1)-th signals, we have,

(11)$$\rho _{^m}^n = \frac{{Cov\left( {p\left( {{{\textbf r}_d}^n,t - \frac{{{r_d}}}{{c_m^n + \Delta c}}} \right),p\left( {{{\textbf r}_d}^{n - 1},t - \frac{{{r_d}}}{{c_m^{n - 1} - \Delta c}}} \right)} \right)}}{{\sqrt {Var\left( {p\left( {{{\textbf r}_d}^n,t - \frac{{{r_d}}}{{c_m^n + \Delta c}}} \right)} \right)Var\left( {p\left( {{{\textbf r}_d}^{n - 1},t - \frac{{{r_d}}}{{c_m^{n - 1} - \Delta c}}} \right)} \right)} }},$$

(12)$$\Delta c_{opt}^n = \arg \max ({\rho_{^m}^n} ),$$

where Cov denotes covariance, Var denotes variance. We apply Eq. (11) to each transducer in turn, and a set of Δc_opt can be determined, then we set $c_{m + 1}^n = c_m^n + \Delta c_{opt}^n$, $c_{m + 1}^{n - 1} = c_m^{n - 1} - \Delta c_{opt}^n$, and take them into Eq. (11), then repeat the process. For practical applications, we use speed error, ${C_{err}}\textrm{ = (1/}n\textrm{)}\sum\nolimits_n {|{c_m^n - c_{m - 200}^n} |}$ with m > 200 to examine the convergence of the computation. When the computation is converged, the sum of SOS $\sum\nolimits_n {|{c_m^n} |}$ is basically stable and C_err should keep within a reasonable range. In this study, we set C_err < 2.5×10⁻⁴ mm/µs as the criteria to stop the iteration. As the iteration is completed, the optimal SOS, c_opt for each transducer can be determined.

3. Numerical simulations

Numerical experiments were first carried out to evaluate the performance of the proposed method. As shown in Fig. 1, the yellow disk is imaging target, which contains four optical absorbers with a diameter of 1 mm, while the centers of the absorbers are located at (−2.5 mm, 2.5 mm), (2.5 mm, 2.5 mm), (−2.5 mm, −2.5 mm) and (2.5 mm, −2.5 mm), respectively. The relative optical energy absorption deposition, i.e., A(r) of four optical absorbers was 1.0, while that of the other region was 0.0. The SOS and density of the target were c = 1.5 mm/µs and ρ = 1000 kg/m³, respectively. The medium in region I had a SOS of c = 1.54 mm/µs and a density of 1000 kg/m³, while the medium in region II had a SOS of c = 1.46 mm/µs and a density of 1000 kg/m³. The ring-shape transducer array consisted of 120 elements, which were arranged on a circle with a diameter of 60 mm.

After the irradiation of a laser beam, PA signals with a center frequency of 2.5 MHz were recorded by the transducers. The PA signals recorded by the transducers at the position of θ = π/2, π, 3π/2 and 2π are shown in Fig. 1(b) with blue solid lines, while the red dash lines are PA signals recorded by adjacent transducers. As we can see, the waveforms of blue lines are similar to that of red lines, but time delays exist between them. This is in coincidence with the theoretical derivation mentioned before.

Before the reconstruction of PA image, the issue on how small of Δθ is suitable for computation needs to be addressed. We set iteration times to 1000, and iterated over the recorded PA signals according to Eqs. (11) and (12), then the calculation SOS, c of each transducer can be got. The calculation results of c with different Δθ are shown in Fig. 2(a), the horizontal axis represents the positions of transducers, while the vertical axis represents the calculation SOS of the transducers. The black dash line is the actual SOS, and the color lines are calculation results with different Δθ. For this numerical model, for we know the actual SOS, we plot the absolute speed errors with different Δθ in Fig. 2(b), the horizontal axis represents the values of Δθ, and the vertical axis represents the absolute speed error between actual and calculation SOS, $C_{err}^{\prime} = {\left. {\frac{1}{n}\sum\limits_{n = 1}^{120} {|{c_m^n - c_0^n} |} } \right|_{m = 1000}}$, where c₀ denotes the actual SOS we set.

Fig. 2. The relationship between Δθ and the accuracy of SOS estimation. (a) The calculation SOS acquired with different Δθ (iteration times m = 1000). (b) Absolute speed errors with different Δθ.

Download Full Size | PDF

As Fig. 2 shows, the accuracy of SOS reconstruction improves as Δθ decreases, however there is no significant difference in calculation values when Δθ is smaller than 4°. Considering that the smaller Δθ is, the higher computational cost is, we set Δθ to 3° in the rest simulations and experiments.

Fig. 3. Simulative results of the model in Fig. 1. (a) The calculation SOS of transducers with different iteration times. (b) The relationship between absolute speed error and iteration times, red line is the result without noise, while the black line is the result with noise. (c) PA image reconstructed by the DAS method with a fixed SOS. (d) PA image reconstructed by the DAS method with an optimized SOS.

Download Full Size | PDF

After the Δθ is determined, the calculation SOS with different iteration times are shown in Fig. 3(a), the horizontal axis represents the positions of transducers, while the vertical axis represents the calculation SOS of the transducers. In order to verify the robustness against noise of the method, a white Gaussian noise was added in the computation, and the signal to noise ratio, i.e. SNR, is 6.6 dB. Figure 3(a) shows that, the computation is more accurate with the increasing of decoding iteration times, meanwhile the computation is insensitive to noise. The relationship of absolute speed error and iteration times is given in Fig. 3(b), where the horizontal axis represents iteration times, and the vertical axis represents the absolute speed error $C_{err}^{\prime} = \frac{1}{n}\sum\limits_{n = 1}^{120} {|{c_m^n - c_0^n} |}$, Fig. 3(b) shows that the speed error decreases as the iteration times increase, and the error tends to be stable when m reaches 1000. The discrepancy between actual and calculation SOS might be caused by the truncation of the higher order. Then, delay-and-sum method (DAS) was carried out to reconstruct PA image, the image reconstructed with a fixed SOS of 1500 m/s is demonstrated in Fig. 3(c), while the image reconstructed with the optimized SOS acquired with iteration times of 2000 [green line in Fig. 3(a)] is shown in Fig. 3(d), where the distortion is restrained as well as the speckle noise, the image quality is improved significantly.

Fig. 4. Numerical model with complex SOS distribution.

Download Full Size | PDF

In order to verify the practicality of the proposed method in biomedical applications, a complex model was introduced, which was made according to a cross-section image of a mouse’s abdomen. The shape of the model was irregular, and the acoustic characteristics of each part in the model were different, the SOS were set according to the animal innards and tissues.

The schematic of the complex model is shown in Fig. 4. Four dark yellow disks with a diameter of 1 mm in Part 1 had an A(r) of 1.0, the rest region in Part 1 had an A(r) of 0.2, and the A(r) of the other parts in the model was zero. The centers of the four disks were located at (−1.5 mm, 1.5 mm), (−1.5 mm, −1.5 mm), (1.5 mm, 1.5 mm) and (1.5 mm, −1.5 mm), respectively. The SOS of Part 1 to 6 was set to 1500 m/s, 1450 m/s, 1540 m/s, 1493 m/s, 1550 m/s and 1560 m/s, respectively. The density of Part 1 to 6 was set to 1000 kg/m³, 950 kg/m³, 1010 kg/m³, 990 kg/m³, 1020 kg/m³ and 1030 kg/m³, respectively. The inhomogeneity is complex and random. The configuration of the transducer array was the same as the model in Fig. 1.

PA signals with a center frequency of 2.5 MHz were recorded by the transducers, then the iterative operation was applied to the signals. The calculation SOS with different iteration times are given in Fig. 5(a). Under the criterion of ${C_{err}}\textrm{ = (1/}n\textrm{)}\sum\nolimits_n {|{c_m^n - c_{m - 200}^n} |}$ < 2.5×10⁻⁴ mm/µs, this iteration was stopped when m reaches about 5000, as shown in Fig. 5(b), and the SOS had been optimized. The PA images reconstructed with a fixed SOS of 1500 m/s and an optimized SOS [red line in Fig. 5(a)] are given in Figs. 5(c) and 5(d), respectively. As the figures show, with the optimized SOS, four disks approaches to their actual shape, and can be distinguished clearly, the image quality is improved in medium with complex SOS distribution.

Fig. 5. Simulative results of the model in Fig. 4. (a) The calculation SOS with different iteration times. (b) The relationship between speed error and iteration times. (c) PA image reconstructed by the DAS method with a fixed SOS. (d) PA image reconstructed by the DAS method with an optimized SOS.

Download Full Size | PDF

4. Experiments

The performance of the method was further investigated by experiments. Figure 6(a) shows the diagram of the experimental system. The irradiation source was a Q-switched Nd:YAG pulse laser with 532 nm wavelength and 8 ns pulse width. The laser beam had a pulse energy of about 80 mJ and a diameter of about 10 mm. PA signals were detected by a cylindrically focused transducer with a central frequency of 5 MHz (V310, Panametrics). Then, signals were amplified (SA-230FS, NF) and recorded at a rate of 60 MHz (PCI-5105, NI). Under the drive of a rotary stepper motor, the transducer scanned the sample with a radius of 38.5 mm and a step of 3°, which is equivalent to an annular array with 120 elements.

Fig. 6. Experimental setup. (a) Sketch map of the experimental system. (b) Photograph of the sample. Part II had a higher SOS than Part I. The distribution of SOS in the whole phantom was inhomogeneous.

Download Full Size | PDF

Figure 6(b) gives a photograph of the phantom, which was made of agar. Five hairs were assigned a vessel-like shape and embedded in the phantom. In order to simulate the variance of the SOS in realistic tissue, the phantom consisted of two parts with different SOS. Part I was a mixture of agar and water with a ratio of about 1.2%. The SOS in Part I was about 1.46 mm/µs [37]. Part II was a mixture of water, agar, and 1-propanol, with a ratio of 100: 1.2: 3. For 1-propanol rose the SOS of Part II to about 1.588 mm/µs [37], Part II had a higher SOS than Part I, the distribution of SOS in the whole phantom was inhomogeneous.

The proposed method was applied to the received PA signals. Figure 7 shows the process of optimizing the SOS by using the proposed method. Figure 7(a) give the waveforms of PA signals captured at four different positions. The black and red lines are signals recorded at two adjacent positions 1 and 2, while the blue and pink lines are recorded at two adjacent positions 3 and 4, as shown in Fig. 6(b). The signals detected at adjacent positions have similar waveforms. However, a tiny time delay exists between the signals detected at position 1 [black line] and 2 [red line]. And there is no obvious time delay between the signals detected at position 3 [blue line] and 4 [pink line]. This is because the SOS along the travel path from the optical absorbers to the position 1 and 2 was different, but the path between the position 3, 4 and the optical absorbers was homogeneous (only Part I). These results agree with the predictions of Eq. (10). The calculation SOS with different iteration times are given in Fig. 7(b). Under the criterion of ${C_{err}}\textrm{ = (1/}n\textrm{)}\sum\nolimits_n {|{c_m^n - c_{m - 200}^n} |}$ < 2.5×10⁻⁴ mm/µs, this iteration was stopped when m reaches about 6000, as shown in Fig. 7(c), and the SOS has been optimized.

Fig. 7. Process of optimizing the SOS. (a) The PA signals received at four different position. (b) The calculation SOS with different iteration times. (c) The relationship between speed error and iteration times.

Download Full Size | PDF

Figure 8 presents the images reconstructed by using DAS method with the preset consistent SOS of 1.5 mm/µs [Fig. 8(a)] and the optimized SOS [Fig. 8(b)], respectively. As shown in Fig. 8(a), the heterogeneity of the medium significantly degrades the image quality. The variance of SOS contributes position deviations and distortions to the image. The shape of the optical absorbers is hard to be distinguished. Then we reconstructed the PA image with the optimized SOS acquired by the proposed method, as illustrated in Fig. 8(b), the distortion is restrained and the false contrast between the absorbers is mitigated, hairs can be distinguished clearly. The image quality has been improved by using the proposed method.

Fig. 8. Reconstructed PA images of the phantom (a) with the preset homogeneous SOS, and (b) with the optimized SOS.

Download Full Size | PDF

5. Summary

In this work, we reveal the “memory effect” in the PA signals propagating through media with inhomogeneous distribution of SOS. We theoretically demonstrate that PA signals received by two adjacent transducers are similar in waveform but different in phase. The phase difference depends on the distributions of SOS along the propagation paths. Based on this physical phenomenon, we propose an iterative method to optimize the SOS used for image reconstruction. Both simulations and experiments verified the proposed method. The results demonstrate that the proposed method is robust against noise and a significant enhancement of PAT has been achieved. In addition, this method does not need to know about the prior knowledge of tissue inhomogeneity, nor a homogeneous detection window.

Although PAT improved by the proposed method, artifacts and distortions still exist and it could be caused by the following reasons. First, the proposed method is utilized to improve the image quality in medium with inhomogeneous SOS distribution, and the axial SOS inhomogeneous is mainly considered, but refraction and reflection exist in both simulations and experiments, which would still bring some noise and distortions in the reconstructions. Second, the truncation of the higher order in theoretical derivation could introduce some minor error in SOS estimation. Finally, for the phantom experiments, the measurement errors, e.g., the misaligned start time of laser fired, could also induced artifacts and distortions. These are the limitations need to be addressed in our following researches.

In summary, although PAT has demonstrated great potential in biomedical applications, it is still facing many challenges, inhomogeneous distribution of SOS in realistic tissue is one of them. The proposed method may have potential value in improving the performance of PAT in the realistic biomedical applications.

Funding

Natural Science Foundation of Jiangsu Province (BK 20181077); The Natural Science Foundation of the Jiangsu Higher Education Institutions of China (16KJB140022); National Natural Science Foundation of China (11834008, 11874217).

Disclosures

The authors declare no conflicts of interest.

References

1. L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nat. Photonics 3(9), 503–509 (2009). [CrossRef]

2. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). [CrossRef]

3. X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica, and L. V. Wang, “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol. 21(7), 803–806 (2003). [CrossRef]

4. D. M. Egolf, R. Chee, and R. J. Zemp, “Sparsity-based reconstruction for super-resolved limited-view photoacoustic computed tomography deep in a scattering medium,” Opt. Lett. 43(10), 2221–2224 (2018). [CrossRef]

5. H. Moradi, S. Tang, and S. E. Salcudean, “Deconvolution based photoacoustic reconstruction with sparsity regularization,” Opt. Express 25(3), 2771–2789 (2017). [CrossRef]

6. D. Wang, Y. Wang, Y. Zhou, J. F. Lovell, and J. Xia, “Coherent-weighted three-dimensional image reconstruction in linear-array-based photoacoustic tomography,” Biomed. Opt. Express 7(5), 1957–1965 (2016). [CrossRef]

7. J. Kang, E.-K. Kim, J. Kwak, Y. Yoo, T.-K. Song, and J. Chang, “Optimal laser wavelength for photoacoustic imaging of breast microcalcifications,” Appl. Phys. Lett. 99(15), 153702 (2011). [CrossRef]

8. L. Xi, S. R. Grobmyer, L. Wu, R. Chen, G. Zhou, L. G. Gutwein, J. Sun, W. Liao, Q. Zhou, H. Xie, and H. Jiang, “Evaluation of breast tumor margins in vivo with intraoperative photoacoustic imaging,” Opt. Express 20(8), 8726–8731 (2012). [CrossRef]

9. L. Xi and H. Jiang, “High resolution three-dimensional photoacoustic imaging of human finger joints in vivo,” Appl. Phys. Lett. 107(6), 063701 (2015). [CrossRef]

10. J. R. Rajian, M. L. Fabiilli, J. B. Fowlkes, P. L. Carson, and X. Wang, “Drug delivery monitoring by photoacoustic tomography with an ICG encapsulated double emulsion,” Opt. Express 19(15), 14335–14347 (2011). [CrossRef]

11. L. Li, L. Zhu, C. Ma, L. Lin, J. Yao, L. Wang, K. Maslov, R. Zhang, W. Chen, J. Shi, and L. V. Wang, “Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high patiotemporal resolution,” Nat. Biomed. Eng. 1(5), 0071 (2017). [CrossRef]

12. T. Guan, W. Shang, H. Li, X. Yang, C. Fang, J. Tian, and K. Wang, “From Detection to Resection: Photoacoustic Tomography and Surgery Guidance with Au@liposome-ICG Nanoparticles in Liver Cancer,” Bioconjugate Chem. 28(4), 1221–1228 (2017). [CrossRef]

13. S. Shang, Z. Chen, Y. Zhao, S. Yang, and D. Xing, “Simultaneous imaging of atherosclerotic plaque composition and structure with dual-mode photoacoustic and optical coherence tomography,” Opt. Express 25(2), 530–539 (2017). [CrossRef]

14. L. Wang, P. Lei, X. Wen, P. Zhang, and S. Yang, “Tapered fiber-based intravascular photoacoustic endoscopy for high-resolution and deep-penetration imaging of lipid-rich plaque,” Opt. Express 27(9), 12832–12840 (2019). [CrossRef]

15. P. K. Upputuri and M. Pramanik, “Dynamic in vivo imaging of small animal brain using pulsed laser diode based photoacoustic tomography system,” J. Biomed. Opt. 22(09), 1–4 (2017). [CrossRef]

16. D. Wu, J. Yang, G. Zhang, and H. Jiang, “Noninvasive in vivo monitoring of collagenase induced intracerebral hemorrhage by photoacoustic tomography,” Biomed. Opt. Express 8(4), 2276 (2017). [CrossRef]

17. P. Zhang, L. Li, L. Lin, P. Hu, J. Shi, Y. He, L. Zhu, Y. Zhou, and L. V. Wang, “High-resolution deep functional imaging of the whole mouse brain by photoacoustic computed tomography in vivo,” J. Biophotonics 11(1), e201700024 (2018). [CrossRef]

18. S. K. Kalva, P. K. Upputuri, and M. Pramanik, “High-speed, low-cost, pulsed-laser-diode-based second-generation desktop photoacoustic tomography system,” Opt. Lett. 44(1), 81–84 (2019). [CrossRef]

19. S. A. Goss, R. L. Johnston, and F. Dunn, “Comprehensive compilation of empirical ultrasonic properties of mammalian tissues,” J. Acoust. Soc. Am. 64(2), 423–457 (1978). [CrossRef]

20. X. Jin, C. Li, and L. V. Wang, “Effects of acoustic heterogeneities on transcranial brain imaging with microwave-induced thermoacoustic tomography,” Med. Phys. 35(7Part1), 3205–3214 (2008). [CrossRef]

21. H. Chao, L. Nie, R. W. Schoonover, L. V. Wang, and M. A. Anastasio, “Special Section on Photoacoustic Imaging and Sensing: Photoacoustic computed tomography correcting for heterogeneity and attenuation,” J. Biomed. Opt. 17(6), 061211 (2012). [CrossRef]

22. Y. Xu and L. V. Wang, “Effects of acoustic heterogeneity in breast thermoacoustic tomography,” IEEE Trans. Ultrason. Ferr. 50(9), 1134–1146 (2003). [CrossRef]

23. B. E. Treeby, T. K. Varslot, E. Z. Zhang, J. G. Laufer, and P. C. Beard, “Automatic sound speed selection in photoacoustic image reconstruction using an autofocus approach,” J. Biomed. Opt. 16(9), 090501 (2011). [CrossRef]

24. C. Yoon, J. Kang, S. Han, Y. Yoo, T. K. Song, and J. H. Chang, “Enhancement of photoacoustic image quality by sound speed correction: ex vivo evaluation,” Opt. Express 20(3), 3082–3090 (2012). [CrossRef]

25. C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, “Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media,” IEEE Trans. Med. Imaging 32(6), 1097–1110 (2013). [CrossRef]

26. J. Xia, C. Huang, K. Maslov, M. A. Anastasio, and L. V. Wang, “Enhancement of photoacoustic tomography by ultrasonic computed tomography based on optical excitation of elements of a full-ring transducer array,” Opt. Lett. 38(16), 3140–3143 (2013). [CrossRef]

27. M. A. Anastasio and T. P. Matthews, “Joint reconstruction of the initial pressure and speed of sound distributions from combined photoacoustic and ultrasound tomography measurements,” Inverse Probl. 33(12), 124002 (2017). [CrossRef]

28. J. Yin, C. Tao, P. Cai, and X. Liu, “Photoacoustic tomography based on the Green’s function retrieval with ultrasound interferometry for sample partially behind an acoustically scattering layer,” Appl. Phys. Lett. 106(23), 503 (2015). [CrossRef]

29. J. Jose, R. G. Willemink, W. Steenbergen, C. H. Slump, T. G. van Leeuwen, and S. Manohar, “Speed-of-sound compensated photoacoustic tomography for accurate imaging,” Med. Phys. 39(12), 7262–7271 (2012). [CrossRef]

30. S. Mandal, E. Nasonova, X. L. Deán-Ben, and D. Razansky, “Optimal self-calibration of tomographic reconstruction parameters in whole-body small animal optoacoustic imaging,” Photoacoustics 2(3), 128–136 (2014). [CrossRef]

31. C. Cai, X. Wang, K. Si, J. Qian, J. Luo, and C. Ma, “Feature coupling photoacoustic computed tomography for joint reconstruction of initial pressure and sound speed in vivo,” Biomed. Opt. Express 10(7), 3447–3462 (2019). [CrossRef]

32. W. Chen, C. Tao, and X. Liu, “Artifact-free imaging through a bone-like layer by using an ultrasonic-guided photoacoustic microscopy,” Opt. Lett. 44(5), 1273–1276 (2019). [CrossRef]

33. W. Rui, Z. Liu, C. Tao, and X. Liu, “Reconstruction of Photoacoustic Tomography Inside a Scattering Layer Using a Matrix Filtering Method,” Appl. Sci. 9(10), 2071 (2019). [CrossRef]

34. Z. Belhachmi, T. Glatz, and O. Scherzer, “A direct method for photoacoustic tomography with inhomogeneous sound speed,” Inverse Probl. 32(4), 045005 (2016). [CrossRef]

35. S. Feng, C. Kane, P. A. Lee, and A. D. Stone, “Correlations and Fluctuations of Coherent Wave Transmission through Disordered Media,” Phys. Rev. Lett. 61(7), 834–837 (1988). [CrossRef]

36. I. Freund, M. Rosenbluh, and S. Feng, “Memory Effects in Propagation of Optical Waves through Disordered Media,” Phys. Rev. Lett. 61(20), 2328–2331 (1988). [CrossRef]

37. Q. Ding, C. Tao, and X. Liu, “Photoacoustics and speed-of-sound dual mode imaging with a long depth-of-field by using annular ultrasound array,” Opt. Express 25(6), 6141–6150 (2017). [CrossRef]

Enhancement of photoacoustic tomography of acoustically inhomogeneous tissue by utilizing a memory effect

Abstract

Corrections

1. Introduction

2. Theory

3. Numerical simulations

4. Experiments

5. Summary

Funding

Disclosures

References

Cited By

Figures (8)

Equations (12)

Optics Express