2.1.

Method for Duplex Photoacoustic and Photoacoustic-Echo Imaging

DCA of clinical PA scans for clutter reduction requires ultrasound motion tracking simultaneously with PA image acquisition, i.e., duplex PA and pulse-echo imaging. Unfortunately our ultrasound scanner (see below) did not permit real-time interleaved acquisition of photoacoustic and pulse-echo frames. In order to achieve duplex imaging, a special method was devised, which is described in detail elsewhere.³³ Briefly, the transmit wave for pulse-echo imaging was generated photoacoustically and thus echo ultrasound could be acquired simultaneously with PA data. This made it possible to operate the scanner in receive-only mode for continuous data acquisition, and thus for freehand clinical scanning. For each laser pulse, light that is absorbed by the tissue generates the PA signal in the usual way, whereas light backscattered from the tissue is absorbed at the probe surface, generating a plane PA transient that propagates into the tissue and acts like a pulse-echo transmit wave. Great care is typically taken to optically shield the probe, in order to avoid troublesome echoes generated by this PA transmit wave. In contrast to this, we have taken advantage of the PA transmit wave, placing a gelatine stand-off pad between the probe and the skin surface for temporal separation of PA signals and photoacoustically generated echo signals (henceforth called PA-echo). For a stand-off pad with sound speed $c$ and thickness $D$, the first echoes were generated at a delay $\tau =D/c$ after the laser pulse, corresponding to the acoustic transit time through the stand-off, while PA transients from the tissue were generated at the same time as the laser pulse. PA-echo signals were thus delayed relative to the PA signals by $\tau$, permitting PA imaging to a depth equal to $D$ without interference of PA-echoes. Because PA signal levels from depths greater than $D$ were much weaker than the PA-echo signals, “pure” pulse-echo data could be assumed after the first detection of echoes. Thus, after a single laser pulse, a PA image can be reconstructed from the signals arriving during the time interval $\tau$ to $2\tau$, and a pulse-echo image from those arriving during the interval $2\tau$ to $4\tau$, both covering an imaging depth equal to $D$. The photoacoustic-echo image is, however, of much poorer quality than a conventional pulse-echo image because it uses an unfocused transmit wave and must therefore rely entirely on receive focusing (which is electronic in the plane of a linear array and mechanical in the orthogonal direction).

2.2.

Photoacoustic Imaging Equipment

Commercial medical ultrasound imaging equipment (z.one™, Zonare Medical Systems Inc., USA) was used for PA and pulse-echo signal acquisition. When operated in a research mode, the z.one™ allows the parallel readout of up to 150 μs of signal on each of 64 transducer elements and thus the acquisition of a large data frame after each single laser pulse. Frame sequences of long ($>1\min$) duration can be stored at high frame rate on internal memory for subsequent read out, enabling acquisition of channel data for real-time sequences of freehand clinical scans. A linear-array transducer (L10-5, Zonare) was used for acoustic detection, featuring 5–10 MHz bandwidth ($-3\mathrm{dB}$) and 128 elements at 0.3 mm pitch.

For PA image generation, a Q-switched Nd:YAG laser (ELEN, Italy) was used at 1064-nm wavelength, delivering 70 mJ per pulse with a 7-ns pulse duration. A 10-Hz repetition rate was chosen for all experiments. The laser light was guided via a bifurcated fiber-optic bundle (Fibreoptic, Switzerland) through two line profile converters, one at each side of the linear probe. The output of the profile converters was projected through the gelatine standoff in front of the probe aperture to create two lines of irradiation adjacent to each other, each having 20-mm length and 5-mm width, thus irradiating a patch of 20 mm by 10 mm of the skin surface below the probe. This resulted in $\text{35}\mathrm{mJ}{\mathrm{cm}}^{-2}$ radiant exposure, well below the ANSI safety limit of 1064 nm ($100\mathrm{mJ}{\mathrm{cm}}^{-2}$).

A 20-mm-thick gelatine standoff was used, allowing duplex PA and PA-echo imaging to a tissue depth of 20 mm. Both PA and PA-echo images were reconstructed using frequency-domain algorithms.^5,34 Although the active (64-channel) receive aperture covered only the center half of the probe aperture, the lateral extent of the reconstructed image was chosen equal to the full aperture (38 mm). This was done so that the image would show steering angles pointing outside the active aperture of the probe on the one hand; on the other hand it would be easier to compare photoacoustic-echo images with conventional pulse-echo B-mode images obtained when the scanner was operating in its normal clinical mode.

2.3.

Data Acquisition Procedure

