Ultrasound-guided photoacoustic image reconstruction: image completion and boundary suppression

Pieter Kruizinga; Frits Mastik; Dion Koeze; Nico de Jong; Antonius F. van der Steen; Gijs van Soest

doi:10.1117/1.JBO.18.9.096017

25 September 2013 Ultrasound-guided photoacoustic image reconstruction: image completion and boundary suppression

Pieter Kruizinga, Frits Mastik, Dion Koeze, Nico de Jong, Antonius F. van der Steen, Gijs van Soest

Author Affiliations +

Journal of Biomedical Optics, Vol. 18, Issue 9, 096017 (September 2013). https://doi.org/10.1117/1.JBO.18.9.096017

Abstract

Photoacoustic (PA) image reconstruction of data acquired with conventional linear arrays suffers from incompleteness due to limited bandwidth and limited view. This problem is compounded by the dominance of boundary signals suppressing other weaker PA sources. In this study, we propose a method that uses a naturally coregistered ultrasound image to enhance the PA image reconstruction, with the dual aim to reconstruct the sources that suffer from incompleteness and reveal weak sources near a strong boundary. In this method, an ultrasound image provides the input for a PA wave field simulation. The simulated PA field can be combined with the measured fields to complete the partial reconstruction or to suppress the dominant boundary signals. Experimental validation of the method was performed with two different phantoms and in vivo data from the lower arm of a healthy volunteer.

1. Introduction

Translating noninvasive photoacoustic (PA) imaging to the clinic likely entails the use of a commercial ultrasound scanner, associated transducer array, and a pulsed laser source.¹^,² This instrument design, providing overlapping fields of view, allows for a natural integration of ultrasonic and PA imaging in a combined modality.³ In this setup, ultrasound (US) can be used to reveal the morphology of the biological structures, while PA provides functional information. This implementation of PA imaging, using conventional ultrasound imaging equipment, is constrained to band-limited signal detection in a limited-view geometry. These two constraints have serious consequences for PA image reconstruction. Here, we consider two consequences: the first is incomplete source reconstruction owing to the limited-view geometry. The second consequence is the suppression of weak PA sources in the vicinity of dominant boundary signals as a result of coherent boundary source build-ups and band-limited detection.

The effect of limited-view detection in PA imaging has been addressed before.⁴^–⁶ The acoustic pressure wave generated by a single PA point source is an outward propagating spherical surface.⁷^,⁸ The source of this pressure wave can therefore be reconstructed by a limited number ( $> 1$ ) of measurements at arbitrary points on this spherical surface. When many PA point sources are packed together, the propagating field generated by this compound source is generally not spherically symmetric. Successful reconstruction of an extended source therefore depends critically on the range of wavenumbers (vectors accounting for frequency and direction of the acoustic waves) that can be detected. In the case of limited-view geometry, only part of the wavenumbers required for full source reconstruction is detected, resulting in a partial source reconstruction. See Fig. 1 for a schematic representation of the limited-view problem. Somewhat related to the limited-view problem is the suppression of weak sources in the proximity of strong boundary signals. As discussed by Guo et al. when PA sources are lined up at a boundary, the individual pressures, which arise after the laser pulse ( $t_{0}$ ), will constructively (due to the similarity in phase and amplitude) build up to a dominant Huygens wavefront.⁹ The sources behind this surface do not necessarily coherently add up to new emerging waves.

Fig. 1

Schematic example of the limited view problem. In (a) a circular pressure distribution at $t_{0}$ with detectors on all four edges of the imaging plane, indicating the full view geometry. The source reconstruction of situation (a) can be found in (b). Subfigure (c) shows the limited view geometry of the same source as in (a), now using detector points along one edge of the imaging plane as is the case for limited view. The reconstruction using limited view geometry is shown in (d).

Interestingly, US imaging suffers less from the limited-view problem and boundary buildup. In US, the acoustic scatterers (sources) exhibit individual time characteristics in their wave emission due to relatively slow propagation of the imaging wave (unlike the optical wave used in PA). The boundary buildup effect and its dominance in PA image reconstruction is a well-known phenomenon. The dominant boundary effect can become increasingly dominant in the case of limited bandwidth detection. Bandwidth limitation is understood as an incomplete acquisition of the full frequency spectrum of the PA signal(s) by the transducer. The part of the PA signal that results from a boundary will most likely exhibit high acoustic pressures, due to the buildup effect, at relatively high frequencies, which is due to the rapid transition in source strength at a boundary interface. In medical imaging, we generally use a transducer with limited bandwidth, sensitive in the few MHz range, acting as a filter receiving the signal. As a result, we are most sensitive to the high-frequency part ( $> 5 MHz$ ) of the PA signal and not so much to the low-frequency part ( $< 5 MHz$ ) of the waves. As an example, an artery with a diameter of 2 mm will, in a cross-sectional view under the condition of homogeneous illumination, emit most acoustic power around 800 kHz ( $F = c_{0} / δ$ , where $F$ denotes the effective frequency, $c_{0}$ the acoustic propagation speed, and $δ$ the effective optical penetration.).¹⁰ The PA signal that we will detect with a commercial 5 MHz transducer is, however, the high pressure–high frequency parts of the PA wave resulting from the 2 mm artery.

