Experimental setup

Figure 1 a shows the setup modeling a flow circuit with acoustic shadows. A gear pump (MGD1000P, RS Components Ltd., Corby, UK) was used to set the flow of degassed water to 4.3 ml/s (on average), monitored with a Coriolis flow meter (CORI-FLOW, Bronkhorst, Veenendaal, The Netherlands). A 1 m-long rigid tube of 6.8 mm inner diameter was used to develop the flow. Thereafter, a vessel phantom was placed in a water-filled Perspex container with an imaging window made of stretch wrap on the side. Imaging was performed using a Vantage 256 system (Verasonics, Kirkland, WA, USA) and the L12-3v probe (Verasonics, set center frequency of 7.3 MHz). SonoVue microbubbles (Bracco Suisse S.A., Geneva, CH) were used as ultrasound contrast agents.

Methods. (a) Experimental setup: degassed water with an ultrasound contrast agent is pumped through a flow phantom, which is imaged using the L12-3v transducer from the side of a water container. (b) Iterative scheme for active attenuation correction (ISAAC). The apodization A _i is obtained based on the error E between the image intensity I _{L , σ} and reference value I _sp . After clipping the integral part E _I , summation with the proportional part E _P and pressure P ₀ gives the requested pressure P _req . At is found by minimizing an error function Y ( A ). After new image acquisition using A _i , the average image intensity ( I _L ) in the marked ROI is smoothed ( I _{L , σ} ). This is the input for the next iteration. The iterative procedure was manually stopped. (c) One-way simulated waveform s(t) . (d) Schematic of the geometry and sampling of s , for a single query point q in the region of interest. For each element, the signal s is sampled at time Δ t _eq , which depends on distances d _eq and z _q (see text). Pressure p _q is obtained by summing the contributions from each element. The dashed lines indicate Δ t _eq <0 and thus s = 0. (e) Example Matrix F (used in eqns [ 3 , 4 ]) for multiple points q . Dividing the result s ( Δ t → e q ) (see [d]) by the distances d → e q , gives one row of this matrix (indicated by the line).

Phantom design

Figure 1 b features an ultrasound image of a typical example of the used flow phantoms (inner diameter 6.8 mm, length 11.4 cm). A stereolithography printer (Form 3+, Formlabs, MA, USA) and Flexible 80a resin (SoS 1760–1800 m/s, attenuation 8.5 dB/mm, at 7.3 MHz [ ]) were used to print the phantoms (no ultraviolet post-curing). To model atherosclerotic plaques, an acoustic shadow was created using a ridge around the 0.5 mm-thick phantom wall. Measurements were performed on nine phantoms with ridge thicknesses, d _ridge , between 0.25 and 2.25 mm, step size of 0.25 mm.

ISAAC: Iterative scheme for active attenuation correction

An iterative algorithm was developed to correct for plaque attenuation. All symbols used in the paper are summarized in Table 1 . Local adjustments to the pressure were achieved by changing the transmit apodization a _e for each transducer element e. A _i at the end of iteration i is a vector composed of the individual a _e , which can be set between 0.2 and 1 on the Verasonics system. The iterations were stopped manually, after which a final adjustment to A led to A _{IS AAC} , which was used for the HFR measurements.

Hardware time-gain compensation was set to its maximum, except for closest to the transducer. An in-house delay-and-sum (DAS) [ , ] implementation was used for image reconstructions. In the reconstruction, no time-gain compensation was applied and a viewing angle [ ] of 90 degrees was used. During ISAAC, the interframe time depends on the time it takes to update the transmit information in the Verasonics system after each iteration. The interframe time was therefore measured using MATLAB’s tic and toc methods.

ISAAC is available as an open-source package on https://github.com/Jelle-Plomp/ISAAC and consists of the following three parts:

Part 1: Obtaining an image profile