Using an integrated probe containing the ultrasound transducer array, the gelatine standoff and the optical-fiber profile converters,³³ volunteers were scanned with freehand probe guidance. A scan always started with the acquisition of a z.one™ internal B-mode image. Then the z.one™ was remotely switched to passive mode, the laser started, and a sequence of 90 data frames (corresponding to 9 s at 10 Hz) was acquired and stored on the internal memory. While acquiring these 90 frames, the tissue was palpated with lateral reciprocating probe motion while keeping a sticky contact between the probe and the skin surface. This led to cyclic dynamic shear deformation of the tissue under investigation. After switching back to z.one™ internal mode, a second B-mode image was acquired. Comparing the internal B-mode images obtained before and after palpation provided a check for out-of-plane motion during palpation. This information was needed to exclude out-of-plane motion as a mechanism for clutter decorrelation, as explained below. Care was therefore taken during scanning to minimize out-of-plane motion.

2.4.

Motion Tracking

Motion tracking for characterizing the tissue deformation produced by palpation, and its subsequent use for compensating for the deformation in photoacoustic images, was based on the PA-echo image sequence. Accuracy of motion tracking strongly depends on image quality. For ultrasound elastography, tracking is typically based on local block-matching between rigid template and search windows, optimizing a matching criterion (e.g., correlation).^35–38 The template window size must be chosen large enough to obtain a statistically significant correlation maximum. On the other hand, with a large window size, dilatation and shear deformation degrade the correlation between rigid windows. With the poor quality of the PA-echo images (see results), rigid-window block matching proved ineffective for motion tracking. Even a deformable mesh approach^39,40 yielded insufficiently stable results, due to strong decorrelation from image to image. Therefore tracking had to be based on a more statistical approach, where the evolution of the PA-echo image over the whole sequence had to be taken into account.

The human eye is well adapted for performing such an analysis and can track pronounced features throughout the image sequence even if they are for part of the time obstructed by noise. Given the proof-of-principle nature of the present study, tracking was therefore carried out by eye, with the aid of an interactive user interface (written in Matlab™) which allowed trajectories to be assigned to sparsely distributed features throughout the image sequence. A model for displacement versus depth was then fitted to the collection of trajectories, assuming the lateral displacement to be a smooth function of depth. This method yielded stable motion tracking even in the presence of severe clutter, albeit with limited accuracy due to the simplified displacement model and trajectory identification errors.

2.5.

Deformation Compensation

The lateral displacement $\mathrm{\Delta }{x}_n(x,z)$ of an image element relative to its pre-deformation position ($x$, $z$) as a function of frame index $n$ was used for deformation compensation of the PA images, where $x$ is the lateral dimension and $z$ is axial dimension (depth) of the image plane. Displacement was assumed to occur entirely parallel to the image plane, such that tissue deformation at frame index $n$ could be modeled as a two-dimensional coordinate transformation $(x_n,z_n)=[x+\mathrm{\Delta }{x}_n(x,z),z]$, where $x_n$ and $z_n$ are the coordinates of the image element in the $n$th frame. Deformation compensation of the PA image at index $n$ was performed by inverting this coordinate transformation based on spatial interpolation. Because of the limited accuracy of the tracking method, an incoherent averaging approach was adopted by applying DCA to the envelope data as opposed to the radio frequency (RF) data. If RF echo-data were to be employed, and if tracking and compensation were not able to reregister correlated signals to within a fraction of the RF phase, this would on the one hand falsify the correlation analysis for clutter identification, and on the other hand result in loss of signal in the DCA images due to phase cancellation. Although coherent averaging of deformation-compensated RF signals may have been more efficient in cancelling the clutter, the use of incoherent averaging was considered an acceptable limitation, given the proof-of-principle character of the present study.

2.6.

Clutter Analysis

Clutter identification and clutter analysis were based on time-lag correlation of the deformation-compensated photoacoustic (DC-PA) image sequence. The time-lag correlation $C_A(n,n+k)$ of an image region $A$ in average-adjusted DC-PA images $u_n$ and $u_{n+k}$ obtained at time indices $n$ and $n+k$ (with time lag $k$) was defined as:

2.7.

Deformation-Compensated Averaging

Clutter reduction is obtained by applying a moving temporal average to the DC-PA sequence. Signals originating from optical absorbers located in the image plane persist throughout the DC-PA sequence, and are therefore not affected by averaging. Averaging of decorrelating clutter on the other hand reduces the clutter level, and thus improves contrast-to-clutter ratio (CCR), the ratio of the difference between two photoacoustic signal levels and the mean clutter level, which is a measure of image contrast resolution. The impact of DCA on CCR depends on both the clutter decorrelation rate and the averaging-window length. The potential gain in CCR is expected to be roughly the square root of the number of uncorrelated measurements, or the ratio of averaging-window length and decorrelation time constant. The averaging length is limited by the maximum probe displacement and the amount of out-of-plane motion. In our experiments, out-of-plane motion was assessed via high quality ultrasound B-scans acquired at the start and the end of the scanning procedure. This permitted the exclusion of out-of-plane motion as an explanation for clutter decorrelation (see discussion). It was not, however, sufficient to quantify slight out-of-plane motion, which occurs over a large number of frames and limits DCA performance owing to decorrelation of true signal. A fixed averaging length of 20 frames was therefore chosen, based on trial and error, as a value that yielded good CCR gain for all scans. The potential CCR improvement was calculated as the square root of the ratio of this averaging length and the clutter decorrelation rate. Assessing the true resulting contrast-resolution improvement by methods other than by eye has proven a challenging task in in vivo scanning, where the true distribution of PA sources is not known. The final conclusions drawn in this preliminary study are therefore mostly based on the visual appearance of DCA results.

