Accounting for speed of sound variations in volumetric hand-held optoacoustic imaging

Abstract

Hand-held implementations of recently introduced real-time volumetric tomography approaches represent a promising path toward clinical translation of the optoacoustic technology. To this end, rapid acquisition of optoacoustic image data with spherical matrix arrays has attained exquisite visualizations of three-dimensional vascular morphology and function deep in human tissues. Nevertheless, significant reconstruction inaccuracies may arise from speed of sound (SoS) mismatches between the imaged tissue and the coupling medium used to propagate the generated optoacoustic responses toward the ultrasound sensing elements. Herein, we analyze the effects of SoS variations in three-dimensional hand-held tomographic acquisition geometries. An efficient graphics processing unit (GPU)-based reconstruction framework is further proposed to mitigate the SoS-related image quality degradation without compromising the high-frame-rate volumetric imaging performance of the method, essential for real-time visualization during hand-held scans.

Key words： speed of sound (SoS); graphics processing unit (GPU)

Introduction

The prospects of clinical translation of optoacoustic imaging have recently been bolstered by the introduction of hand-held systems employing concave tomographic geometries for optimal optoacoustic data collection [1–6]. This is in contrast to earlier approaches based on incorporating light delivery elements into standard planar or linear pulse-echo ultrasonography probes [7–9], commonly leading to sub-optimal optoacoustic imaging performance due to the limited-view data acquisition [10]. Yet, concave acquisition geometries do not permit close contact between the detection aperture and tissue surface. Thus, an acoustic coupling medium between the two surfaces is required [11,12], which may introduce an additional source of acoustic mismatches and alter the accuracy of tomographic reconstructions.

In general, the propagation of optoacoustically-excited waves is known to be affected by various acoustic phenomena, including scattering and reflections [13,14], acoustic attenuation [15–17], as well as speed of sound (SoS) variations in the medium [18–24]. Scattering and attenuation are generally negligible in soft tissues for the low MHz frequency range employed in most array-based optoacoustic scanners. On the other hand, the SoS in soft tissues can vary up to 10% with respect to that of water at room temperature [25], which may lead to significant reconstruction inaccuracies if the characteristics of the acoustic propagation path substantially differ for the different source-detector pairs. In small animal imaging systems with a broad or full tomographic coverage, the SoS distribution in the medium can be readily estimated tomographically by measuring the time-of-flight of ultrasound waves [26,27]. In this case, it has been shown that the quality of optoacoustic reconstructions can be improved by accounting for the heterogeneous SoS distribution instead of assuming a uniform SoS [19]. Proper selection of the SoS may alternatively be done with an autofocusing approach [28,29], which however may not be applicable to accurately reconstruct a relatively large region. However, the SoS distribution cannot be accurately measured in the case of a hand-held human imaging system having very partial tomographic coverage. In this work, we investigate on the effects of the SoS variations in volumetric hand-held optoacoustic imaging based on curved (concave) tomographic arrays and suggest a graphics processing unit (GPU)-based implementation of a back-projection algorithm that allows real-time operation for this configuration.

Materials and methods

Tomographic reconstruction in concave-array-based hand-held scanners

Fig.1 Propagation of an optoacoustic wave generated at an absorbing point within the imaged tissue to the detection point located in the coupling medium. Acoustic refraction takes place at the interface due to a difference in the speed of sound (SoS)

Full size|PPT slide