ISAAC was initialized by acquiring an image using the apodization A ₀ , where a _e = 0.2 ∀ e ∈[1, N ], with N = 128, the number of transducer elements. A rectangular region of interest (ROI) was selected in the vessel lumen (blue outline in Fig. 1 b). At each iteration, the pixel-wise intensity in the ROI was averaged column-wise to obtain the line profile I _L . As the speckle pattern is not of interest, I _L was smoothed to obtain I _{L , σ} using a Gaussian filter (standard deviation [SD] 1 mm, approximately equal to twice the speckle size [ ]). The median value in the highest 10% of values in I _{L , σ} of the first iteration was empirically selected as a measure for the signal outside the shadow and thus used as the set point value I _sp .

Part 2: Linear acoustic model

A linear acoustic model was used to homogenize the image intensity by adapting A _i . We assumed spherical emission (with 1/ r amplitude decay, where r is the traversed distance) of the acoustic signal s ( t ) by each of the elements, where the amplitude is modulated by a _e . The transmit waveform simulated by the Verasonics system ( Fig. 1 c) was used as s ( t ) for t ∈ [ 0 , t m a x ] , where t m a x = 0.724 μ s is the duration of s . Outside of this interval, s ( t ) = 0.

The pressure in a medium could be obtained using an angular spectrum approach [ ], integrating the contributions from each element over space and time, and subsequently taking the envelope. However, this would be computationally too expensive. Therefore, inspired by DAS, we calculated the pressure at a query point q ( x _q , z _q ) at a single moment in time τ _q , by linear summation of the contributions (each with their own phase). Figure 1 d (left) shows the model geometry for a single point q . At time τ _q , the peak positive pressure of the plane wavefront arrives at q , corresponding to the distance z _q between q and the transducer surface. The phase of the wave from element e at point q at time τ _q is determined by the travel time t _eq over the distance d _eq from e to q . The signals were thus sampled at Δ t e q = τ q − t e q + τ p , where τ _p is the time between t = 0 and the peak of the simulated waveform. A single cycle waveform was used, with the default on-time of 0.67. The on-time is the fraction of a half-cycle during which the tri-state pulse is non-zero, represented by the B parameter of the transmit waveform struct (TW.B) in the Verasonics code. The standard value set to minimize the generation of harmonics by the system is 0.67. This results in a nearly linear relation between a _e and the amplitude p ¯ e of s for element e (see Appendix C ). We thus assumed p ¯ e = γ a e , where γ is an unknown constant relating pressure to apodization. The acoustic pressure p ~ q at q is then ( eqn [1] ):

For simplicity, we continue with a normalized pressure, p q = p ~ q / γ ( eqn [2] ):

Figure 1 d describes the sampling of s for a single q and illustrates the sum in eqn [2] .

For the entire line of queried points (a total of Q points), this generalizes into the system of equations ( eqn [3] ):

The matrix F (see also Fig. 1 e) was defined using the mid-line of the ROI as the query points. The system describes how the (relative) pressure at the query points is mathematically related to the apodization A. See Appendix A for more details on the validity of this model.

Part 3: Obtaining the required pressure and apodization

A proportional integral controller (see Fig. 1 b) adjusted the apodization based on the relative error E between the obtained and requested image intensity I _{L , σ} and I _sp . The proportional error E _p was summed with the integral error E _I and pressure profile P ₀ (from eqn [3] using A ₀ ) to obtain the requested transmit pressure P _req . To prevent large error accumulation, the integral error E _I was first clipped such that P _req – E _p remained approximately between P ₀ and P ₁ (obtained from eqn [3] ) using a e = 1 ∀ e ∈ [ 1 , N ] .

The ‘Minimize’ step in Fig. 1 b used a non-linear least-squares method (MATLAB’s lsqcurvefit ) to find A such that the function Y below was minimized. Minimization was used because F is not necessarily invertible ( eqn [4] ):

A schematic of the implementation of ISAAC in the Verasonics system is presented in Figure 2 .