3.1.

Duplex Photoacoustic and Photoacoustic-Echo Imaging, and Motion Tracking

Figure 1 shows example duplex PA and PA-echo ultrasound images of the first author’s neck in transaxial section close to the thyroid (at the position indicated in the pictorial diagram inset in Fig. 2). Figure 1(a) is the PA image reconstructed from the upper part of a frame of ultrasound signals, containing the PA data, and Fig. 1(b) is the echo ultrasound image reconstructed from the lower part of the same frame, containing PA-echo data. The PA-echo image is shown with an overlay of the PA image after thresholding to show only the signals with the highest amplitude from Fig. 1(a). This illustrates the registration between the PA image and the PA-echo image, which is inherent because the two images were generated with the same laser pulse. For comparison, Fig. 1(c) is the conventional pulse-echo image obtained with the internal B-mode facility of the z.one™. The PA-echo images can therefore be used as an intermediate step for registering PA images with the conventional pulse-echo images, based on corresponding features such as those indicated by the white arrows in Fig. 1. In the present study, however, the PA-echo sequence was used for tracking tissue motion during the period of palpation. As explained above, the poor quality of the PA-echo images meant that motion tracking had to be performed visually using envelope data. For this purpose, trajectories such as that shown in Fig. 1(d) (white dotted line) were interactively chosen, each being drawn by eye on a constant-depth M-mode (CM-mode) image, as illustrated. CM-mode images at various depths were used and trajectories were drawn starting at various lateral positions. The displacement model described above was then fitted to the collection of trajectories, which had a regularizing effect, reducing tracking errors. This gave what was believed to be a robust estimate of the displacement as a function of position and frame index, which could then be used to determine the deformation-compensation function that was to be applied to the PA image sequence.

Fig. 1

(a) PA image, reconstructed from the PA part of the signal data. (b) Echo ultrasound image reconstructed from the PA-echo part of the signal data, overlaid with PA data after thresholding, to show only the highest amplitudes from (a). (c) The z.one™ internal B-mode image is shown. The similarity between the z.one™ image and the PA-echo image permits registration of PA features with the anatomical context shown in the z.one™ image (white arrows). (d) Example of manual motion tracking based on the PA-echo image sequence. A tissue trajectory (dotted arrow) is drawn on a constant-depth M-mode (CM mode) image, obtained at 10-mm depth indicated by the horizontal dotted line in (b).

Fig. 2

PA image from the neck, (a) prior to DCA processing, and (b) post DCA. Full white arrows denote image features with clearer appearance after DCA, demonstrating contrast improvement. Empty white arrows denote features which are distinguished only after DCA, demonstrating both contrast improvement and increased imaging depth. The dashed oval encircles a region of generally improved contrast. [Inset in (a): black bar denotes the probe position.]

3.2.

Image Analysis and Evidence for the Existence of Clutter

Image analysis to demonstrate the presence or absence of clutter in various image regions was based on the time-lag correlation of the DC-PA sequence as defined in Eq. 1. Results are shown for the neck scan in Fig. 2, and for the forearm scan in Fig. 3. Analysis regions were chosen at locations where DCA processed scans suggested low level of true PA signal. One region was chosen at depth, where large tissue displacement relative to the skin surface was expected to result in strong clutter decorrelation [see Figs. 4(a) and 5(a), region $A$]. A second region was chosen close to the skin surface, where the rate of change of displacement was low, but the shear rate was large [Figs. 4(a) and 5(a), region $B$]. White overlays to the right of Figs. 4(a) and 5(a) show the lateral-displacement depth profile at the maximum probe displacement, which reaches a maximum of 5 and 8 mm displacement, respectively, at the largest depth where the shear rate appears to be low, and suggests a large shear rate close to the surface. Figures 4(b) and 5(b) show the displacements obtained at the depths of regions $A$, as a function of frame index, illustrating the oscillatory motions of the palpation. Figures 4(c) and 5(c) show the time-lag correlation functions $C(n,n+k)$ evaluated at frame indices $n$ for both large and zero displacement rate [indicated in Figs. 4(b) and 5(b) by vertical lines], for the image regions $A$. Figures 4(d) and 5(d) show the time-lag correlation functions for the same frame indices, but for regions $B$. It may be seen that $C(k)$ strongly depends on the displacement rate. Zero displacement rate results in slow decorrelation [solid line in Figs. 4(c) and 5(c)]. After a small dip at lag $k=1$, caused by uncorrelated noise, correlation stays fairly constant. This demonstrates a high level of systematic signal at zero displacement rate. A steady decrease in correlation for lags larger than five is related to growing displacement rate at those lags. In contrast, at maximum displacement rate, full decorrelation takes place even for small lags ($k<\text{five}$) within region $A$. The interdependence between displacement rate and decorrelation is further substantiated in Figs. 4(b) and 5(b), where the decorrelation rate, defined as the slope of the exponential fit to the correlation function, is plotted together with the displacement curve. Decorrelation rate is smallest for zero displacement rate and largest for large displacement rates.