The acoustic propagation problem is schematically depicted in Fig. 1. As a first order approximation, let us assume that the SoS in tissue c_t is uniform and different from the SoS in the coupling medium (water) c_w. Furthermore, the interface between both media is considered to be flat, which is a realistic assumption when e.g. using a flat membrane pressed against the skin surface [11,30]. The distortion of the signals with respect to the case of a uniform acoustic propagation medium relates to a time-shift of the signals associated with the variations in their time-of-flight. The exact calculation of the time-of-flight involves calculating the propagation path by means of the Fermat’s principle or, equivalently, by considering the Snell’s law at the interface. Fermat’s principle states that the actual path of a wave propagating between two points is that for which the time of flight is minimum. If two different values of the SoS are considered, the propagation path between r′_j (the j-th voxel in the imaged region of interest) and r_i (the i-th detection element) is characterized by the intersection point at the interface between the two media (defined by the distance x in Fig. 1) According to Fermat’s principle, the value of x can be calculated by minimizing the time of flight between r′_j and r_i, i.e., t_ij satisfies

(1){\displaystyle \frac{\partial }{\partial x}}({{\displaystyle t}}_{ij})={\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{{{\displaystyle d}}_{\mathrm{t}}}{{{\displaystyle c}}_{\mathrm{t}}}}+{\displaystyle \frac{{{\displaystyle d}}_{\mathrm{w}}}{{{\displaystyle c}}_{\mathrm{w}}}}\right)=\mathrm{0},

where t_ij is the time-of-flight of ultrasound waves between r′_j and r_i along the path defined in Fig. 1, being d_t and d_w the propagating distances in tissue and water, respectively. The time-of-flight for the calculated values of x for all points of the region of interest (ROI) and the detection locations can be fed into a back-projection-based reconstruction method, which calculates the optical absorption at r′_j via [31,32]

(2)H({{\displaystyle r\prime }}_j)={{\displaystyle \Sigma }}_i\hspace{0.17em}{{\displaystyle p}}_{\mathrm{f}\mathrm{i}\mathrm{l}\mathrm{t}}({{\displaystyle r}}_i,{{\displaystyle t}}_{ij}).

The filtered pressure p_filt(r_i,t_ij) in Eq. (2) represents the measured pressure variations after implementing all the relevant corrections, such as the finite impulse response (FIR) filter combining deconvolution with the impulse response of the sensing elements, noise-reduction filtering and differentiation [31]. GPU-based implementation of the three-dimensional (3D) back-projection formula has recently enabled real-time volumetric visualization of structures at a frame rate determined by the pulse repetition frequency of the laser, i.e., tens of Hz [31]. However, direct solution of Eq. (1) to account for the SoS variations is associated with a fourth-order polynomial equation, whose analytical solution involves a large number of computations. Alternatively, Eq. (1) can be iteratively solved with Newton’s method via