Schematic representation of iterative scheme for active attenuation correction (ISAAC) as an event sequence (events in light blue) in the Verasonics Vantage system. After acquisition, synchronization and reconstruction events, in the next event (compensation) an external process (a MATLAB function) was called, which contained all sub-processes within the dark blue borders. After the parameters were updated, both a wait and a synchronization event were implemented. Wording corresponds to the Verasonics coding language; the reader is referred to their programming manual for more information. ROI, region of interest.

Partial compensation

ISAAC aims to homogenize I _{L , σ} . However, as attenuation happens both on the forward and return paths, a homogeneous I _{L , σ} means that the acoustic pressure in the shadow region is higher than in the non-shadow region, which could lead to bubble destruction. To assess how much compensation is sufficient to achieve reliable blood flow quantification, HFR measurements were also performed using the following set of apodization profiles A _pc ( ψ ) ( eqn [5] ):

These profiles correspond to partial compensations for signal loss in the ultrasound image.

HFR acquisition, pre-processing, reconstruction and PIV

All HFR recordings consisted of 1000 frames taken at 5000 fps, using unfocused plane waves. The other settings were the same as during ISAAC.

Data analysis was performed in MATLAB (R2023b). First, clutter filtering using singular value decomposition was applied on the raw HFR data [ ] using all 1000 frames as one ensemble, rejecting the first 20 singular values. Next, reconstruction was performed using the in-house DAS algorithm and echoPIV was performed on the resulting IQ data (without log compression). The echoPIV code was adapted from Voorneveld et al. [ ]. Six iterations were used for PIV analysis, with kernel sizes of 64 × 64, 32 × 32 and 16 × 16 pixels, each applied twice for consecutive iterations. Each frame was compared with the fourth next frame, ensuring sufficient displacement between frames. An average of two correlation maps was used for each vector. No temporal smoothing was applied. To illustrate this, the first PIV frame would be constructed based on the average cross-correlation between kernels taken from frame 1 and 5, and 2 and 6, continuing from frame 3 onwards for the second PIV frame. Before averaging the results, two PIV frames were removed at each of the nine transitions between image data buffers (occurring after every 100th image acquisition). This was because at these transitions, the time to the next acquisition was slightly larger (by an unknown amount) than the set frame rate, thereby preventing the recovery of velocity estimates between these frames. The average velocities were obtained based on the remaining 479 vector fields.

Evaluation

Measurements were performed on all nine phantoms. To additionally show the influence of the driving voltage, measurements were performed at 11.2V and 30.4V. The efficacy of the compensation strategy was evaluated by comparing the averaged echoPIV-derived velocity profiles in the region behind the ridge U ¯ B R and the control upstream U ¯ U P , for different d _ridge . The flow profiles at the different locations were assumed to be identical due to the constant tube diameter. In the results, the estimated parabolic flow profiles based on the data of the flow rate sensor are also shown for comparison. However, due to the error margin of this estimate, as well as potential flow disturbances at the connection between the entry tube and the phantom, the upstream velocity profile was chosen as the reference velocity profile. The tube was positioned parallel to the transducer and thus the x -axis (see Fig. 1 b). U ¯ B R and U ¯ U P were computed by averaging the u component of the vectors (flow in the x -direction) over time and over 4 mm in the x -direction. The profiles U ¯ B R were aligned with U ¯ U P using cross-correlation to correct for a perceived shift in z position caused by SoS differences. The relative difference Δ U ¯ was calculated as eqn (6) :

Related posts:

Ambient Pressure Sensitivity of Subharmonic Vibrating Single Microbubbles Comparative Assessment of Real-Time and Offline Short-Lag Spatial Coherence Imaging of Ultrasound Breast Masses Ultrasound-based Velocity Vector Imaging in the Carotid Bifurcation: Repeatability and an In VivoComparison With 4-D Flow MRI Hepatic Portal Venous Perfusion Imaging Using Vessel-Labeling Super-Resolution Ultrasound The Potential for Blood–Brain Barrier Disruption During Transcranial Ultrasound Super-Resolution Imaging Delivering Gd-Labeled IgG Antibodies Into the Mouse Brain Following Focused Ultrasound Treatment

Stay updated, free articles. Join our Telegram channel

Join