Coherent-weighted three-dimensional image reconstruction in linear-array-based photoacoustic tomography.

While the majority of photoacoustic imaging systems used custom-made transducer arrays, commercially-available linear transducer arrays hold the benefits of affordable price, handheld convenience and wide clinical recognition. They are not widely used in photoacoustic imaging primarily because of the poor elevation resolution. Here, without modifying the imaging geometry and system, we propose addressing this limitation purely through image reconstruction. Our approach is based on the integration of two advanced image reconstruction techniques: focal-line-based three-dimensional image reconstruction and coherent weighting. We first numerically validated our approach through simulation and then experimentally tested it in phantom and in vivo. Both simulation and experimental results proved that the method can significantly improve the elevation resolution (up to 4 times in our experiment) and enhance object contrast.

1. Introduction

As an emerging imaging modality, photoacoustic (PA) tomography (PAT) has evolved rapidly over the past few years [1-3]. Different types of transducer arrays have been used in PAT, including ring-shaped arrays, liner arrays and spherical arrays. Among these, linear arrays possess several advantages, such as low price, handheld convenience and wide clinical recognition. However, the three-dimensional (3D) imaging capability of linear arrays is very limited due to the poor elevation resolution. Over the past few years, several methods have been proposed to address this limitation. These methods can be classified into two categories based on their principles. The first category improves the elevation resolution by converting the elevation direction into axial or lateral directions. The implementation normally involves complicated scanning geometries [4,5] or the use of additional linear arrays [6]. The second category, introduced recently, is based on the acoustic diffraction through a thin slit [7]. While this method does not require modifying the scanning geometry, the signal amplitude is degraded because most transmitted ultrasound waves will be blocked by the metal slit.

In this study, we proposed the first method to improve the elevation resolution for linear transducer array purely through image reconstruction. This method is based on a novel combination of coherent weighting (CW) and focal-line (FL) 3D image reconstruction. CW is an adaptive weighting technique that calculates the coherence of received photoacoustic signals and assigns a weighting factor to the delay-and-summed signals. The weighted image highlights coherently summed signals, which subsequently improves the spatial resolution [8,9]. While the CW concept was proposed almost two decades ago, it has never been applied to improve elevation resolution in a linear transducer array due to the lack of an efficient 3D reconstruction algorithm. The FL reconstruction algorithm, which was originally proposed for a ring-shaped transducer array with concave elements [10], addressed this issue, because the effect of an acoustic lens is equivalent to a curved concave element. The algorithm takes into account the detection sensitivity in three dimensions to improve elevation resolution up to the acoustic focal size (focal size of L7-4 is ~1.5 mm). Integrating CW and FL algorithms allows for significant improvement in elevation resolution through a single translational scan along the elevation direction. Table 1 shows a comparison of the methods mentioned in [4-7] and our CWFL method. In this comparison, we assume that all methods used a 10 Hz pulse repetition laser to image over 12 mm region in elevation. It can be seen that our method possesses the shortest acquisition time and the simplest scanning geometry. A common limitation for all methods mentioned in [4-7] is that their depth of field is limited by either the scanning geometry [4-6] or the slit [7]. Therefore, they only improve the elevation resolution in a confined region. In comparison, CWFL has a large depth of field and improves the elevation resolution throughout the axial axis. To better quantify the efficiency in resolution improvement, we defined a speed to resolution factor (SRF) calculated by dividing the imaging frame rate (frames/minute) by the elevation resolution. A higher SRF value indicates faster imaging speed and finer elevation resolution. It can be seen the CWFL method has the highest SRF among methods that do not require probe modification. While slit-PAT [7] exhibits a better SRF, the probe was modified by adding a thin metal slit, which degrades SNR (slit blocks more than half of the photoacoustic signal) and limits the ultrasound imaging capability. Considering all these factors, CWFL represents the most simple and efficient approach to improve the elevation resolution.

Table 1

Comparison of different methods

	Single side access system [4]	Whole body system [5]	Bi-direction scanning [6]	Slit PAT [7]	CWFL
Scanning method	Continuous and simultaneous rotation and translation scan	180-degree rotational scan + translational scan at each 1.5 degree interval	2 translational scan perpendicular to each other	1 translational scan	1 translational scan
Acquisition time (min)	5.12	80.00	2.67	0.20	0.20
Modification of transducer array	No	No	No	Yes (slit)	No
SNR	High	High	High	Low	High
Depth of field	Small	Small	Small	Medium	Large
SRF	1.14	0.0366	3.48	16.2	4.56

