Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Sub-Nyquist sampling-based high-frequency photoacoustic computed tomography

Open Access Open Access

Abstract

High-frequency (greater than 30 MHz) photoacoustic computed tomography (PACT) provides the opportunity to reveal finer details of biological tissues with high spatial resolution. To record photoacoustic signals above 30 MHz, sampling rates higher than 60 MHz are required according to the Nyquist sampling criterion. However, the highest sampling rates supported by existing PACT systems are typically within the range of 40–60 MHz. Herein, we propose a novel PACT imaging method based on sub-Nyquist sampling. The results of numerical simulation, phantom experiment, and in vivo experiment demonstrate that the proposed imaging method can achieve high-frequency PACT imaging with a relatively low sampling rate. An axial resolution of 22 μm is achieved with a 30-MHz transducer and a 41.67-MHz sampling rate. To the best of our knowledge, this is the highest axial resolution ever achieved in PACT based on a sampling rate of not greater than 60 MHz. This work is expected to provide a practical way for high-frequency PACT imaging with limited sampling rates.

© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

Photoacoustic (PA) computed tomography (PACT) is a fast-evolving biomedical imaging modality in the last two decades [1]. Based on the PA effect, PACT can achieve high spatiotemporal resolution at large imaging depths and thus has vast applications in a range of biomedical fields [2], such as neurology [3], oncology [4], rheumatology [5], surgical navigation [6], and other areas of medical imaging problems. High-resolution PACT provides the opportunity to reveal finer details of biological tissues [79]. The spatial resolution of a PACT system can be affected by multiple system factors and is closely related to the characteristics of ultrasonic detectors, including the center frequency, bandwidth, aperture size, and view angle [10]. To enhance the spatial resolution and thus improve the image quality of PACT, one can simply increase the center frequency of ultrasonic detectors [9]. According to the Nyquist sampling criterion, high-frequency PACT imaging demands data acquisition (DAQ) systems with analog-to-digital converters (ADC) of high sampling rates. In practice, the highest sampling rates supported by existing PACT systems are typically within the range of 40–60 MHz [11], which leads to a Nyquist frequency of not greater than 30 MHz.

To achieve high-frequency PACT imaging with bandwidths covering frequencies above 30 MHz, researchers usually employ DAQ systems with sampling rates higher than 60 MHz [8,1214]. However, such DAQ systems are usually not cost-effective. Inspired by a delayed-excitation DAQ method previously reported for high-frequency ultrasound imaging [15], Fu and Jokerst proposed a new method called interleave-sampled PA imaging that enables high-frequency PACT imaging with a relatively low sampling rate, e.g., a 41.67-MHz sampling rate with a 30-MHz transducer [16]. The proposed method requires two acquisitions that are precisely shifted with each other by half of the sampling period and are sampled at the same sampling rate (e.g., 41.67 MHz). Interleaving these two acquisitions forms a virtual acquisition that is equivalent to the acquisition sampled at a doubled sampling rate (e.g., 83.33 MHz). The authors achieved a 63-μm axial resolution and a 91-μm lateral resolution via the proposed imaging method.

In this Letter, we propose another novel imaging method based on sub-Nyquist sampling that can also achieve high-frequency PACT imaging with a relatively low sampling rate. Sub-Nyquist sampling, also known as bandpass sampling, allows for sampling bandpass signals with a sampling rate less than twice the highest frequency content [17] and has been widely used in ultrasound imaging for reducing computational complexity [18]. The idea of the proposed method is based on the fact that the frequency response of high-frequency ultrasonic detectors is usually band-limited.

Suppose we have a bandpass signal with a center frequency of fc and a bandwidth of B. According to the theory of sub-Nyquist sampling [17], we can sample the bandpass signal without aliasing at a rate of fs that is constrained to