To formalize the imaging setting, let us consider the inverse problem for PA image reconstruction where we want to reconstruct $s$ being the initial source pressure at $t_{0}$ , from our measured signals $m$ . We describe propagation of $s$ by the linear operator $M$ (accounting for speed of propagation, dispersion, etc.), and let $R$ be the observation operator that converts the pressure field into a measurement. Putting these together, we obtain Eq. (1) for our measured signals and Eq. (2) for the inverse reconstruction.

Eq. (1)

R [M (s)] = m,

Eq. (2)

s = M^{- 1} [R^{- 1} (m)],

where we assumed the operators

R

and

M

to be invertible.

A real-world measurement cannot acquire all data. For example, there is a finite number of transducers and bandwidth limitation in each experiment. We introduce the approximate observation operator $\tilde{R}$ , which takes these practical measurement issues into account. This implies that our actual measured signals $m_{real}$ resulting from our source $s$ obey

Eq. (3)

m_{real} = \tilde{R} [M (s)] .

Now, reconstruction of the original source requires the inverse of $M$ and $\tilde{R}$ . The inverse of $M$ may be found by making assumptions on the properties of the medium (speed of sound, attenuation, etc.) and by reversing time in the Green’s function describing propagation of $s$ . This propagation model will be called $\tilde{M}$ . Due to the limitations in bandwidth and wavenumber sampling, some frequencies and propagation vectors that are supported by $M (s)$ are not observed in the measurement. These components cannot be used for reconstruction, and as a result, the inverse of $\tilde{R}$ generally does not exist.

An estimate $\hat{s}$ of the source $s$ based on the measurement $m_{real}$ can be written as follows:

Eq. (4)

\hat{s} = {\tilde{M}}^{- 1} [{\tilde{\tilde{R}}}^{- 1} (m_{real})],

Eq. (5)

\hat{s} = {\tilde{M}}^{- 1} ({\overset{\approx}{R}}^{- 1} {\tilde{R} [M (s)]}),

where we introduce

\overset{\approx}{R}

as the modified observation operator, containing our knowledge about the measurement process. It accounts for practical limitations of the measurement and is subject to a regularization that renders it invertible.

In a physical measurement, we will never be able to fully reconstruct the original source, as ${\overset{\approx}{R}}^{- 1}$ is not the inverse of $\tilde{R}$ , despite efforts to approximate $\tilde{R}$ ; bandwidth limitation is a property of a physical measurement. To improve the quality of $\hat{s}$ , we can modify one or more of the three approximated operators. Improving $\tilde{R}$ requires more measurements of the signal in the first place (larger bandwidth, adding sensor points, etc.). Greater accuracy of the estimated physical parameters of the propagation medium translates into an improved $\tilde{M}$ . Finally, more detailed information about the measurement itself (including instrument and object) may improve $\overset{\approx}{R}$ .

Previously, it has been shown that an acoustic velocity map derived from US measurements can be used to improve the reconstruction of PA images.¹¹^,¹² These studies used the fact that both US and PA reconstruction rely on a shared physical parameter, the speed of sound propagation. The improved reconstruction results from a better estimation of the $M$ operator. Another use of US information for improved PA source reconstruction is based on deformation tracking with US imaging to separate real PA sources from background clutter as was shown by Jaeger et al.¹³^,¹⁴ Here ${\overset{\approx}{R}}^{- 1}$ is improved to generate a better estimate $\hat{s}$ by discarding all the components of ${\tilde{M}}^{- 1} [{\overset{\approx}{R}}^{- 1} ({\tilde{m}}_{real})]$ that do not move with the tissue structure as imaged by ultrasound. In this paper, we will address the limited-view and limited-bandwidth reconstruction problems discussed above by augmenting $\tilde{R}$ with simulated detectors, taking advantage of morphological information derived from an US image to improve image reconstruction in cases where limited-view geometry and band limitation applies.

The technique we introduce here was inspired by the observation that large-scale features in a PA image often resemble those seen by US. This suggests that differences in optical and acoustic properties of structures perceived in the images are tissue properties and thus often appear together. In other words, structures perceived with US will also contribute to the PA field, albeit by a different contrast mechanism. In our method, we take an ultrasound image of the same scene as the source input for a PA wave field simulation. To improve source reconstruction hampered by limited view, we use the predicted field to acquire the missing wavenumbers and use them to improve PA image reconstruction. To suppress the dominant boundary signals, we use the predicted field to mask-out the strong boundary signals perceived with the band and view limited transducer.

We explicitly write the limited-view measurement as

Eq. (6)

m_{l} = {\tilde{R}}_{l} [M (s)],