The feasibility of our method was first numerically tested using Field II [11,12], and then experimentally validated through three phantoms and two human wrist experiments. Both simulated and experimental results proved that the method can significantly improve the elevation resolution (up to 4 times in our experiment) and enhance the object contrast.

2. Method

2.1 CWFL algorithm

For a thin cylindrically focused transducer, there exists a focal line that goes across the element focus and is perpendicular to the x (axial)-z (elevation) plane (Fig. 1) of the element. Due to its special position, the wavefront of a photoacoustic wave emitted from a point on this line will reach the entire surface of the transducer element simultaneously, maximizing the receiving sensitivity. For an arbitrary point in 3D space, while the emitted wavefront has multiple paths to reach the transducer element, only the one that interacts with the FL results in the maximum received amplitude. Using this unique path to calculate the time of arrival in image reconstruction, we can improve the elevation resolution up to the height of the focal zone [10]. The principle of FL reconstruction can be formulated as where I is the total number of elements in the transducer array, J is the total frame number along elevation direction and c is the speed of sound. denotes the photoacoustic signal received by element i at time t and frame j.

Fig. 1

Schematic of the 2D reconstruction and 3D FL reconstruction concepts. A: point of reconstruction. A′: projection of A on the imaging plane. B: an artifact point induced by 2D reconstruction (AF = BF). O: center of the transducer element. F: intersection of focal line and OA′.

For a point A located at (Fig. 1), the acoustic time of arrival is calculated by To obtain and , we need to first project A to the 2D imaging plane to form A'. Connecting A' to the center of transducer element O forms a crossing point F with the focal line. then equals and equals . As discussed in Ref [10], this travel path has better receiving sensitivity than simply connecting points A and O in a conventional 3D reconstruction algorithm.

The coherent weighting factor (CWF) is calculated during the image reconstruction process. CWF is defined as Here all parameters hold the same definition as in Eq. (1). Based on this definition, the CWF value will increase during constructive summation when the projected signals are perfectly in phase, while it will decrease during destructive summation when the projected signals are random in phase. The value of CWF ranges from 0 to 1, depending on the degree of in-phase addition among the photoacoustic signals [13]. Applying CWF to the reconstructed image allows for enhancing the coherently summated signals while suppressing the out-of-phase signals and randomized noise.

2.2 Simulation

To test the CWFL algorithm, we first used Field II numerical simulation [11,12]. In the simulation, we generated three point sources, which were distributed along the axial direction at 35 mm, 40 mm, and 45 mm from the array, respectively. The parameters of the transducer were set to be the same as L7-4, which was a 128-element linear transducer array with 5 MHz center frequency, 25 mm focal distance, 7 mm element height and 0.25 mm element width [14]. For simplicity, in the simulation we used a curved element instead of an acoustic lens, to form the acoustic focus. The 3D data set was acquired by numerically scanning the array along the elevation direction.

2.3 Phantom and human wrist experiment

Experimental validations were first performed by imaging two phantoms: a three-tube phantom and a complex tube phantom. Both phantoms were made of 0.5-mm inner-diameter silicon tube filled with P-Pc (phosphorus phthalocyanine), an organic dye with a strong absorption at 1064 nm [15,16]. For the three-tube phantom, the tubes formed three straight lines along the lateral direction and they were placed at 37 mm, 46 mm and 55 mm away from the transducer, respectively. This phantom can be used to verify the resolution improvement at different axial positions. For the other phantom, the tube formed a more complex curved structure at 42 away from the transducer. The tube distances were chosen to better demonstrate improvement in elevation resolution (the elevation resolution of transducer array degrades with increasing distance from transducer focus). For both phantom experiments, 3D imaging was performed by moving the phantoms at a speed of 1 mm/s over 30 mm (step size 0.1 mm) elevation distance. A schematic of the experimental setup is shown in Fig. 2. Photoacoustic excitation was provided by a Nd:YAG laser (10 ns pulse duration) with 1064 nm output (near the peak absorption wavelength of the dye) and 10 Hz pulse repetition rate. The light was routed to the imaging region through a 1.25-cm-diameter fiber bundle. The light intensity was around 17.5 mJ/cm2 over the imaging plane. Photoacoustic signals were detected by a 128-element linear transducer array (5 MHz central frequency ATL/Philips L7-4). The detected signals were first amplified by 54 dB and then digitized by a 128-channel ultrasound data acquisition system (Vantage, Verasonics Inc.) with a sampling rate of 20 MS/s. After each laser pulse our imaging system simultaneously acquires photoacoustic data from all 128 transducer elements.

Fig. 2

Schematic of the PAT system and imaging geometry.