(3){{\displaystyle x}}_{n+\mathrm{1}}={{\displaystyle x}}_n-{\displaystyle \frac{f(x)}{f^{\text{'}}(x)}},

where f(x)={\displaystyle \frac{\partial }{\partial x}}({{\displaystyle t}}_{ij}). The number of iterations in Eq. (3) determines the accuracy and the reconstruction time, although a few iterations are generally sufficient for convergence. For the first iteration, x0 was taken by considering a straight line between r_iand r′_j(Fig. 1). This first guess represents the straight acoustic ray model (SARM), which enables a good estimation of the time-of-flight for small SoS variations and a relatively small angle of incidence at the interface [19].

Numerical simulations

Numerical simulations were first done to assess the image distortions produced when assuming a uniform SoS. Without loss of generality, a two dimensional case was considered for the ease of computations. Five point absorbers were positioned at depths ranging from 0 to 2 cm within a tissue, as depicted in Fig. 2(a). The SoS within the coupling medium (water) was taken as 1500 m/s whereas a 10% increase (1650 m/s) was considered in the tissue as also indicated in Fig. 2(a). The optoacoustic signals at 91 equally-spaced locations along an arc covering an angle of 90° (Fig. 2(a)) were calculated with the k-wave toolbox [23] and band-pass filtered between 0.1 and 7.5 MHz. In fact, this sensor configuration is analogous to a recently-developed 3D optoacoustic hand-held probe [33].

Experimental measurements

A phantom experiment was subsequently done to validate the simulation results. For this, eight 200 mm diameter polyethylene microspheres were embedded in an agar matrix (1.3% agar powder by weight). The agar solution contained 33% of glycerine (by volume) to mimic an approximately 10% increase in the SoS with respect to water [19]. Optoacoustic imaging was done by illuminating the phantom with nanosecond-duration light pulses from an optical parametric oscillator (OPO)-based laser (Innolas, Krailling, Germany) set to 750 nm wavelength. The generated signals were collected by a matrix array detector (Imasonic SaS, Voray, France) that consisted of 256 individual piezocomposite elements distributed along a spherical surface covering an angle of 90° [33]. The phantom and the matrix array were immersed in a water tank for acoustic coupling. Prior to reconstruction, the acquired signals were deconvolved with the impulse response of the transducer elements and band-pass filtered between 0.1 and 7.5 MHz. Tomographic reconstruction was then performed using Eq. (2) for a volume of interests of 4×4×20 mm³ (80×80×400 voxels).

The performance of the proposed reconstruction algorithm was finally demonstrated by 3D hand-held imaging of a human palm. The laser wavelength was set to 800 nm in this case, corresponding approximately to the near-infrared isosbestic point of hemoglobin. Much like in the phantom experiment, the acquired signals were deconvolved with the impulse response of the transducer elements and band-pass filtered between 0.1 and 7.5 MHz prior to reconstruction. Optoacoustic image volumes of 12×12×12 mm³ containing 240×240×240 voxels were reconstructed for each laser shot without applying signal averaging.

Results

The simulation results are shown in Fig. 2. Specifically, the reconstructed images within the dashed area in Fig. 2(a) are shown in Figs. 2(b)−2(e). Figures 2(b) and 2(c) were obtained by assuming a uniform value for the SoS in the medium of 1500 and 1590 m/s, respectively. While the superficial absorber can be accurately reconstructed when assuming a uniform SoS equal to that in water (Fig. 2(b)), shapes of the deeper absorbers are distorted. It is indeed possible to heuristically choose a SoS value to optimize the reconstruction for a particular absorber at a given depth (Fig. 2(c)). However, an accurate reconstruction in the entire ROI cannot be achieved in this way. On the other hand, Figs. 2(d) and 2(e) show the reconstructions obtained by considering the actual SoS distribution for 0 and 10 iterations in Eq. (3). Here all the five absorbers are accurately reproduced even when considering a random initial guess x₀ (Fig. 1), confirming usefulness of the SARM approximation.

Maximum intensity projections (MIPs) of the 3D reconstructions corresponding to the phantom experiment are displayed in Fig. 3. In particular, Figs. 3(a)−3(c) display images reconstructed by assuming a homogeneous medium with SoS values of 1505, 1540 and 1575 m/s, respectively. Much like in the simulation, it is not possible to simultaneously resolve both shallow and deeply embedded absorbers. On the other hand, Fig. 3(d) showcases reconstructions rendered by considering two different SoS in the media (1505 and 1650 m/s for the water and phantom, respectively). In this case, all absorbers can be resolved with the fast and simple SARM approach, which was implemented on a GPU as described in Ref. [31]. The computational times achieved with this implementation are displayed in Fig. 3(e) as a function of the number of reconstructed voxels. The computation times for the homogeneous case are also displayed as a reference. Considering common laser pulse repetition rates of tens of Hz and typical reconstructed volumes containing about one million voxels, the SARM approximation can be confidently used for rendering 3D optoacoustic reconstructions in real time.

A comparison of the results obtained with homogeneous and heterogeneous SoS reconstruction approaches for the hand-held human scan is shown in Fig. 4. Specifically, Figs. 4(a)−4(c) show the MIPs of a reconstructed volume for depths ranging from 0 to 4.5 mm within the tissue and Figs. 4(d)−4(f) display the MIPs for a depth range between 4.5 and 9 mm. The images in the first and second columns of Fig. 4 were rendered by assuming a uniform SoS of 1510 and 1535 m/s, respectively. One notes that superficial structures (e.g. vessels V1 and V2 in Fig. 4(a)) are more accurately resembled for a SoS of 1510 m/s whereas SoS value of 1535 m/s better fits for the deeper areas (e.g. vessel V3 in Fig. 4(e)). On the other hand, better vascular contrast is attained for both superficial and deeper structures when the two SoS reconstruction is applied, namely, SoS value of 1510 m/s for the coupling medium versus 1575 m/s within the tissue, as displayed in the third column of Fig. 4. Since the images shown in Fig. 4 present localized features (blood vessels) on a large uniform background, their contrast was quantified according to Weber fraction (I−I_b)/I_b [34], where I was taken as the maximum of an image and the background I_b was taken as the average. The values obtained were 15.2, 14.0 and 16.0 for Figs. 4(a)−4(c) and 9.6, 9.2 and 9.6 for Figs. 4(d)−4(e), respectively, i.e., the highest values of the Weber fraction are rendered with the two SoS reconstruction procedure.

Fig.2 Simulated optoacoustic reconstructions for point absorbers embedded at different depths within the imaged tissue. (a) Location of the five absorbers within the tissue and the detection geometry (blue dots). The values of the speed of sound (SoS) in tissue c_t and water c_w are indicated. Scalebar – 10 mm. Reconstructed images of the absorbers when assuming a uniform SoS of 1500 and 1590 m/s in the entire medium are shown in (b) and (c), respectively. The region of interest is labeled by a dashed rectangle in (a). Reconstructions obtained by considering the actual SoS distribution with 0 and 10 iterations in Eq. (3) are shown in (d) and (e), respectively

Full size|PPT slide

Fig.3 Experimental imaging results for an agar-glycerine phantom with a higher SoS than water containing 200 mm absorbing microspheres. (a)−(c) Reconstructions obtained by considering a uniform SoS of 1505, 1540 and 1575 m/s, respectively. (d) Reconstructions obtained by considering SoS values of 1505 and 1650 m/s in water and the phantom material, respectively. Scalebar – 2 mm. (e) Comparison of the computational time when considering homogeneous (blue triangles) or heterogeneous (red squares) SoS distributions as a function of the number of reconstructed voxels

Full size|PPT slide

Fig.4 3D optoacoustic images acquired from a human palm. (a)−(c) Maximum intensity projections along the depth direction for imaging depths ranging from 0 to 4.5 mm. (d)−(f) Maximum intensity projections along the depth direction for imaging depths ranging from 4.5 to 9 mm. First column – reconstructions assuming a uniform SoS of 1510 m/s. Second column – reconstructions assuming a uniform SoS of 1535 m/s. Third column – reconstructed image assuming different SoS in water (1510 m/s) and tissue (1575 m/s). Scalebar – 4 mm

Full size|PPT slide

Conclusions

The results showcased indicate that the spatial resolution and image contrast of volumetric optoacoustic reconstructions can be optimally maintained across different tissue depths by considering different SoS in tissue and the coupling medium. The proposed methodology serves for image quality enhancement but may also contribute to increasing the effective imaging depth by better focusing on the deeper structures. More sophisticated reconstruction procedures considering a heterogeneous SoS distribution within the tissue may further help improving the images. However, the main advantage of the currently proposed dual SoS approach lies in its computational simplicity that does not compromise the real-time image rendering performance, which guarantees its applicability in clinical practice.

Acknowledgement

The authors acknowledge support from the European Research Council Consolidator grant ERC-2015CoG-682379.

References

1 Dean-Ben X L, Razansky D. Adding fifth dimension to optoacoustic imaging: volumetric time-resolved spectrally enriched tomography. Light, Science & Applications, 2014,3(1): e137

2 Buehler A, Kacprowicz M, Taruttis A, Ntziachristos V. Real-time handheld multispectral optoacoustic imaging. Optics Letters, 2013, 38(9): 1404–1406 

PMID

3 Taruttis A, Ntziachristos V. Advances in real-time multispectral optoacoustic imaging and its applications. Nature Photonics, 2015, 9(4): 219–227 

4 Wang L V, Yao J. A practical guide to photoacoustic tomography in the life sciences. Nature Methods, 2016, 13(8): 627–638 

PMID

5 Stoffels I, Morscher S, Helfrich I, Hillen U, Leyh J, Burton N C, Sardella T C, Claussen J, Poeppel T D, Bachmann H S, Roesch A, Griewank K, Schadendorf D, Gunzer M, Klode J. Metastatic status of sentinel lymph nodes in melanoma determined noninvasively with multispectral optoacoustic imaging. Science Translational Medicine, 2015, 7(317): 317ra199 160;

PMID

6 Deán-Ben X L, Merčep E, Razansky D. Hybrid-array-based optoacoustic and ultrasound (OPUS) imaging of biological tissues. Applied Physics Letters, 2017, 110(20): 203703 

7 Kim C, Erpelding T N, Jankovic L, Pashley M D, Wang L V. Deeply penetrating in vivo photoacoustic imaging using a clinical ultrasound array system. Biomedical Optics Express, 2010, 1(1): 278–284 160;

PMID

8 Fronheiser M P, Ermilov S A, Brecht H P, Conjusteau A, Su R, Mehta K, Oraevsky A A. Real-time optoacoustic monitoring and three-dimensional mapping of a human arm vasculature. Journal of Biomedical Optics, 2010, 15(2): 021305 160;

PMID

9 Piras D, Grijsen C, Schütte P, Steenbergen W, Manohar S. Photoacoustic needle: minimally invasive guidance to biopsy. Journal of Biomedical Optics, 2013, 18(7): 070502 

PMID

10 Deán-Ben X L, Razansky D. On the link between the speckle free nature of optoacoustics and visibility of structures in limited-view tomography. Photoacoustics, 2016, 4(4): 133–140 160;

PMID

11 Neuschmelting V, Burton N C, Lockau H, Urich A, Harmsen S, Ntziachristos V, Kircher M F. Performance of a multispectral optoacoustic tomography (MSOT) system equipped with 2D vs. 3D handheld probes for potential clinical translation. Photoacoustics, 2016, 4(1): 1–10 160;

PMID

12 Mercep E, Dean Ben X L, Razansky D. Combined pulse-echo ultrasound and multispectral optoacoustic tomography with a multi-segment detector array. IEEE Transactions on Medical Imaging, 2017, doi: 10.1109/TMI.2017.2706200

PMID

13 Deán-Ben X L, Ntziachristos V, Razansky D. Artefact reduction in optoacoustic tomographic imaging by estimating the distribution of acoustic scatterers. Journal of Biomedical Optics, 2012, 17(11): 110504 160;

PMID

14 Deán-Ben X L, Ma R, Rosenthal A, Ntziachristos V, Razansky D. Weighted model-based optoacoustic reconstruction in acoustic scattering media. Physics in Medicine and Biology, 2013, 58(16): 5555–5566 160;

PMID

15 Treeby B E. Acoustic attenuation compensation in photoacoustic tomography using time-variant filtering. Journal of Biomedical Optics, 2013, 18(3): 036008 160;

PMID

16 Huang C, Nie L, Schoonover R W, Wang L V, Anastasio M A. Photoacoustic computed tomography correcting for heterogeneity and attenuation. Journal of Biomedical Optics, 2012, 17(6): 061211 160;

PMID

17 Deán-Ben X L, Razansky D, Ntziachristos V. The effects of acoustic attenuation in optoacoustic signals. Physics in Medicine and Biology, 2011, 56(18): 6129–6148 160;

PMID

18 Modgil D, Anastasio M A, La Rivière P J. Image reconstruction in photoacoustic tomography with variable speed of sound using a higher-order geometrical acoustics approximation. Journal of Biomedical Optics, 2010, 15(2): 021308 160;

PMID

19 Deán-Ben X L, Ntziachristos V, Razansky D. Effects of small variations of speed of sound in optoacoustic tomographic imaging. Medical Physics, 2014, 41(7): 073301 160;

PMID

20 Wurzinger G, Nuster R, Paltauf G. Combined photoacoustic, pulse-echo laser ultrasound, and speed-of-sound imaging using integrating optical detection. Journal of Biomedical Optics, 2016, 21(8): 086010 160;

PMID

21 Haltmeier M, Nguyen L V. Analysis of iterative methods in photoacoustic tomography with variable sound speed. SIAM Journal on Imaging Sciences, 2017, 10(2): 751–781 160;

22 Li L, Zhu L R, Ma C, Lin L, Yao J J, Wang L D, Maslov K, Zhang R Y, Chen W Y, Shi J H, Wang L V. Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution. Nature Biomedical Engineering, 2017, 1(5): Art. No. 0071

23 Treeby B E, Cox B T. k-wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields. Journal of Biomedical Optics, 2010, 15(2): 021314 160;

PMID

24 Huang C, Wang K, Nie L, Wang L V, Anastasio M A. Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media. IEEE Transactions on Medical Imaging, 2013, 32(6): 1097–1110160;

PMID

25 Szabo T L. Diagnostic ultrasound imaging inside out. Burlington, MA: Elsevier Academic Press, 2004, p. xxii, 549 p

26 Jose J, Willemink R G, Resink S, Piras D, van Hespen J C, Slump C H, Steenbergen W, van Leeuwen T G, Manohar S. Passive element enriched photoacoustic computed tomography (PER PACT) for simultaneous imaging of acoustic propagation properties and light absorption. Optics Express, 2011, 19(3): 2093–2104 160;

PMID

27 Xia J, Huang C, Maslov K, Anastasio M A, Wang L V. Enhancement of photoacoustic tomography by ultrasonic computed tomography based on optical excitation of elements of a full-ring transducer array. Optics Letters, 2013, 38(16): 3140–3143 160;

PMID

28 Treeby B E, Varslot T K, Zhang E Z, Laufer J G, Beard P C. Automatic sound speed selection in photoacoustic image reconstruction using an autofocus approach. Journal of Biomedical Optics, 2011, 16(9): 090501  

PMID

29 Mandal S, Nasonova E, Deán-Ben X L, Razansky D. Optimal self-calibration of tomographic reconstruction parameters in whole-body small animal optoacoustic imaging. Photoacoustics, 2014, 2(3): 128–136 160;

PMID

30 Deán-Ben  X L, Fehm T F, Gostic M, Razansky D. Volumetric hand-held optoacoustic angiography as a tool for real-time screening of dense breast. Journal of Biophotonics, 2016, 9(3): 253–259  

PMID

31 Deán-Ben X L, Ozbek A, Razansky D. Volumetric real-time tracking of peripheral human vasculature with GPU-accelerated three-dimensional optoacoustic tomography. IEEE Transactions on Medical Imaging, 2013, 32(11): 2050–2055 160;

PMID

32 Xu M, Wang L V. Universal back-projection algorithm for photoacoustic computed tomography. Physical Review E:  Statistical, Nonlinear, and Soft Matter Physics, 2005, 71(1 Pt 2): 016706

33 Deán-Ben X L, Razansky D. Portable spherical array probe for volumetric real-time optoacoustic imaging at centimeter-scale depths. Optics Express, 2013, 21(23): 28062–28071 160;

PMID

34 Peli E. Contrast in complex images. Journal of the Optical Society of America A, Optics and Image Science, 1990, 7(10): 2032–2040 160;

PMID