Fig. 3

Analogous to Fig. 2, but showing results from the left forearm.

Fig. 4

Image statistics for the neck scan. (a) PA image with windows for correlation analysis (white dashed rectangles, $A$ and $B$) and a depth profile of the maximum lateral displacement (white vertical line, and vertical graph, to the right). (b) Displacement at the depth of region $A$ (solid graph, arrow points to left scale), decorrelation rate in window $A$ and window $B$ as function of frame index $n$ (black and grey dots, respectively, arrows point to right scale). (c) Time-lag correlation graphs for window $A$, for zero displacement rate (solid line), starting at the time indicated in (b) (solid vertical line), and for rapid displacement (dashed line), starting at the time indicated in (b) (dashed vertical line). At zero displacement rate, a small dip from $k=0$ to $k=1$ indicates stochastic noise (gray arrow), while after that, correlation persists up to about $k=10$, suggesting systematic signal. In contrast, with rapid displacement, decorrelation is much faster. The spacing between the dotted and the solid line indicates that clutter is present (black arrow). (d) Same as in (c), but for window $B$. Decorrelation is slower owing to slower displacement close to the skin surface.

A possible explanation for decorrelation with tissue deformation could be out-of-plane motion of diffusely distributed PA sources. This would however require significant displacement perpendicular to the image plane, of around 2 mm (thickness of the image plane) per five displacement steps. However, the B-mode images acquired prior to and following the palpation scan indicated only small out-of-plane motion, around 2 mm over the whole 90-frame scan sequence, i.e., about 0.1 mm per five displacement steps. Out-of-plane motion can therefore be excluded as an explanation for decorrelation. A further possible explanation for decorrelation is speckle reformation associated with the spatial rearrangement of subresolution photoacoustic sources upon tissue shear. This effect would be proportional to the shear rate⁴¹ as opposed to the displacement rate. However, the time-lag correlation evaluated in region $B$, where the shear rate was large, indicates very slow decorrelation [Figs. 4(d) and 5(d)], compared to region $A$ where the shear rate was low. The evidence collected therefore strongly supports the hypothesis that clutter was present in the clinical scans.

3.3.

Deformation-Compensated Averaging

For DCA, a moving temporal average was applied to the DC-PA sequence, with window length spanning 20 frames, or 2-s scan time. This interval length was chosen as a compromise between contrast improvement and blurring. Due to inaccurate motion tracking, blurring of real PA features over long time intervals could not be avoided. In addition, out-of-plane motion, while not sufficient to explain short-time decorrelation of the diffuse PA background, still affected the persistence of PA structures on longer time scales. A moving average over 20 frames yielded good results with all scans.

Figures 2 and 3 illustrate the effect of DCA on the two PA image sequences acquired during palpation. They illustrate that DCA improves contrast resolution in the clinical scans. This is mainly noticed around some high amplitude features already visible prior to DCA. DCA led to a more pronounced appearance of these features, and a smoothing of speckled PA background leading to a cleaner appearance of their vicinity (full white arrows). In addition, however, there were features made visible at depth which were not visible prior to DCA (empty white arrows). This indicates that DCA not only improved contrast resolution, but also led to increased imaging depth. Indeed, it can be seen in Figs. 2 and 3 that there is a pronounced increase in the clutter-reducing effect of DCA for depths greater than about 8 mm.

Fig. 5

Image statistics for the forearm. (a), (b), (c), and (d) analogous to Fig. 3.

The contrast-resolution improvement is in agreement with the theoretical expectation described above. Based on the averaging length of 20 frames and a decorrelation rate of 0.5 per frame (taken from Figs. 4 and 5), the maximum potential CCR improvement was expected to be three (the square root of the product of averaging length and clutter-decorrelation rate). Quantitatively assessing the achieved CCR improvement was however challenging, because the “true” spatial distribution of PA sources was not known.