where

{\tilde{R}}_{l}

denotes our limited-view signal acquisition. The reconstructed source based on the measured data

m_{l}

,

{\hat{s}}_{\exp}

, is obtained by

Eq. (7)

{\hat{s}}_{\exp} = {\tilde{M}}^{- 1} [{\overset{\approx}{R}}_{l}^{- 1} (m_{l})] .

We propose to fill in the missing wavenumbers in $m_{l}$ by estimating them based on the tissue morphology.

We construct a simulated source $s_{sim}$ based on the ultrasound image and ${\hat{s}}_{\exp}$ by an algorithm $A$ ,

Eq. (8)

s_{sim} = A ({\hat{s}}_{\exp}; s_{U S}),

where

s_{U S}

is an estimate of

s

, assuming a constant or smoothly varying pressure distribution within its support; its borders are derived from a segmentation of the ultrasound image. With a PA field simulation we can now compute the simulated signals

m_{sim}

.

Eq. (9)

m_{sim} = {\tilde{R}}_{sim} [\tilde{M} (s_{sim})] .

In the simulation, we consider ${\tilde{R}}_{sim}$ as a closed surface signal acquisition as depicted in Figs. 1(a) and 1(b). We choose ${\tilde{R}}_{sim} = {\tilde{R}}_{sim, l} + {\tilde{R}}_{sim, c}$ as the augmented observation operator including the missing wavenumbers not captured in the actual experiment, though retaining the bandwidth limitation of $\tilde{R}$ . The set of complementary observations is represented by ${\tilde{R}}_{sim, c}$ , as observed by virtual detectors on the right, bottom, and left border of the imaging region in Figs. 1(a) and 1(b). Based on $s_{sim}$ , we can obtain the limited-view signals $m_{sim, l}$ and the complementary signals $m_{sim, c}$ by

Eq. (10)

m_{sim, l} = {\tilde{R}}_{(sim, l)} [\tilde{M} (s_{sim})],

Eq. (11)

m_{sim, c} = {\tilde{R}}_{(sim, c)} [\tilde{M} (s_{sim})] .

Minimization of $m_{l} - m_{sim, l}$ may serve to calibrate the procedure $A$ . The above formulation of the PA measurement problem will form the basis for the image reconstruction algorithm we describe in this work.

In order to complement an image with missing wavenumbers due to limited-view geometry, we add the simulated measurement to the actual experimental data to estimate ${\hat{s}}_{cmp}$ .

Eq. (12)

{\hat{s}}_{cmp} = {\tilde{M}}^{- 1} [{\overset{\approx}{R}}^{- 1} (m_{l} + m_{sim, c})] .

For the boundary suppression reconstruction ${\hat{s}}_{bs}$ , we subtract the boundary signals resulting from large-scale structures as seen in the ultrasound image.

Eq. (13)

{\hat{s}}_{bs} = {\tilde{M}}^{- 1} [{\overset{\approx}{R}}^{- 1} (m_{l} - m_{sim, l})] .

As $s_{sim}$ , resulting from the ultrasound segmentation, is homogeneous or smoothly varying, the simulated signals $m_{sim, l}$ will contain only signals that correspond (in case of band limitation) to the boundary of $s_{US}$ . When $m_{sim, l}$ is subtracted from $m_{l}$ , we obtain a set of signals that includes all signals that do not belong to the boundary of $s$ , thus resulting in a boundary supressed image.

Figure 2 shows a schematic overview of the proposed method. We validated this method for both applications by simulations and experimentally using two different vessel mimicking phantoms and in vivo data obtained from the lower arm of a healthy volunteer.

Fig. 2

Schematic overview of the proposed ultrasound-guided photoacoustic (PA) image reconstruction method.

2. Methods and Materials

2.1.

Simulation

The technique of image completion and background subtraction can be tested by simulating PA and US imaging. All simulations in this paper were performed using the freely available acoustic field simulation Matlab toolbox k-Wave¹⁵^,¹⁶ in two spatial dimensions. We defined five unique PA sources:

1. a homogeneous filled circle, mimicking a blood-filled vessel lumen;
2. a homogeneous ring, mimicking the source expected from a vessel wall;
3. a square with a source strength that is discontinuously decaying with depth;
4. a square with continuously decaying source strength and four point sources near the upper boundary;
5. a square with continuously decaying source strength.

All sources, except the last one, have an acoustic impedance difference with respect to the surrounding medium. Acoustic contrast was induced by a random scattering mask with an average impedance of 1.55 MRayl compared to 1.48 MRayl of water. Detectors were placed at all four edges of the simulation grid. Ultrasound plane wave imaging was simulated with the top-row detector points and a short two-cycle 8 MHz transmission pulse. Image reconstruction of the received ultrasound field was done through beamforming in the Fourier domain.¹⁸ An 8 MHz, 60% bandpass filter was applied to mimic band limitation normal for conventional PA acquisition. PA image reconstruction was performed with use of the time reversal k-Wave implementation. More details on this simulation can be found in Table 1.