The human wrist experiment was performed in compliance with the University at Buffalo IRB protocol. The experimental setup differed from the phantom experiments in a way that the transducer array was moving (at a speed of 1.0 mm/s over 40 mm) while the objects (wrist) were stationary. The first subject held the wrist at 40 mm away from the transducer surface, while the second subject held the wrist at 47 mm away from the transducer surface. To ensure light coverage, we used a 5 cm-length line-fiber here and the fiber bundle moved together with the transducer array. The same 1064 nm laser was used in the experiment for excitation. The excited area at the skin surface was around 5 cm by 2 cm in size and the maximum intensity was 15 mJ/cm2, far below the maximum permissible exposure limits (100 mJ/cm2) defined by American National Standards Institute (ANSI).

3. Result and discussion

For comparison, the simulated point sources were reconstructed with three different methods: two-dimensional (2D) reconstruction (and stack of 2D images to form a 3D image), FL reconstruction and CWFL reconstruction. The results are shown in Fig. 3. The 2D reconstruction assumes that all photoacoustic signals come from the reconstruction (imaging) plane (Fig. 2), resulting in false projection of out-of-plane signals (Fig. 1, point B). This causes strong artifacts along the elevation direction (Fig. 3(A)) and subsequently degrades the elevation resolution to 2.6 mm (averaged resolution of the three points). Compared to 2D reconstruction, FL reconstruction calculates the appropriate time of flight of the signal in three dimensions. This procedure significantly suppresses the arc-shaped artifacts shown in Fig. 3(A) and improves the elevation resolution to 1.6 mm (Fig. 3(B)). For CWFL, we further improve the elevation resolution to 0.9 mm, which is almost two times improvement compared to the FL reconstruction result. The reconstructed points look very similar to ideal point sources (Fig. 3(C)).

Fig. 3

Reconstructed images of numerical simulation. (A) 2D reconstructed image. (B) FL reconstructed image. (C) CWFL reconstructed image.

For better comparison, the experimental data were also reconstructed by the three methods used in the simulation study. Figure 4(A) shows a photograph of the three-tube phantom. Figures 4(B)-4(D) show the depth encoded maximum amplitude projection (MAP) of the reconstructed images. For better illustration, the same threshold was applied to all images. The elevation resolution was quantified by averaging the full width at half-maximum (FWHM) of the tube along the elevation direction. For 2D reconstruction, because the elevation resolution degrades as the axial distance increases, the FWHM of the tube sequentially increases from 2.8 mm (tube 1) to 3.5 mm (tube 2) and to 4.0 mm (tube 3). This resolution change can be clearly seen in Fig. 4(B), where the tubes look broader from 1 to 3. The FL reconstruction method renders a more uniform elevation resolution at different axial positions. As shown in Fig. 4(C), the three tubes now exhibit similar width. The elevation resolution was quantified to be 1.8 mm for tube 1 and tube 2 and 1.7 mm for tube 3. Figure 4(D) is the CWFL reconstructed image, in which the three tubes show even narrower width than Fig. 4(C). The elevation resolution is 1.1 mm for tube 1 and tube 2 and 1.0 mm for tube 3. Compared to the 2D reconstruction, CWFL reconstruction improved the elevation resolution of tube 3 up to 4 times. This result successfully proves that our method performs well for objects placed at different distances from the transducer array. This imaging capability is essential for deep tissue PA imaging. Note that because of experimental noise and the small variation in the individual elements’ focusing, the experimental elevation resolutions are slightly worse than their simulation counterparts.

Fig. 4

Photograph and depth-encoded MAP of reconstructed images for the three-tube experiment (MAP was performed along the axial direction). (A) Photograph of the three-tube phantom. (B) PA image reconstructed with 2D reconstruction method. (C) PA image reconstructed with the FL reconstruction method. (D) PA image reconstructed with the CWFL reconstruction method.

The tube in the complex phantom was twisted to contain both horizontal and vertical structures (Fig. 5(A)). Figure 5(B) shows the depth-encoded image from 2D reconstruction. The background noise in the image is obvious, and the horizontal structures are blurry due to the poor elevation resolution. Figure 5(C) shows the depth-encoded image reconstructed with the FL method. Here the tube structures become thinner and the overall profile of the phantom can now be identified. However, the tube boundaries are still unclear and the background noise is still perceptible. Figure 5(D) shows the depth-encoded image reconstructed with the CWFL method. Here, the structure of the phantom is clearly shown with negligible noise. The elevation resolution for different methods was quantified along two dashed lines (Fig. 5(B)). The average elevation resolution of the five tube cross sections (labeled with 1 to 5 in Fig. 5(B)) was 3.5 mm for 2D reconstruction, 1.8 mm for FL reconstruction, and 1.1 mm for the CWFL reconstruction. We also quantified SNR by dividing the average signal within box 7 and 9 (Fig. 5(D)) with standard deviation of background within box 6 and box 8 (Fig. 5(D)). The SNR for the 2D reconstruction was only 18 (average of the sampling box set). FL reconstruction improved the SNR to 120, and the combined algorithm further improved the SNR to 245. Both the elevation resolution and the SNR for the complex tube experiment were similar to the single tube experiment. This experiment again successfully validated the superiority of the CWFL reconstruction method.