$$\frac{{2{f_\textrm{c}} - B}}{m} \ge {f_\textrm{s}} \ge \frac{{2{f_\textrm{c}} + B}}{{m + 1}},$$
where m is an arbitrary, positive integer ensuring that fs  ≥ 2B. We can rearrange Eq. (1) as
$$B \le \frac{{2{f_\textrm{c}}}}{{2m + 1}} \le \frac{{2{f_\textrm{c}}}}{3}.$$
Smaller bandwidth brings larger m and thus smaller fs. Since bandwidth is essential for high-quality PACT imaging [9], we should choose smaller m, e.g., m = 1 in this Letter. In this case, the highest bandwidth supported by sub-Nyquist sampling is 2fc/3, and the corresponding fs is 4fc/3. Figure 1(a) illustrates the basic principle of the proposed sub-Nyquist sampling-based high-frequency PACT imaging method (m = 1, B = 2fc/3, fs = 4fc/3). Firstly, we employ a bandpass filter (passband: [2/3fc, 4/3fc]) to eliminate the noise and undesired spectra outside the band of interest [Figs. 1(b) and 1(c)]. Secondly, we sample the filtered signal at 4fc/3. The desired spectrum in the high-frequency band ([2fc/3, 4fc/3]) is replicated in the low-frequency band ([0, 2fc/3]) with spectral inversion after ADC digitization [Fig. 1(d)]. Thirdly, we fold the replicated spectrum back into the original high-frequency band ([2fc/3, 4fc/3]) along 2fc/3 [Fig. 1(e)]. The sampling rate is doubled in this step. The recovered spectrum shown in Fig. 1(e) is identical to the original spectrum of interest shown in Fig. 1(c). Finally, the recovered signal is used for PACT image reconstruction. Specifically, the third step (recovery of aliased signal) can be implemented as follows. Suppose F(ω) is the Fourier spectrum of the aliased digital signal, f(t). To fold F(ω) back into the original high-frequency band, we construct another spectrum, G(ω), as
$$G(\omega ) = \left\{ {\begin{array}{@{}ll@{}}F({{4{f_{\rm c}}} / 3} - |\omega |),&{{{\textrm{if}}\;2{f_{\rm c}}} / 3}< |\omega |< {{4{f_{\rm c}}} / 3}\\ 0, & {\textrm{otherwise}}.\end{array}} \right.$$
Taking the inverse Fourier transform of G(ω) gives the recovered digital signal g(t), the sampling rate of which is twice that of f(t).

 figure: Fig. 1.

Fig. 1. Schematic illustrating the basic principle of the proposed sub-Nyquist sampling-based high-frequency PACT imaging method. (a) Flow chart of the proposed method. (b) Original analog signal. (c) Filtered analog signal. (d) Replicated digital signal with spectral inversion. (e) Recovered digital signal.

Download Full Size | PDF

It is critical to ensure that the original analog signals are filtered with an analog bandpass filter before being sampled by the ADC to eliminate the noise and undesired spectra outside the band of interest. Otherwise, these noise and undesired spectra would be added to the replicated spectrum after digitization. Under such circumstances, the original spectrum of interest would be contaminated and could not be recovered accurately.

To demonstrate the performance of the proposed imaging method, numerical simulation, phantom experiment, and in vivo experiment were conducted. The numerical simulation was performed by imaging a microsphere. Considering that smaller microspheres yield higher frequency content, the diameter of the microsphere was set to 20 μm to represent a point target while matching the detection bandwidth of the transducer to maximize detection sensitivity as much as possible. An analytical model [19] was employed for generating PA signals of the microsphere, which were then received by a 30-MHz linear transducer array with a bandwidth of 20 MHz. The microsphere was placed 8 mm from the transducer array, as shown in Fig. 2(a). The received PA signals were first digitized at 1 GHz to serve as the original analog signals. To compare the performance of the proposed method with other methods, the digitized signals were separately filtered with multiple filtering configurations and down-sampled to different sampling rates (Table 1). Herein, we suppose the highest sampling rate supported by the DAQ is 62.5 MHz. Sampling rates of 83.33 and 125 MHz were used for simulating interleaved sampling based on 41.67 and 62.5 MHz, respectively. A highpass filter and a lowpass filter were combined to act as a bandpass filter for Config. 1. For Configs. 2 to 5, a lowpass filter configured with different cutoff frequencies was employed to avoid aliasing resulting from frequency content above half of the sampling rates. Varying sampling rates for Configs. 2 to 5 leads to varying cutoff frequencies of the lowpass filter. The cutoff frequencies and spectral responses of the above filters can be found in Table 1 and Fig. S1 (Supplement 1). PACT image reconstruction was performed based on the filtered back-projection algorithm [19].

 figure: Fig. 2.

Fig. 2. Numerical simulation confirming the feasibility of the proposed imaging method. (a) Schematic of the simulation setup showing a 20-μm microsphere placed 8 mm from a linear transducer array. (b) Ground truth for the simulation obtained with Config. 6. (c)–(g) Imaging results obtained with Configs. 1 to 5. (h) Summation of (c) and (d). (i) and (j) Axial and lateral resolution comparisons.

Download Full Size | PDF

Tables Icon

Table 1. Sampling Configurations for Performance Comparison

Figure 2(b) shows the imaging result obtained based on Config. 6. Due to the limited detector bandwidth [see Fig. S2(f), Supplement 1] and view angle [9], the reconstructed microsphere shown in Fig. 2(b) is not round-shaped. Nevertheless, Fig. 2(b) is regarded as the ground truth in this simulation. Figures 2(c)–2(g) are the imaging results obtained with Configs. 1 to 5. Negativity artifacts shown in Figs. 2(b)–2(g) result from the limited detector bandwidth and view angle [20]. Intensity profiles along the red dashed lines shown in Fig. 2(b) are plotted in Figs. 2(i) and 2(j) with the FWHM (full width at half maximum) quantified as axial (z axis) and lateral (x axis) resolution. Due to spectral aliasing and the loss of low-frequency content [see Fig. S2(a), Supplement 1] [9], Fig. 2(c) shows clearer ringing artifacts (white arrow) when compared with the ground truth. The hollow structure of the microsphere shown in Fig. 2(d), the moderately broadened microsphere shown in Fig. 2(e), and the mildly broadened microsphere shown in Fig. 2(f) both result from spectral aliasing and the loss of high-frequency content [see Figs. S2(b)–S2(d), Supplement 1]. Considering that Configs. 5 and 6 have similar patterns of frequency content [see Figs. S2(e) and S2(f), Supplement 1], Fig. 2(g) shows similar imaging result to the ground truth. Figure 2(c) is sharper than Fig. 2(e) owing to the sampling of a higher frequency content. The result obtained based on the proposed imaging method [Fig. 2(c)] is as sharp as those obtained based on the interleaved imaging method [Figs. 2(f) and 2(g)]. However, Fig. 2(c) suffers from severe ringing artifacts that can be well suppressed by adding Fig. 2(d) to Fig. 2(c), which yields Fig. 2(h). The above findings are well supported by the axial and lateral resolution comparisons shown in Figs. 2(i) and 2(j).

A fine hair (diameter: 10–20 μm) was embedded in an agar phantom and imaged to evaluate the performance of the proposed method in practice. The specifications of the imaging system can be found in [6] and will not be detailed here. The highest sampling rate supported by the imaging system is 62.5 MHz. A commercial 256-channel linear transducer array (MS400; Visualsonics, Inc., Canada) with a center frequency of 30 MHz and a bandwidth of 20 MHz was used for signal receiving. The imaging settings are the same as those described in the numerical simulation. Figures 3(a)–3(e) are the imaging results obtained with Configs. 1 to 5. As expected, the proposed imaging method [Fig. 3(a)] provides finer structure but worse ringing artifacts (white arrows) compared with the conventional imaging method [Fig. 3(c)]. The ringing artifacts shown in Fig. 3(a) are well suppressed by adding Fig. 3(b) to Fig. 3(a), as shown in Fig. 3(f). The interleaved imaging method also provides shaper results [Figs. 3(d) and 3(e)] compared with the conventional imaging method [Fig. 3(c)]. Figure 3(a) shows a smaller FWHM [Fig. 3(g)] and severe ringing artifacts, both resulting from spectral aliasing and the loss of low-frequency content, compared with Fig. 3(d). On the contrary, Fig. 3(b) shows a larger FWHM [Fig. 3(g)], resulting from spectral aliasing and the loss of high-frequency content, compared with Fig. 3(d). Theoretically, the combination of Configs. 1 and 2 should cancel each other out and produce a similar FWHM [Fig. 3(f)] as that obtained with Config. 4 [Fig. 3(d)]. However, due to spectral aliasing, Fig. 3(f) shows a slightly smaller FWHM. The above findings are well supported by the axial intensity profiles [Fig. 3(g)] and the digital spectra (see Fig. S3, Supplement 1).

 figure: Fig. 3.

Fig. 3. Phantom experiment of a fine hair evaluating the performance of the proposed method in practice. (a)–(e) Imaging results obtained with Configs. 1 to 5. (f) Summation of (a) and (b). (g) Axial intensity profiles along the red dashed line shown in (e).

Download Full Size | PDF

Finally, we conducted an in vivo experiment on a mouse ear to explore the potential of the proposed method in imaging complex biological tissues. Approval of all ethical and experimental procedures and protocols was granted by the Institutional Animal Care and Use Committee (IACUC) of the University of Science and Technology of China under No. USTCACUC1803065. Blood vessels of the mouse ear shown in Fig. 4(a) were imaged based on Configs. 1 to 5, with results presented in Figs. 4(b)–4(f). As expected, the proposed imaging method [Fig. 4(b)] provides higher spatial resolution but suffers from severe ringing artifacts when compared with the conventional imaging method [Fig. 4(d)]. The ringing artifacts shown in Fig. 4(b) are well suppressed by adding Fig. 4(c) to Fig. 4(b), as shown in Fig. 4(g). The performance of the combined imaging method [Fig. 4(g)]is comparable to that of the interleaved imaging method [Figs. 4(e) and 4(f)]. The findings are well supported by the characterized FWHM values [Figs. 4(h) and 4(i)] of the blood vessel shown in insets and the digital spectra (see Fig. S4, Supplement 1).

 figure: Fig. 4.

Fig. 4. In vivo experiment of a mouse ear exploring the potential of the proposed method in imaging complex biological tissues. (a) Photograph of the mouse ear showing lots of microvessels. (b)–(f) Imaging results of blood vessels (white arrows) obtained with Configs. 1 to 5. (g) Summation of (b) and (c). Insets are the zoomed-in images of the vessel in the yellow dashed boxes. (h) and (i) Intensity profiles of the blood vessel shown in insets along z axis (h) and x axis (i) with FWHM values characterized.

Download Full Size | PDF

Due to hardware limitation, the cutoff frequency of the highpass filter in Config. 1 was configured as 20 MHz, not 20.83 MHz. For consistency, the 20-MHz cutoff frequency was used in both the simulation and experiments. A specially designed numerical simulation shows that the impact of this tiny deviation (20.83 MHz versus 20 MHz) on imaging results is insignificant (see Fig. S5, Supplement 1). In practice, the analog filters are usually nonideal, having limited suppression ratios beyond the band of interest, which will lead to spectral aliasing and affect the imaging quality to some extent. Another specially designed numerical simulation indicates that high-performance analog filters are expected to reduce imaging artifacts resulting from the limited suppression ratios of the nonideal filters (see Fig. S6, Supplement 1).

The proposed method has the same level of complexity as the interleaved sampling method when two acquisitions are necessary (Configs. 1 and 2). However, under some circumstances, the proposed method may be superior to the interleaved sampling approach. Specifically, the former only needs one acquisition (Config. 1) to capture a high-frequency content, while the latter principally needs two. Although the results shown here suggest that two acquisitions are required to suppress ringing artifacts, we believe there may be scenarios (e.g., machine learning is used to suppress the ringing artifacts resulting from the loss of low-frequency components) where only one acquisition is sufficient in the future.

Herein, we demonstrate the concept based on a 30-MHz transducer. But actually, the proposed method can be used for signals with cutoff frequencies beyond 40 MHz when appropriate bandpass filters are available. High-frequency PACT primarily targets small structures (e.g., microvessels) whose frequency content matches the detection bandwidth of the transducer. Most low-frequency content produced by large structures will be lost due to the limited detection bandwidth of the transducer, but not the inappropriate setting of the sampling configuration. The loss of low-frequency content is also inevitable for the interleaved sampling method and only impacts the image quality of large structures, not small structures.

In summary, we proposed a novel imaging method based on sub-Nyquist sampling to achieve high-frequency PACT imaging with a relatively low sampling rate. The results show that the proposed imaging method provides better spatial resolution when compared to the conventional imaging method. It is necessary to add another acquisition that samples a low-frequency content to the proposed acquisition to suppress ringing artifacts induced by the proposed imaging method. The performance of the combined imaging method is comparable to that of the interleaved imaging method. We achieve an axial resolution of 22 μm based on a relative low sampling rate (i.e., 41.67 MHz) with the proposed imaging method. To the best of our knowledge, this is the highest axial resolution ever achieved in PACT based on a sampling rate of not greater than 60 MHz. This work is expected to provide a practical way for high-frequency PACT imaging with limited sampling rates.

Funding

National Key Research and Development Program of China (2022YFA1404400); National Natural Science Foundation of China (12174368, 61705216, 62122072); Anhui Provincial Department of Science and Technology (18030801138, 202203a07020020); University of Science and Technology of China (YD2090002015); Institute of Artificial Intelligence at Hefei Comprehensive National Science Center (23YGXT005).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

REFERENCES

1. L. V. Wang and S. Hu, Science 335, 1458 (2012). [CrossRef]  

2. L. V. Wang and J. Yao, Nat. Methods 13, 627 (2016). [CrossRef]  

3. S. Gottschalk, O. Degtyaruk, B. Mc Larney, et al., Nat. Biomed. Eng. 3, 392 (2019). [CrossRef]  

4. L. Lin, P. Hu, J. Shi, et al., Nat. Commun. 9, 2352 (2018). [CrossRef]  

5. J. Jo, C. Tian, G. Xu, et al., Photoacoustics 12, 82 (2018). [CrossRef]  

6. S. Liu, H. Wang, C. Zhang, et al., IEEE Trans. Biomed. Eng. 67, 2033 (2020). [CrossRef]  

7. G. Ku, X. Wang, G. Stoica, et al., Phys. Med. Biol. 49, 1329 (2004). [CrossRef]  

8. H. Zafar, A. Breathnach, H. M. Subhash, et al., J. Biomed. Opt. 20, 051021 (2015). [CrossRef]  

9. C. Tian, M. Pei, K. Shen, et al., Phys. Rev. Appl. 13, 014001 (2020). [CrossRef]  

10. C. Tian, C. Zhang, H. Zhang, et al., Rep. Prog. Phys. 84, 036701 (2021). [CrossRef]  

11. W. Choi, E.-Y. Park, S. Jeon, et al., Biomed. Eng. Lett. 8, 139 (2018). [CrossRef]  

12. J. Pan, Q. Li, Y. Feng, et al., Nat. Commun. 14, 3250 (2023). [CrossRef]  

13. L. Vionnet, J. Gateau, M. Schwarz, et al., IEEE Trans. Med. Imaging 33, 535 (2014). [CrossRef]  

14. G. Li, L. Li, L. Zhu, et al., J. Biomed. Opt. 20, 066010 (2015). [CrossRef]  

15. W. Qiu, J. Xia, Y. Shi, et al., IEEE Trans. Biomed. Eng. 65, 15 (2018). [CrossRef]  

16. L. Fu and J. Jokerst, Opt. Lett. 47, 3503 (2022). [CrossRef]  

17. R. G. Lyons, in Understanding Digital Signal Processing, 1st ed. (Addison-Wesley, 1996), pp. 23–47.

18. J. Kang, H. Yoon, C. Yoon, et al., IEEE Trans. Ultrason., Ferroelectr., Freq. Contr. 69, 2001 (2022). [CrossRef]  

19. M. Xu and L. V. Wang, Phys. Rev. E 71, 016706 (2005). [CrossRef]  

20. K. Shen, S. Liu, T. Feng, et al., J. Phys. D: Appl. Phys. 54, 074001 (2021). [CrossRef]  

Supplementary Material (1)

NameDescription
Supplement 1       Supplement 1

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)