Table 1

Input parameters for photoacoustic field simulation with k-Wave.

Simulation parameter	Simulation only	Boundary suppression phantom experiment	All other experiments
Grid size	$1024 \times 1024$	$760 \times 860$	$886 \times 906$
Voxel size (μm)	30	37	35
Maximum frequency (MHz)	25	20	22
Speed of sound ( $m / s$ )	1480	1480	1510
Density ( $kg / m^{3}$ )	1000	1000	1000
Acoustic absorption [ $dB / (MHz / cm)$ ]	0.5	0.5	0.5
Simulation duration (s)	137 (av.)	84	129

2.2.

Phantom

For experimental validation of the boundary suppression technique, we prepared a polyvinyl alcohol (PVA) mimicking phantom with three inclusions in the wall. For the preparation of PVA, we dissolved 20 g PVA grains in 200 ml water at 80°C. To induce scattering of light and sound, we added 1 g of SiC (800 mesh) and 1 g of ${SiO}_{2}$ (1 to 10 μm) to 200 ml 9% PVA mixture. Two batches of 40 ml each with different optical absorption coefficient ( $μ_{a}$ ) were prepared by adding Ecoline® 508 Prussian blue water based ink (Royal Talens, Apeldoorn) in two concentrations: 2 ml ( $μ_{a 1}$ ) in solution 1 and 0.5 ml ( $μ_{a 2}$ ) in solution 2. The vessel wall was made using solution 1, and the three inclusions were made using solution 2.¹⁹ The inner and outer diameters of the vessel were, respectively, 6 and 10 mm. The inclusions measured $1 \times 1$ , $1 \times 2$ , and $1 \times 3 mm$ in cross-section. The length of the phantom was 40 mm. The PVA solution was poured into a mold and put through six freeze-thaw cycles to create a stiff gel phantom.

For the image completion experiment, we mimicked the PA source emerging from lumen of a carotid artery with a diameter of 6 mm. For this purpose, we used a slightly modified version of the recipe of Teirlinck et al.²⁰ In short, 18 g of agar was dissolved in a mixture of 0.5l demineralized water, 7 ml 50 wt.% aqueous solution benzalkoniumchloride (ACROS Organics), and 8 ml glycerol (VWR International, Amsterdam, The Netherlands). To induce optical and acoustic scattering, we added 2 g of ${TiO}_{2}$ (Sigma-Aldrich). For optical absorption, we added 1 ml of Prussian blue ink. This mixture was put in a pressure chamber to clear any air trapped during the dissolving process. The mixture was then heated to 90°C. After heating, the dissolved agar mixture was poured into a cylindrical plastic tube with an inner diameter of 6 mm and left to cool down and solidify.

2.3.

Imaging

We used a wavelength tuneable laser (OPOTEK Vibrant B/355-II) generating 5 ns pulses at 10 Hz repetition rate at 650 nm for PA signal generation. The laser was coupled to a one-to-two fiber bundle linear array that could be interfaced with the linear array ultrasound transducer. For PA signal detection and US imaging, we used a 128-element linear array (pitch: 245 μm, pulse-echo $- 6 dB$ bandwidth: 4 to 9 MHz) (Vermon, Tours, France) connected to an open 128-channel US system (Lecoeur Electronique, Chuelles, France), capable of digitizing signals with 80 MHz sampling at 12 bits. For the US image formation, we compounded multiple lines beamformed plane wave insonifications over a $- 7$ to $+ 7$ degree steering angle. Beamforming of both the PA and US signals was done in the Fourier domain.¹⁷^,¹⁸^,²¹ For the image completion phantom experiment, we performed image reconstruction by time reversing the signals. The in vivo images were obtained from the lower arm of a healthy volunteer. The arm was kept stable in water with a temperature of 38°C. In all experiments, we averaged 60 frames to obtain a practical signal-to-noise ratio.

2.4.

Ultrasound-Guided Image Reconstruction

The procedure for ultrasound-guided PA image reconstruction is outlined in Eqs. (1) to (13). Several choices need to be made in the practical implementation of $\tilde{M}$ , ${\tilde{M}}^{- 1}$ , $\tilde{R}$ , ${\overset{\approx}{R}}^{- 1}$ , and particularly $A$ . We present here a minimal scheme to demonstrate the functionality of the reconstruction methods, which can be adapted to specific applications.

All propagating ultrasound fields, whether for simulation or reconstruction, were implemented in the k-Wave toolbox referenced above, used in Matlab (The Mathworks, Natick, MA) version R2012a. This models $\tilde{M}$ (forward) and ${\tilde{M}}^{- 1}$ (backward propagation). Realistic assumptions on the speed of sound, acoustic absorption, dispersion, etc. for this operator have to be made in order to obtain an image. In our experiments, we consider no specific prior knowledge of $\tilde{M}$ , so we assume the medium to be homogeneous and used the same parameters for all reconstructions, save for the speed of sound, which is higher at elevated temperature.