Fig. 5

Photograph and depth-encoded MAP of reconstructed images for the complex tube experiment (MAP was performed along the axial direction). (A) Photograph of the complex tube phantom. (B) PA image reconstructed with 2D reconstruction method. (C) PA image reconstructed with the FL reconstruction method. (D) PA image reconstructed with the CWFL reconstruction method.

To demonstrate the feasibility of CWFL algorithm for in vivo imaging, we imaged one right wrist and one left wrist of two healthy male volunteers, respectively, and reconstructed the data using all three methods (Fig. 6). Photographs of the wrists are shown in Figs. 6(A)-6(B), in which blood vessels cannot be seen clearly due to light scattering. MAP of reconstructed images are shown in Figs. 6(C)-6(H). By comparing the images reconstructed with different methods, we found that the CWFL algorithm provided the best elevation resolution and revealed the vessels most clearly, with the thinnest diameter. Along dashed-line marked blood vessel cross sections (Fig. 6(G) and Fig. 6(H)), we quantified the elevation resolution for the CWFL reconstruction to be 1.0 mm (subject 1) and 1.1 mm (subject 2). These results are better than 2D reconstruction (3.3 mm for subject 1 and 3.5 mm for subject 2) and FL reconstruction (1.8 mm for subject 1 and 2.0 mm for subject 2). Besides improvement in elevation resolution, our method significantly suppressed the background noise, and even revealed vessels (labeled with 1 and 2 in Fig. 6(G) and 3 and 4 in Fig. 6(H)) that were not shown in the 2D reconstructed images (Fig. 6(C) and Fig. 6(D)). These results successfully validated that our method works well for in vivo deep tissue PA imaging.

Fig. 6

MAP of reconstructed images from human-wrist experiments (MAP of C-H was performed along the axial direction over 12 mm range). Subject 1’s wrist was place at 40 mm away from the transducer surface. Subject 2’ wrist was placed at 47 mm away from the transducer surface. The vessels shown in C to H are 2~14 mm underneath the skin surface. (A-B) Photographs of the wrists of subject 1 and subject 2, respectively (the red dashed box indicates the imaging region). (C-D) Images reconstructed with the 2D method. (E-F) Images reconstructed with the FL method. (G-H) Images reconstructed with the CWFL method. The elevation resolution was calculated along the dashed lines.

4. Conclusion

In summary, we proposed a coherent-weighted 3D reconstruction technique for linear transducer array to improve the elevation resolution. By testing the feasibility of this method through simulation, phantom experiments and in vivo human-wrist experiments, we proved that the method can significantly improve the elevation resolution by up to 4 times in our experiment. In fact, because the intrinsic elevation resolution of a linear array degrades with increasing axial distance, the improvement could be even larger if the object is placed further away from the acoustic focus. Besides the improvement in elevation resolution, the method also suppresses noise and offers images with higher SNR. Compared to existing approaches to improve elevation resolution [4-7], our method features the simplest imaging geometry, retains both ultrasound and photoacoustic imaging capability, provides the largest depth of field, and possesses the fastest imaging speed. These advantages make our method very convenient for clinical imaging, as the subject only needs to stay stationary for several seconds. By using a higher speed laser (1064 nm lasers with up to 150 kHz pulse repetition rate are available from EdgeWave GmbH), our method can potentially achieve real-time 3D imaging at video rate, allowing photoacoustic imaging to be implemented in clinical 3D ultrasound probes.

One potential limitation of the technique is that the CWF may alter the quantitative information among pixels and render inaccurate results in functional imaging (e.g., quantification of oxygen saturation). However, this can be addressed by applying CWF after the functional analysis. Another interesting topic is to apply the CWF for more than one times on the FL reconstructed image to further improve elevation resolution. While this approach works in noise-free simulated result, we found that, in experiment, CWF may degrade the image features. Applying it for only once provides a good balance between image quality and elevation resolution. Nevertheless, our approach can be readily applied to current linear array PAT systems to provide 3D imaging with high elevation resolution, thereby, advancing the application of PAT.