Fig. 1.
Fig. 1. Schematic illustrating the basic principle of the proposed sub-Nyquist sampling-based high-frequency PACT imaging method. (a) Flow chart of the proposed method. (b) Original analog signal. (c) Filtered analog signal. (d) Replicated digital signal with spectral inversion. (e) Recovered digital signal.
Fig. 2.
Fig. 2. Numerical simulation confirming the feasibility of the proposed imaging method. (a) Schematic of the simulation setup showing a 20-μm microsphere placed 8 mm from a linear transducer array. (b) Ground truth for the simulation obtained with Config. 6. (c)–(g) Imaging results obtained with Configs. 1 to 5. (h) Summation of (c) and (d). (i) and (j) Axial and lateral resolution comparisons.
Fig. 3.
Fig. 3. Phantom experiment of a fine hair evaluating the performance of the proposed method in practice. (a)–(e) Imaging results obtained with Configs. 1 to 5. (f) Summation of (a) and (b). (g) Axial intensity profiles along the red dashed line shown in (e).
Fig. 4.
Fig. 4. In vivo experiment of a mouse ear exploring the potential of the proposed method in imaging complex biological tissues. (a) Photograph of the mouse ear showing lots of microvessels. (b)–(f) Imaging results of blood vessels (white arrows) obtained with Configs. 1 to 5. (g) Summation of (b) and (c). Insets are the zoomed-in images of the vessel in the yellow dashed boxes. (h) and (i) Intensity profiles of the blood vessel shown in insets along z axis (h) and x axis (i) with FWHM values characterized.

Tables (1)

Tables Icon

Table 1. Sampling Configurations for Performance Comparison

Equations (3)

Equations on this page are rendered with MathJax. Learn more.

2 f c B m f s 2 f c + B m + 1 ,
B 2 f c 2 m + 1 2 f c 3 .
G ( ω ) = { F ( 4 f c / 3 | ω | ) , if 2 f c / 3 < | ω | < 4 f c / 3 0 , otherwise .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.