The implementation of the observation operators such as ${\tilde{R}}_{sim}$ and ${\overset{\approx}{R}}^{- 1}$ in our case is straightforward. In all simulations, we defined detector points at every grid cell along four sides of the simulation grid. All detectors record the acoustic wave field for every time step, stored as a pressure trace per detector. The actual measurement is emulated by ${\tilde{R}}_{sim}$ integrating the simulated detector traces over the physical transducer element surfaces. Band limitation, characteristic for a conventional US transducer, was introduced by filtering the received signals with a fourth-order Butterworth bandpass filter (4 to 10 MHz). In the simulation, pressure traces were filtered with the bandpass filter supplied with the k-Wave toolbox. The operator ${\overset{\approx}{R}}^{- 1}$ upsamples the experimentally recorded radio frequency signal $m_{l}$ or $m_{sim}$ to the simulation grid. The same detector points were used to re-emit the acoustic field for time-reversed image reconstruction.

The algorithm $A$ constructs a PA source $s_{sim}$ based on information from US image and an initial PA source estimate ${\hat{s}}_{\exp}$ . It consists of two steps: segmentation of the relevant structures in the US image and assignment of a source pressure to these structures in order to match the experimentally observed PA signal. In the present study, aiming to demonstrate the principle of filling in missing wavenumbers in PA using information from the US image, we opt for manual segmentation (Matlab function roipoly.m). A variety of methods have been proposed for automated US segmentation and border detection.²² These tend to be application specific, however, and their performance depends strongly on US image quality. Development of such an automatic segmentation technique goes beyond the scope of this paper. After selection, the region of interest (ROI) was discretized on the simulation grid and a source pressure $s_{sim}$ was assigned, scaled to the intensity observed in the original PA image reconstruction. A possible small mismatch between the selected ROI and the real PA source may be compensated for through a channel-dependent cross correlation between the experimental data and the simulated data.

The processing of image completion and boundary suppression of the phantom data was done on the raw pre-beamformed experimental signals. In vivo PA signal quality was too poor to allow direct summation and time reversal of $m_{\exp}$ and $m_{sim}$ . For this reason, we performed the ultrasound-guided image reconstructions using beamformed signals.

For image completion, we segmented features in the ultrasound image corresponding to a PA source of interest that could potentially be hampered by incomplete reconstruction due to limited view. Using $s_{sim}$ as defined in Eq. (8), $m_{l}$ as measured, and Eqs. (11) and (12), we computed the ultrasound-guided PA reconstructions.

For boundary suppression, we selected features in the US image corresponding to anatomical structures that may also be recognized in the PA image, attributed a source pressure, and performed PA wave field simulation. A mask was created based on the normalized envelope of the simulated signals, setting to zero all data points in $m_{l}$ at which $m_{sim, l}$ exceeded a threshold value as a practical approximation to Eq. (13). This procedure reduces the sensitivity of ${\hat{s}}_{bs}$ , which is much smaller than $s$ in general, to alignment and calibration errors. The result was filtered and beamformed to obtain the boundary-free PA image.

3. Results

3.1.

Simulation

Image completion and boundary suppression was studied using simulated targets. Figures 3(a) and 3(b) show the simulated plane wave ultrasound scan and the PA pressure map at $t = 0$ , respectively. Image reconstruction of the PA source using the top-row detectors (the same as used for US imaging) can be found in Fig. 3(c). The reconstructed limited-view image obtained from the US estimated PA source $s_{sim, l}$ is displayed in Fig. 3(d). The result obtained with the image completion method is shown in Fig. 3(e): the missing arc segments in the top two objects and the missing sides of the square objects have been filled in with information derived from the US image. The absence of ultrasound contrast in the rightmost square means that these missing wave vectors cannot be filled in. The effect of boundary suppression is shown in Fig. 3(f): the circular structures, lacking internal contrast, disappear. The leftmost square object only shows the internal boundary, which has PA but not US contrast. The small sources just below the top of the middle square are now clearly resolved. Again the square target on the right is unmodified.

Fig. 3

(a) Image obtained from the simulated sources after plane wave US imaging. (b) PA source pressures at zero time. (c) PA source reconstruction by time reversing the signals obtained by detector points at the top row of the simulation grid. (d) PA source reconstruction with data obtained from the estimated PA source. (e) Result obtained with the image completion method. (f) Result obtained with the boundary suppression method.

3.2.

Image Completion Phantom Experiment

Figure 4 presents the experimental realization of PA imaging completion on the cylindrical phantom. The US image of the phantom obtained with plane wave imaging in Fig. 4(a) and the direct PA image ${\hat{s}}_{\exp}$ in Fig. 4(b) are combined to generate the augmented PA source ${\hat{s}}_{cmp}$ shown in Fig. 4(e). The simulated limited-view PA image $\tilde{M} [{\overset{\approx}{R}}^{- 1} (m_{sim, l})]$ in Fig. 4(c) and the full simulation reconstruction $s_{sim}$ in Fig. 4(d) are intermediate steps. Figure 4(f) shows an overlay image of (a) and (e). The sides of the circular PA source are completed by the US-derived information.

Fig. 4

(a) US image and (b) PA image of the cylindrical phantom. (c) Reconstructed PA image from simulated signals obtained in limited view geometry. (d) Reconstructed PA image from simulated signals obtained in full view geometry. (e) Complemented PA source based on data in (b) and (d). (f) Fusion of the ultrasound image shown in (a) and the completed PA image (e).

3.3.

Boundary Suppression Phantom Experiment

The boundary suppression method was demonstrated experimentally on a vessel phantom with low-absorbing inclusions in the wall. A microscopy cross-section of the object is shown in Fig. 5(a). Neither the US [Fig. 5(b)] nor the PA image [Fig. 5(c)] of the phantom has adequate contrast for the inclusions. The PA signal generated by the large-scale structure, which obscures the inclusions, can be simulated based on the US image, according to ${\tilde{M}}^{- 1} [{\overset{\approx}{R}}^{- 1} (m_{sim, l})]$ . It is shown in Fig. 5(d). Subtracting this from the measured PA data ${\hat{s}}_{\exp}$ yields the PA signal of the inclusions ${\hat{s}}_{bs}$ , shown in Fig. 5(e). This image clearly shows the heterogeneous structure of the phantom wall. Note that all images are normalized to their respective maxima.

Fig. 5

(a) Microscopy cross-sectional image of the vessel phantom. (b) US image and (c) PA image of the vessel phantom. (d) Simulated PA image. (e) Boundary-suppressed PA image. (f) Overlay of US and boundary-suppressed PA image.

3.4.

In Vivo Imaging

To demonstrate that this method can also be applied in vivo, we conducted another experiment where we imaged a cross-section of the arm of a volunteer. In this experiment, we identified the skin as a dominant boundary and selected a blood vessel for image completion. Figure 6(a) shows the ultrasound image that was used to identify the water–skin boundary and the blood vessels perceived as darker regions, nonscattering regions in the image. The original PA image can be seen in Fig. 6(b). The skin and the partially imaged blood vessels are clearly identifiable. The image displayed in Fig. 6(d) is the result after boundary suppression in the image domain. The final image (f) is the result of both boundary suppression and image completion in the case of two arteries.

Fig. 6

Ultrasound-guided PA image reconstruction in vivo, imaging the arm of a volunteer. (a) Cross-sectional ultrasound image. (b) PA image of the same cross-section. (c) Zoomed region, containing two blood vessels. (d) PA image after boundary suppression of the signal from the skin. (e) Complementary simulated PA source ${\tilde{M}}^{- 1} ({\overset{\approx}{R}}^{- 1} (m_{sim, c}))$ . (f) Image of the two blood vessels after boundary suppression and image completion.

4. Discussion

In this paper, we propose ultrasound-guided PA image reconstruction. This method uses information derived from US and PA images to simulate the full PA pressure field. We apply this concept to two techniques for US-augmented PA image reconstruction: image completion and boundary suppression. Both techniques, under certain conditions, are applicable to in vivo PA imaging, as we demonstrated in this paper.

For the image completion technique, we used the simulated signals to complete the experimental data with wavenumbers that were not captured due to limited view, i.e., those perpendicular to the transducer. Image completion may be especially useful for the interpretation of images containing blood vessels, as the relevant anatomy can be more readily recognized. This also facilitates automated image analysis.²³ Blood vessels are of interest to PA imaging because they strongly absorb light and may be used for monitoring of tissue oxygenation and local metabolism.

Our approach for image completion shows similarities to the technique called iterative time reversal, included in k-Wave. Both methods rely on the simulation of missing wavenumbers at virtual detectors. The main difference is that iterative time reversal requires the initial acquisition of the transverse k-vectors by a second, perpendicular, transducer to start the optimization. We propose, in contrast, to generate an estimate of $m_{sim, c}$ in a single step based on independent image data.

For the boundary suppression technique, we used the simulated signals to suppress the experimental signals that emerged from strong boundary interfaces. Boundary suppression is applicable to a variety of cases where a small PA source is located close to a dominant interface, which would make the object of interest difficult to identify. One aspect of the technique that contributes to improved visualization is that the signal originating from the object of interest now constitutes the full dynamic range of the image. The US-based simulation of the PA background selectively removes the signal originating from the large-scale structure. This cannot be achieved by the simple application of a threshold or mask. Exact colocation also enables the separation of multiple sources (clutter and object) within the signal envelope, relating to the limited bandwidth situation discussed in the Introduction.

We demonstrate US augmented PA reconstruction in two dimensions in the present paper. The methods can easily be generalized to three spatial dimensions. The computational tools for this extension are available, and the assumptions we make hold equally valid.

4.1.

Assumptions and Limitations

The ultrasound-guided image reconstruction method proposed in this paper relies strongly on the assumption that structures that appear in US also contribute to the PA field, as was illustrated by the simulation experiment (Fig. 3). This assumption is not necessarily valid in all imaging settings, and proper application of the method requires adequate judgment of the user. Still, in many cases, a large-scale similarity between US and PA images of a scene is observed. This can be explained by the fact that image contrast in both modalities depends on tissue-specific properties. Optical absorption is the main contrast mechanism for PA imaging, but its variation may coincide with changes in echogenicity, density, and speed of sound, which are responsible for US contrast. The complementarity between the modalities is reinforced by the fact that the image reconstruction is almost identical.

The quality of the US image has a strong impact on the accuracy of the PA source reconstruction. The simulated source $s_{sim}$ depends on the quality of the input $s_{US}$ . High-quality US images of a known anatomy may allow automated segmentation by a customized computer-based algorithm $A$ . Likewise, the initial estimate ${\hat{s}}_{\exp}$ requires accurate knowledge of the medium and imaging system parameters in order to compute ${\tilde{R}}^{- 1}$ , ${\overset{\approx}{R}}^{- 1}$ , and ${\tilde{M}}^{- 1}$ . The data presented in the present study could benefit from better signal-to-noise ratio and resolution in the US data, offered by high-end commercial imaging systems. More realistic modeling of transducer characteristics (directionality, impulse response) will improve $m_{sim, l}$ . Sophisticated signal and image processing techniques can be applied to ensure optimal registration between experimental and simulated data.

We require a constant or smoothly varying pressure distribution within the support of $s_{sim}$ on $t = 0$ . Strong variations in the source pressure without acoustic contrast will generate a spurious PA signal that may or may not be of interest, but is not affected by our reconstruction. An example is the horizontal line in the bottom-left object in Fig. 3.

5. Conclusions

We propose a new method for image completion and boundary suppression in PA imaging. We show that morphological information derived from an US image can be used to predict the PA wave field by simulation. We use these modeled signals to either complete the initial PA reconstruction or subtract features that are common in both images, thereby enhancing the specific contrast offered by normal PA. This goes beyond simply overlaying the data obtained by US and PA separately. Rather, we use information from one modality (US) to simulate the expected signal in the other (PA). The simulated measurement in the second modality can be used to constrain or augment the—often ill-posed—problem of image reconstruction. We demonstrated the validity of the method with phantom and in vivo measurements.

Acknowledgments

The authors would like to thank Victor Stoev for his assistance during the experimental work. Also thanks to Geert Springeling and Michiel Manten for constructing excellent setups and phantom moulds to perform the experiments mentioned in this paper.

References

1.

C. Kimet al., “Deeply penetrating in vivo photoacoustic imaging using a clinical ultrasound array system,” Biomed. Opt. Express, 1 (1), 278 –284 (2010). http://dx.doi.org/10.1364/BOE.1.000278 BOEICL 2156-7085 Google Scholar

2.

P. Beard, “Biomedical photoacoustic imaging,” Interface Focus, 1 (4), 602 –631 (2011). http://dx.doi.org/10.1098/rsfs.2011.0028 Google Scholar

3.

S. Emelianovet al., “Synergy and applications of combined ultrasound, elasticity, and photoacoustic imaging,” in Proc. of the 2006 IEEE Int. Ultrasonics Symp., 405 –415 (2006). Google Scholar

4.

Y. Xuet al., “Reconstructions in limited-view thermoacoustic tomography,” Med. Phys., 31 (4), 724 –733 (2004). http://dx.doi.org/10.1118/1.1644531 MPHYA6 0094-2405 Google Scholar

5.

M. XuL. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 (4), 041101 (2006). http://dx.doi.org/10.1063/1.2195024 RSINAK 0034-6748 Google Scholar

6.

G. Paltaufet al., “Experimental evaluation of reconstruction algorithms for limited view photoacoustic tomography with line detectors,” Inverse Probl., 23 (6), S81 –S94 (2007). http://dx.doi.org/10.1088/0266-5611/23/6/S07 INPEEY 0266-5611 Google Scholar

7.

G. DieboldT. SunM. Khan, “Photoacoustic monopole radiation in one, two, and three dimensions,” Phys. Rev. Lett., 67 (24), 3384 –3387 (1991). http://dx.doi.org/10.1103/PhysRevLett.67.3384 PRLTAO 0031-9007 Google Scholar

8.

B. CoxP. Beard, “Fast calculation of pulsed photoacoustic fields in fluids using k-space methods,” J. Acoust. Soc. Am., 117 3616 –3627 (2005). http://dx.doi.org/10.1121/1.1920227 JASMAN 0001-4966 Google Scholar

9.

Z. GuoL. LiL. V. Wang, “On the speckle-free nature of photoacoustic tomography,” Med. Phys., 36 (9), 4084 –4088 (2009). http://dx.doi.org/10.1118/1.3187231 MPHYA6 0094-2405 Google Scholar

10.

A. OraevskyS. JacquesF. Tittel, “Measurement of tissue optical properties by time-resolved detection of laser-induced transient stress,” Appl. Opt., 36 (1), 402 –415 (1997). http://dx.doi.org/10.1364/AO.36.000402 APOPAI 0003-6935 Google Scholar

11.

X. JinL. Wang, “Thermoacoustic tomography with correction for acoustic speed variations,” Phys. Med. Biol., 51 (24), 6437 –6448 (2006). http://dx.doi.org/10.1088/0031-9155/51/24/010 PHMBA7 0031-9155 Google Scholar

12.

S. Manoharet al., “Concomitant speed-of-sound tomography in photoacoustic imaging,” Appl. Phys. Lett., 91 (13), 131911 (2007). http://dx.doi.org/10.1063/1.2789689 APPLAB 0003-6951 Google Scholar

13.

M. Jaegeret al., “Reduction of background in optoacoustic image sequences obtained under tissue deformation,” J. Biomed. Opt., 14 (5), 054011 (2009). http://dx.doi.org/10.1117/1.3227038 JBOPFO 1083-3668 Google Scholar

14.

M. Jaegeret al., “Improved contrast deep optoacoustic imaging using displacement-compensated averaging: breast tumour phantom studies,” Phys. Med. Biol., 56 (18), 5889 –5901 (2011). http://dx.doi.org/10.1088/0031-9155/56/18/008 PHMBA7 0031-9155 Google Scholar

15.

. http://www.k-wave.org/ Google Scholar

16.

B. Coxet al., “k-space propagation models for acoustically heterogeneous media: application to biomedical photoacoustics,” J. Acoust. Soc. Am., 121 3453 –3464 (2007). http://dx.doi.org/10.1121/1.2717409 JASMAN 0001-4966 Google Scholar

17.

B. TreebyB. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., 15 (2), 021314 (2010). http://dx.doi.org/10.1117/1.3360308 JBOPFO 1083-3668 Google Scholar

18.

P. Kruizingaet al., “Plane-wave ultrasound beamforming using a nonuniform fast Fourier transform,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 59 (12), 2684 –2691 (2012). http://dx.doi.org/10.1109/TUFFC.2012.2509 ITUCER 0885-3010 Google Scholar

19.

A. Kharineet al., “Poly (vinyl alcohol) gels for use as tissue phantoms in photoacoustic mammography,” Phys. Med. Biol., 48 (3), 357 –370 (2003). http://dx.doi.org/10.1088/0031-9155/48/3/306 PHMBA7 0031-9155 Google Scholar

20.

C. J. Teirlincket al., “Development of an example flow test object and comparison of five of these test objects, constructed in various laboratories,” Ultrasonics, 36 (1), 653 –660 (1998). http://dx.doi.org/10.1016/S0041-624X(97)00150-9 ULTRA3 0041-624X Google Scholar

21.

J. Lu, “2D and 3D high frame rate imaging with limited diffraction beams,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 44 (4), 839 –856 (1997). http://dx.doi.org/10.1109/58.655200 ITUCER 0885-3010 Google Scholar

22.

J. A. NobleD. Boukerroui, “Ultrasound image segmentation: a survey,” IEEE Trans. Med. Imaging, 25 (8), 987 –1010 (2006). http://dx.doi.org/10.1109/TMI.2006.877092 ITMID4 0278-0062 Google Scholar

23.

T. OrugantiJ. G. LauferB. E. Treeby, “Vessel filtering of photoacoustic images,” Proc. SPIE, 8581 85811W (2013). http://dx.doi.org/10.1117/12.2005988 PSISDG 0277-786X Google Scholar

Citation Download Citation

Pieter Kruizinga, Frits Mastik, Dion Koeze, Nico de Jong, Antonius F. van der Steen, and Gijs van Soest "Ultrasound-guided photoacoustic image reconstruction: image completion and boundary suppression," Journal of Biomedical Optics 18(9), 096017 (25 September 2013). https://doi.org/10.1117/1.JBO.18.9.096017

Published: 25 September 2013

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 6 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Image restoration

Ultrasonography

Photoacoustic spectroscopy

Image segmentation

Acoustics

Sensors

In vivo imaging

1.

Introduction

Fig. 1

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Fig. 2

2.

Methods and Materials

2.1.

Simulation

Table 1

2.2.

Phantom

2.3.

Imaging

2.4.

Ultrasound-Guided Image Reconstruction

3.

Results

3.1.

Simulation

Fig. 3

3.2.

Image Completion Phantom Experiment

Fig. 4

3.3.

Boundary Suppression Phantom Experiment

Fig. 5

3.4.

In Vivo Imaging

Fig. 6

4.

Discussion

4.1.

Assumptions and Limitations

5.

Conclusions

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years