Ultrasound contrast enhanced imaging has seen widespread uptake in research and clinical diagnostic imaging. This includes applications such as vector flow imaging, functional ultrasound and super-resolution ulm. All of these require testing and validation during development of new algorithms with ground truth data. In this work we present a comprehensive simulation platform bff that generates contrast enhanced ultrasound images in vascular tree geometries with realistic flow characteristics and validation algorithms for ulm. bff allows complex micro-vascular network generation of random and user-defined vascular networks. Blood flow is simulated with a fast cfd solver and allows arbitrary input and output positions and custom pressures. The acoustic field simulation is combined with non-linear mb dynamics and simulates a range of point spread functions based on user-defined mb characteristics. The validation combines both binary and quantitative metrics. BFF's capacity to generate and validate user-defined networks is demonstrated through its implementation in the ultra-sr Challenge at the ius 2022 of the ieee.
The ability to produce ulm images, and the availability of a ground truth in localisation and tracking enables objective and quantitative evaluation of the large number of localisation and tracking algorithms developed in the field. bff can also benefit deep learning based methods by automatically generating datasets for training. bff is a fully comprehensive simulation platform for testing and validation of novel ulm techniques and is open source.

§ INTRODUCTION

Contrast enhanced ultrasound has seen a widespread uptake in research and clinical diagnostic imaging. This includes emerging imaging modalities such as vector flow imaging, super-resolution ulm and functional ultrasound. ulm is of particular interest as it is the only imaging modality that can visualize vascular structures and flow velocity information at microscopic resolution in deep tissue in vivo <cit.>. This makes ultrasound valuable for clinical applications where tissue microvascular flow is of interest. For example, cancer research has shown that increased angiogenesis is an early event in the development of tumours <cit.>. In the case of atherosclerosis, which can lead to strokes and cardiac arrest, angiogenesis is the predominant form of neovascularization <cit.> and ulm could have the potential to be used for screening.

Most ultrasound contrast agents are highly compressible gas mb surrounded by a lipid monolayer. mb behave as scatterers due to the difference in acoustic impedance between the gas and the surrounding blood. But mb also respond non-linearly to the incident pressure wave due to their compressible gas component. Taking advantage of their non-linear response, contrast sequences such as pi or am allow their easy separation from regular tissue. ulm localizes isolated mb and tracks them over many frames as they flow through the vascular network. Through the localization ulm creates sub wavelength resolution vessel representations and flow velocity maps.

The requirement for isolated bubbles poses a sparsity constraint that generates a series of limitations on the mb concentration that can be used and the total acquisition time that is needed. Improving the performance of ulm towards the goal of real-time ulm is desired. This puts the focus on providing better bubble separation and increasing the localization precision. Efforts have been made to increase the speed of ulm for example by localizing mb at higher densities through deconvolution and multi-feature tracking <cit.>, through sparsity-based methods <cit.>, by splitting the bubble signal in the 3D Fourier domain <cit.> or through the use of phase change contrast agents <cit.>. However, to evaluate the large number of localisation and tracking algorithms properly, realistic datasets with ground truth are required. For ulm, such datasets should comprise realistic mb psf that coherently interact when overlapped, and the flow should reflect the structures and behavior of the clinical end-application. Datasets also need to be large and contain tens of thousands of mb to aid in the design and network training of the increasingly popular deep learning based ulm methods <cit.>.

There are many tools for simulating the propagation of ultrasound wave fields. Two prominent examples are Field II <cit.> and k-Wave <cit.>. There are also equations that model the behavior of mb, and commercial tools for computing flow such as StarCCM <cit.>. But they are not specifically designed for ulm and they are difficult to combine. Random bivariate Gaussian curves have been used as psf for training cnn models <cit.>. But these are limited due to the lack of sidelobes or bubble tails, and the interaction between overlapping psf is not coherent. Womersley flow has been used to perform fast simulation of flow on large vessels <cit.>. But this method yields results for a single cylindrical section of a vessel, it is not a realistic model for microvasculature, and is computationally expensive to solve. There have been several attempts to model the fluctuation of a mb's radius when exposed to an incident pressure field. The most basic form is the Rayleigh-Plesset equation <cit.> that models a gas bubble in an infinite pool of liquid, with the limitation of its assumption that the bubble boundary moves slower than the speed of sound. Further models do consider the acoustic emission of the bubble <cit.>, but not the behavior of the lipid monolayer. ULM normally insonifies mb with pressure high enough to produce significant non-linear response. A more recent model added shell buckling behavior with compression and rupture on large oscillations, which considers the physical properties of the lipid monolayer and its possible buckled, elastic and broken states <cit.>. A more general approach was demonstrated by combining k-Wave and the Marmottant model <cit.>. This method included reflections and aberrations produced by the non-linear propagation on the tissue, but it is highly demanding on computing power and memory, and requires a spatial grid to be defined which affects the accuracy of the mb locations. None of the efforts so far are designed for large dataset generation, and comprise all four essential components: bubble response, acoustic fields, flow behavior and validation.

To address the shortcomings we created a fully comprehensive simulation platform for ulm, called bff, that incorporates arbitrary micro-vessel network generation with a quick cfd solver and linear acoustic propagation with non-linear mb dynamics and subsequent binary and quantitative evaluation for ulm algorithms. In the following we will describe how bff uses a custom tool to generate organic microvessel structures, how the Hagen–Poiseuille model for quick cfd simulation is calculated, how Field II as the acoustic propagation backend is combined with the Marmottant ode to simulate the response of the mb and how the evaluation of ulm is achieved by means of the implementation of bff in the Ultrasound Localisation and TRacking Algorithms for Super Resolution (ULTRA-SR) Challenge at the ius 2022 of the ieee.

§ METHODS

§.§ Microvascular Flow

To model microvascular structures we first need to briefly define microvascular flow. The Reynolds number describes the ratio of inertial forces to viscous forces. In small vessels of the microvasculature the Reynolds number is very low and radial components of flow are zero considering that flow is laminar for $Re < 2300$

\begin{equation}\label{eqn:reynolds}
Re = \frac{u D}{\mu}
\end{equation}

Where $u$ is the flow velocity, $D$ the vessel diameter and $\mu$ the kinematic viscosity of blood. Furthermore, such flow is quasi steady, fully developed and vessels segments are assumed to be rigid and straight. Subsequently, the axial momentum equation of the flow field derived from the Navier Stokes equation for incompressible Newtonian fluids reduces to

\begin{equation} \label{eqn:axial_momentum}
\frac{1}{r}\frac{\delta}{\delta r}\left( r \frac{\delta u_z}{\delta r} \right) = - \frac{1}{\mu}\frac{\delta p}{\delta z}
\end{equation}

where $\delta p/\delta z$ is the pressure as a function of the axial coordinate $z$ and $r$ is the radial coordinate. Imposing a no-slip condition on the wall ($u=0$ at $r=R$) a parabolic shape of the velocity profile can be determined

\begin{equation} \label{eqn:parabolic_profile}
u_z(r) = u_{z,max} \left( 1 - \frac{r^{2}}{R^{2}}\right)
\end{equation}

where $u_{z,max}$ is the maximum velocity, which is twice the mean velocity ($u_{z,mean} = 1/2 u_{z,max}$) for a parabolic (cross-sectional) flow profile

\begin{equation}\label{eqn:max_velocity}
u_{z, max} = \frac{R^{2}}{4\mu}\frac{\Delta p}{l}
\end{equation}

The volume flow in a vessel with finite length $l$ can be calculated with the Hagen–Poiseuille law using the hydrodynamic resistance $\xi$ and the pressure difference $\Delta p$.

\begin{equation} \label{eqn:hagen_poiseuile}
Q = \frac{\pi r^{4}}{8\mu l} \Delta p = \frac{\Delta p}{\xi}
\end{equation}

These basic equations of flow and pressure will be used in the network solver as a system of linear equations.

§.§ Network Generation

bff implements an organic vessel tree generation tool based broadly in a recurrent process of vasculogenesis and sprouting angiogenesis. Vessel structures are approximated as weighted directed graphs with nodes and edges. Each edge represents a cylindrical section of a vessel, with characteristics of radius, length and orientation; and each node represents the connection between two edges. The final output of this tool is a randomized binary tree structures that can be fully customized to mimic what is seen in vivo.

The process of network creation starts with an initial state defined by: initial position, orientation, radius, and recursion level. A vessel is created by iteratively appending new edges. Each new edge will have new randomized properties based on the current state: the rotation both in elevation and asimuth can be tuned to different vessel tortuosities, and the radius can be changed as new segements are added, creating a radius decay along the vessel. Random locations along the created vessel are chosen as 'sprouts' or locations for bifurcations, and the process is run recursively, creating new vessels at those locations. The whole process is highly customizable. This is achieved by defining parameters as functions of the current status of the generation algorithm: $param=f(n, d, e1, e2, r, lvl)$ with:

* $n$: the current node

* $d, e1, e2$: the current edge orientation (orthonormal reference system)

* $r$: the current radius

* $lvl$: the current recursion level

The network generation can be constraint by a global set of parameter as shown in Table <ref>. Organs usually have specific shapes, and that shape limits the extent of its vessel structure. For example the $inside\_f$ function is used for creating specific organ shapes, by limiting the generation of vessels only to positions inside the shape.

Parameter     Description

edge_step_f     step size in meters

inside_f     when to stop generating

rot_f     new orientation to use on next section of current vessel

r_decay_f     new radius to use on next section of current vessel

bif_occurs_f     whether or not a bifurcation occurs at this place

bif_r_decay_f     new radius to use on new branch vessel

bif_rot_f     new orientation to use on new branch vessel

The global network generation parameters can be used to constrain the network generation. For example, edge_step_f describes the spatial frequency of edges and is inside_f determines the outlining 3D shape of the network.

§.§ Network Solver

The flow and pressure in a vessel network can be solved for every position as a system of linear equations. To construct the equations, the node-edge incidence matrix, input and output pressures, and the hydrodynamic conductance of all the edges are required. In addition, the following assumptions are made:

* fluid is incompressible

* inputs and outputs are hanging nodes

* a hanging node is connected to only one edge

Vessel networks are represented as a weighted directed graph, containing nodes and edges illustrate in Figure <ref>.

0.7 !

vertex/.style=circle, fill=black!25, minimum size=20pt,
highlighted_vertex/.style=vertex, fill=red!24,
edge/.style=draw, thick, ->,
highlighted_edge/.style=draw, line width=5pt, -, red!10,

background, main

[scale=1.8, auto,swap]

[vertex] () at $\name$;

/ /in
[edge] () -> node[weight] $\weight$ ();

node[highlighted_vertex] at () $\vertex$;

/ in
[highlighted_edge] (.center) – (.center);

Simple sample network to illustrate the relation between node (a-f) and edges (1-7). Nodes connected to only one edge are hanging nodes. Their corresponding edge is highlighted.

The incidence matrix I encodes the structure of the vessel network. The matrix I can be used to express a system of linear equations that relate the edge pressure differences and the node pressures, and another system of linear equations that relates the current at each node. The edge-node incidence matrix for the sample network is

\begin{equation*}
\mathbf{I}=
\begin{blockarray}{r ccccccc}
\begin{block}{r ccccccc}
& a & b & c & d & e & f & g \\
\end{block}
\begin{block}{r [ccccccc]}
e_1 &  1  & -1  &  0  &  0  &  0  &  0  &  0 \\
e_2 &  0  &  1  & -1  &  0  &  0  &  0  &  0 \\
e_3 &  0  &  1  &  0  & -1  &  0  &  0  &  0 \\
e_4 &  0  &  0  &  1  & -1  &  0  &  0  &  0 \\
e_5 &  0  &  0  &  1  &  0  & -1  &  0  &  0 \\
e_6 &  0  &  0  &  1  &  0  &  0  & -1  &  0 \\
e_7 &  0  &  0  &  0  &  1  &  0  &  0  & -1 \\
\end{block}
\end{blockarray}
\end{equation*}

The pressure at the hanging nodes is known from the start, but the pressure at non-hanging nodes are unknowns that need to be solved. The incidence matrix is separated into two: $I_h$ and $I_{nh}$, containing the columns of $I$ corresponding to hanging nodes and non-hanging nodes:

\begin{equation*}
\mathbf{I_h}=
\begin{blockarray}{r cccc}
\begin{block}{r cccc}
& a & e & f & g \\
\end{block}
\begin{block}{r [cccc]}
e_1 &  1  &  0  &  0  &  0  \\
e_2 &  0  &  0  &  0  &  0  \\
e_3 &  0  &  0  &  0  &  0  \\
e_4 &  0  &  0  &  0  &  0  \\
e_5 &  0  &  1  &  0  &  0  \\
e_6 &  0  &  0  &  1  &  0  \\
e_7 &  0  &  0  &  0  &  1  \\
\end{block}
\end{blockarray}
\mathbf{I_{nh}}=
\begin{blockarray}{r ccc}
\begin{block}{r ccc}
& b & c & d \\
\end{block}
\begin{block}{r [ccc]}
e_1 & -1  &  0  &  0  \\
e_2 &  1  & -1  &  0  \\
e_3 &  1  &  0  & -1  \\
e_4 &  0  &  1  & -1  \\
e_5 &  0  &  1  &  0  \\
e_6 &  0  &  1  &  0  \\
e_7 &  0  &  0  &  1  \\
\end{block}
\end{blockarray}
\end{equation*}

The pressure difference across edges can be expressed as:

\begin{equation} \label{eqn:edge_node_pressure}
P_e = I_h P_0 - I_{nh} P_n
\end{equation}

with $P_0$ the known pressure at each hanging node, and $P_n$ the unknown pressure at non-hanging nodes. Given the assumption that the fluid is incompressible, the sum of flows in each non-hanging node must be zero. This can be expressed as a system of linear equations:

\begin{equation} \label{eqn:flow_node}
I_{nh}^{T} Q_e = 0
\end{equation}

Equation <ref> relates the pressure difference at the ends of an edge to the volume flow going through it. This can be expanded into a system of equations for each edge as:

\begin{equation} \label{eqn:edge_pressure_flow}
Q_e = C P_e
\end{equation}

with $P_e$, edge pressure difference vector; $Q_e$, Edge volumetric flow vector; and $C$, Edge Flow conductance Matrix

The edge conductance matrix is a diagonal matrix whose values be obtained by calculating the flow conductance at each edge, using its radius and the viscosity of the liquid on <ref>:

\begin{equation*}
\mathbf{C}=
\begin{blockarray}{r ccccccc}
\begin{block}{r ccccccc}
& e_1 & e_2 & e_3 & e_4 & e_5 & e_6 & e_7 \\
\end{block}
\begin{block}{r [ccccccc]}
e_1 &   \frac{1}{\xi_1}  &  0  &  0  &  0  &  0  &  0  &  0 \\
e_2 &  0  &   \frac{1}{\xi_2}  &  0  &  0  &  0  &  0  &  0 \\
e_3 &  0  &  0  &   \frac{1}{\xi_3}  &  0  &  0  &  0  &  0 \\
e_4 &  0  &  0  &  0  &   \frac{1}{\xi_4}  &  0  &  0  &  0 \\
e_5 &  0  &  0  &  0  &  0  &   \frac{1}{\xi_5}  &  0  &  0 \\
e_6 &  0  &  0  &  0  &  0  &  0  &   \frac{1}{\xi_6}  &  0 \\
e_7 &  0  &  0  &  0  &  0  &  0  &  0  &   \frac{1}{\xi_7} \\
\end{block}
\end{blockarray}
\end{equation*}

Given <ref>, <ref>, and <ref> we can replace <ref> in <ref> to obtain:

\begin{equation}
I_{nh}^{T} C P_e = 0
\end{equation}

by replacing <ref> we get:

\begin{equation}
I_{nh}^{T} C  (I_h P_0 - I_{nh} P_n) = 0
\end{equation}

\begin{equation}
I_{nh}^{T} C I_h P_0 - I_{nh}^{T} C I_{nh} P_n = 0
\end{equation}

\begin{equation}
I_{nh}^{T} C I_{nh} P_n = I_{nh}^{T} C I_h P_0
\end{equation}

which can be expressed as a linear system of equations in the form of $M x = b$ with:

* $M = I_{nh}^{T} C I_{nh}$

* $x = P_n$

* $b = I_{nh}^{T} C I_h P_0$

This will solve for the pressure at each non-hanging node, and, using equations <ref>, <ref>, the volumetric flow at each edge can be derived. Note, a hanging node can be placed anywhere within the network.

§.§ Track generation

The bubble positions and their tracks through the network implemented through a fast track generation algorithm based three conditions:

* the probability of a bubble taking a path at a bifurcation is proportional to the flow in that path.

* bubbles move only in a streamline (they maintain their relative radial position through the whole network according to the rules of laminar flow).

* the number of bubbles is much greater than the number of all possible tracks

All possible tracks or root-to-leaf paths can be extracted from the network by computing all the combinations of paths at each vessel bifurcation. At each bifurcation, a bubble can randomly go through one of the paths. The probability of choosing one is calculated based on flow conservation and fluid incompressibility. The volume flow rate of the inlet $Q_1$ is the sum of its outlet $Q_2$ and $Q_3$:

\begin{equation}
Q_1 = Q_2 + Q_3
\end{equation}

The fraction of the incoming volume that will go to branch two will then be $\frac{Q_2}{Q_1}$, so we can approximate the probability of taking that branch as:

\begin{equation}
p = \frac{Q_2}{Q_1}
\end{equation}

The probability of a bubble traversing a whole track is calculated as the product of all bifurcation probabilities in a track:

\begin{equation} \label{eqn:trac_probability}
p_{track} = \prod_{i \in T}^{} \frac{edge_i.Q}{edge_{i-1}.Q}
\end{equation}

To generate the ground truth information, any number of bubbles are moved with the network flow. For each bubble, a track is randomly chosen using the track probability calculated with <ref>. A starting position is chosen for the bubbles comprising both radial, angular and axial position in the starting vessel. From the starting position, a particle will move through edges at a constant speed of $edge.vel * ( 1-r^2 )$, with $edge.vel$ the (max) velocity at the center of the vessel, and $r$ is the fractional radius of the streamtube the particle is traveling in. The particle position is updated using its velocity and the $\Delta t$. Once the particle reaches the end of the current edge, the overshoot is corrected and the starting position for the next edge defined:

\begin{equation} \label{eqn:overshoot}
pos = next\_edge.start + overshoot * \frac{next\_edge.vel}{cur\_edge.vel}
\end{equation}

where $next\_edge.start$ is the starting location of the next edge, $next\_edge.vel$ is its velocity and $cur\_edge.vel$ the current edge velocity. This process is repeated until the bubble reaches the end of the track. The result of this simulation is a table of events with in the format $[frame, bubble_id, x, y, z]$. Additional parameters such as particle velocity or tube radius can be added to the table as ground truth if desired. For the final output table, all the tables from different bubbles are stacked vertically and sorted by frame number.

§.§ Non-linear Bubble Simulation

bff can solve for any ode to obtain the response of a mb, as a default, we chose to use a modified Rayleigh–Plesset equation as follows <cit.>:

\begin{multline}  \label{eqn:marmottant_ode}
\rho_l (R \ddot{R}+ \frac{3}{2}\dot{R}^2) = [ P_0 + \frac{2 \sigma (R_0)}{R_0} ](\frac{R}{R_0})^{-3\kappa} (1-\frac{3 \kappa}{c} \dot{R}) \\
- P_0 - \frac{2\sigma(R))}{R} - \frac{4\mu\dot{R}}{R} - \frac{4\kappa_s\dot{R}}{R^2} - P_{ac}(t)
\end{multline}

\begin{equation} \label{eqn:marmottant_sgm}
\sigma(R) =
\begin{cases}
0 & \text{if} R \leq R_{buckling} \\
\chi (\frac{R^2}{R_{buckling}^2} - 1) & \text{if} R_{buckling} \leq R \leq R_{break} \\
\sigma_{water} & \text{if ruptured and } R \geq R_{ruptured}       \\
\end{cases}
\end{equation}

This ode is solved using Runge-Kutta methods at a higher sampling rate than the ultrasound simulation, with an user-defined oversampling factor. Solving the ode yields the bubble radius over time and its first two derivatives: $R$, $\dot{R}$, and $\ddot{R}$. The scattered pressure at a distance $d$ from the center of the bubble is defined <cit.>:

\begin{equation} \label{eqn:bub_scatered}
P=\frac{\rho_l}{d}(R^2 \ddot{R} + 2R\dot{R}^2)
\end{equation}

bff has implementations of both a custom bubble that users can setup, and a commercial agent such as SonoVue® (sulphur hexafluoride) that is pre-configured with fitted parameters <cit.> as shown in Table <ref>).

Param     Description     SonoVue

$\rho_l$     Density of surrounding liquid     $10^3 \frac{Kg}{m^3}$

$\sigma_l$     Surface tension of surrounding liquid     $0.073 \frac{N}{m}$

$\mu_l$     Viscosity of surrounding liquid     $2.0 10^{-3} P s$

$\kappa$     Polytropic gas exponent     $1.095$

$\kappa_s$     Surface dilatational viscosity from the mono-layer     $ 7.2 10^{-9} N$

$\chi$     Elastic compression modulus of the mono-layer     $1.0 \frac{N}{m}$

$R_0$     Equilibrium radius of the bubble     $0.975 \mu m$

$R_{buckle}$     Lower radius limit of the elastic state     $0.975 \mu m$

Microbubble Simulation Parameters for a bubble of radius $3\mu m$

§.§ Integration into Linear Acoustic Simulation

Field II is a platform that simulates pressure fields from arbitrarily shaped transducers <cit.>.
It models the ultrasound system as a Linear Time Invariant system which is characterized by it's impulse response; the response to a Dirac excitation. This assumption of linearity makes Field II fast, allowing it to divide the problem into parallelizable parts, by splitting each transducer element into small mathematical sub-elements, and making every scatterer independent. Over the years, Field II's results have been shown to be accurate with ultrasound experiments <cit.>. The assumption of linearity is broken when incorporating the non linear dynamics of the contrast agent. To integrate bubble simulation into the FieldII's acoustic simulation, a multi step process is required. It is assumed that bubbles are independent and there is no coupling between them.

First, the transmitted pressure signal at the location of the bubble is calculated. This is done by simulating the transducer and calculating the pressure over time at the bubble location. Second, using the pressure signal calculated at step 1 as an input, we calculate the bubble response signal by using equation <ref> and then equation <ref>. Finally, the bubble response signal is propagated back to the transducer. Field II does not have a direct way to calculate the signal received with the transducer when transmitting from a point source in space. This can be done by convolving the bubble response and the spatial impulse response from the bubble position to the transducer, which, because of the reversibility property of ultrasound, is equal to the impulse response from the transducer to a mb position.

\begin{equation}
RX(t) = \text{bub{\_}resp}(t) \ast H(t, \text{bub{\_}pos})
\end{equation}

§.§ Evaluation

The evaluation metrics of bff present estimations of both the fractions of correct events and the error for localization and tracking of an ulm algorithm.

Localization Evaluation
A localization is considered a tp if it lies close enough to a ground truth position; the localisation is inside a defined search radius centered around the ground truth location. Following this logic, a ground truth mb with no localisation in its search radius is considered a fn, and a localization that is not inside the search radius of any ground truth fp location is a fp. Two metrics, precision and recall are used to indicate the fractions of the correct localisation:

\begin{equation} \label{eqn:loc_precision}
Precision = \frac{TP}{TP + FP}
\end{equation}

\begin{equation} \label{eqn:loc_recall}
Recall = \frac{TP}{TP + FN}
\end{equation}

A third metric which describes the average distance between the true position of a mb (B) and the localisation position (L) is calculated to reflect the root mean square error (RMSE) of all true positive localisation ($N_{TP}$):

\begin{equation} \label{eqn:track_mean_err}
RMSE = \frac{\sum_{i\in TP} dist(B_i \leftrightarrow L_i)}{N_{TP}}
\end{equation}

Tracking Evaluation

Tracking is a process that analyzes localizations over consecutive frames and determines pairs of localizations that correspond to the same bubble. Tracking is strongly dependent on the localization performance. A bubble that was not localized will never be paired, and a fp localization should never be paired. If one were to consider all localizations for tracking, fn and fp localizations would generate a double penalty on the erroneous localizations. To separate the tracking evaluation from the localization evaluation, only true positive localizations must be considered, i.e. only pairings of correctly localized bubbles in both frames are analyzed. This diminishes the total number of pairings to work with. As the total number of localisation increases the impact of this systematic bias is reduced.

Across two consecutive frames, two paired localizations that have the same localization id are considered tp if their associated ground truth bubbles have the same bubble id. Any bubble pair that doesn't have an associated localization pair is considered a fn. Any localization pair that is not associated with a bubble pair is considered a fp. Tracking is evaluated both in fraction of correct pairs and in effective pair distance. Precision and Recall are used to indicate the fractions of correct pairings:

\begin{equation} \label{eqn:track_precision}
Precision = \frac{TP}{TP + FP}
\end{equation}

\begin{equation} \label{eqn:track_recall}
Recall  = \frac{TP}{TP + FN}
\end{equation}

There is no direct way of measuring distance error for a pairing process.
The correct pair distance can be used, i.e. the distance a bubble travels from one frame to the other.
It is a fair metric, bubbles that travel further away are more difficult to track, which would yield a better result in the metric.
Given that the total number of pairs is dependent on the localization process, the fraction of correct paired distance. This is equivalent to a weighted Jaccard index with the weights being the distance of each pair.

\begin{equation} \label{eqn:track_jaccaard}
J = \frac{TP_d}{TP_d + FP_d + FN_d}
\end{equation}

\begin{equation}
TP_d = \sum_{i\in TP} dist(pair_i)
\end{equation}

\begin{equation}
FP_d = \sum_{i\in FP} dist(pair_i)
\end{equation}

\begin{equation}
FN_d = \sum_{i\in FN} dist(pair_i)
\end{equation}

To provide better clarity, a remapping of this metric of the form $L=2*J-1$ is used, so the range of the metric goes from $[0,1] \rightarrow [-1, 1]$. A negative number means track is mostly wrong, and positive number is mostly right. The final metric has then the following expression:

\begin{equation} \label{eqn:track_metric}
J_{map} = \frac{TP_d - FP_d - FN_d}{TP_d + FP_d + FN_d}
\end{equation}

§.§ BFF Validation Experiments

§.§.§ Network Generation

Multiple networks were generated with random parameters. In addition, a kidney and liver were mimicked to illustrate the ability of generating organic like structures with little effort. The network generator was manually constrained only by the $inside\_f$ to yield a similar geometric shape.

§.§.§ Bubble Simulation

A simulation was performed to validate bff's non linear simulation capabilities and replicated in vitro. The purpose was to test for psf realism and variety. The in-vitro experiment was done with an L11-5v 128-element linear array transducer with a tansmitted frequency of 7 MHz. Perfluorobutane mb where created following the process described in <cit.>, and placed in a beaker with de-gassed, filtered water at room temperature. Images were acquired with a Verasonics Vantage 256 platform (Verasonics Inc., Redmond, WA). Plane wave imaging compounding using three angles and a transmission amplitude corresponding to a mi of 0.05 was used. The simulation replicated the transducer, and mb with random parameters and locations. The resulting images where visually assessed and psf of similar morphology were identified and extracted.

§.§.§ ultra-sr Challenge

bff's qualities were demonstrated at the ultra-sr challenge at the ius 2022 of the ieee. Contrast enhanced ultrasound videos from four randomly generated networks with ground truth were produced for this competition. The proposed evaluation criteria were used for objective and quantitative evaluation of the participants localisation and tracking algorithms. The synthetic data was created using networks with varying density on different spatial sections. To minimize the predictability of the path of a mb the vessel tree generation parameters e.g., bifurcation probability, angle between consecutive segments or mb seeding probability were randomized. A training dataset was also provided with a single branch with a low mb concentration. The competition data combined multiple branching vessel trees into a single network with high mb concentration.

§ RESULTS

§.§ Network Generation

Figure <ref> (A) shows randomly generated networks with a different number of maximum allowed recursion levels. The vessels are color coded according to their size. Different amount of tortuosity and three-dimensionality are demonstrated with complexity of the network increasing from left to right. Networks are constrained e.g. by the $inside\_f$ function to a rectangular shape. Figure <ref> (B, C) shows the ability of bff to generate physiological shapes based on such simple constraints. bff is used to generate a network closely resembling the corony microvascularture (B) and the shape of a kidney (C). The parameterised generation of networks can generate any arbitrary shape.

A: Randomly generated networks obtained using different maximum allowed recursion levels. Recursion level from left to right: lvl=2, lvl=3, lvl=4, lvl=5. Rectangular constraint on network shape. B: Superresolution Ultrasound Image of a Rabbit Kidney<cit.>. C: Generated Network mimicking a the ramifications of one of the interlobar arteries of the rabbit's kidney.

§.§ Network Seeding

<ref> shows the seeding of a randomly generated network along the tracks with a bubble radial position r=0.

Particles seedded on the network

§.§ Bubble Simulation

bff's parameterised bubble generation results in a variety of shapes of the point spread function associated with a mb. Figure <ref> qualitatively shows the difference in the B-Mode image between simulated mb (B, D) and mb from an in vitro beaker acquisition (A, C). The mono-lobe (A, B) and multi-lobe (C, D) morphological shapes are distinct from each other. The side lobe signal is stronger in the simulation compared to the in vitro experiment.

Comparison between psf morphology seen in simulation and in vitro. A & C: Mono-lobe and multi-lobe morphologies observed at in vitro beaker experiments. B & D: Mono-lobe and multi-lobe morphologies can also be observed when recreating the experiment using bff.

§.§ ultra-sr Challenge

Two randomly generated networks of mb seeded vessel trees (<ref>) were simulated using bff as shown in Figure <ref>. The final networks were generated from multiple smaller networks to make localization and tracking difficult. No apparent structure within the networks is recognisable. <ref> shows the B-Mode images of both networks. It demonstrates the effect of adding Additive Colored Noise and Gaussian White Noise to the beamformed radio frequency signals. Through a Time Gain Compensation the signal to noise ratio decreases with depth. The signal to noise ratio is higher in the low frequency dataset (<ref>, B) compared to the higher frequency dataset (<ref>, A).

<ref> shows the B-Mode frames, superimposed with both the ground truth locations and the predicted locations. The two expanded regions show in detail how accurate the predictions are. The prescence of fn (isolated green cross) locations can be appreciated in both the hf and lf images. Additionally, both images (C & D) show two neighbouring bubbles being identified as one.

The final superresolution image and its velocity map are shown in <ref> and <ref>. While the observed network resembles the one presented in <ref>, the density of tracks is not the same, and large black gaps can be seen. This result is expected when the localization and tracking results (<ref>) are taken into consideration. For the hf dataset, a precision$=0.51$ indicates that around 50% of the true localizations was not captured, and of the possible pairings to be made with the correctly localized positions, only 43% were actually paired. On the other hand, the results also show a high tracking precision of $0.99$ for the hf dataset, which means that only 1% of the tracks present in the final image will be inexisting lines created by the algorithm.
The performance of the algorithm is worse for the lf dataset despite the higher signal to noise ratio. This was expected because of its much larger psf.

Structures generated with bff for the ultra-sr. A: The vessel network for the hf simulation. B: The vessel network for the lf simulation

Generated B-Mode frames. A: Frame N149 of the hf video simulation. B: Frame N4 of the lf video simulation

Superposition of a B-Mode frame and the bubble localizations. The ground truth is shown as a green cross, and the algorithm output is shown as a red cross. A & C: hf video frame N149 and detail of a region. B & D: lf video frame N4 and detail of a region.

Resulting SR Images obtained from localization and tracking. A: hf network. B: lf network

Resulting Velocity Map Image obtained from tracking. A: hf network. B: lf network

Precision     0.48

Recall     0.51

RMSE     $ 135.79 \mu m$

Precision     0.99

Recall     0.43

$J_{map}$     -0.71

Low Frequency

Precision     0.74

Recall     0.66

RMSE     $ 47.34 \mu m$

Precision     0.56

Recall     0.44

$J_{map}$     -0.57

High Frequency

Evaluation results

§ DISCUSSION

We developed bff, a fully comprehensive simulation platform for ulm algorithm development and evaluation. The framework consists of four key components: a microvascular structure generator, a flow simulator, an acoustic field simulator coupled with nonlinear microbubble dynamics and an evaluation pipeline for binary and quantitative assessment. The framework code is object oriented which makes it modular and easy to extend. bff allows the user to create short, concise and powerful scripts that simulate dynamic mb excited by common ultrasound transducers and is open source.

The microvessel structure generator can be customized to yield networks that resemble realistic vascular structures, <ref>.A shows how parameters like the recursion level, can generate models with increasing vessel densities. Organic structures like the vessels structures seen on kidneys, can be replicated. <ref>.B shows how the ramifications of interlobar arteries of the rabbit's kidney can mimicked. To this point, only binary trees <ref> are capable of being generated, but the implementation can be extended to generate other networks. The addition of loop connections could help mimic intussusceptive angiogenesis for even more realistic structures.

The incorporation of non-linear simulation of mb via Modified Rayleigh–Plesset ode has proven to generate realistic psf (<ref>) and enable simulation of coded transmission techniques such as am. The simulation not only is realistic in terms of the shapes and sizes of the resulting psf, but also realistic rf data is generated. This means that the datasets generated using bff can be used to assess algorithms and techniques that require that type of input, for example beamforming.

bff helps address one important problem present in ulm algorithm creation, evaluation and validation: the creation of realistic microvascular phantoms and data with ground truth. The realism of bff’s simulations greatly facilitate the evaluation of ulm algorithms, and the datasets generated by bff framework has been used by the ultra-sr challenge (<https://ultra-sr.com>) at the ieee ius 2022 to evaluate a number of localisation and tracking algorithms.

Even though this work focuses on ulm as a proof of concept, bff is not limited to this application only. Any contrast imaging related application can greatly benefit from datasets with ground truth for evaluation. Additionally, the accurate flow simulation is of special interest for applications like Ultrasound Vector Flow Imaging or Perfusion Imaging.

Data driven algorithms such as deep learning methods will benefit from this framework thanks to its compute time which is adequate for large dataset generation. The most time consuming step in the simulation framework is the acoustic field simulation. For the purpose of keeping it as open-source bff uses the free version of FieldII<cit.>, further speed up can be achieved by using the commercial version. It should be noted that Field II does not simulate non-linear propagation of ultrasound. However, given the low mechanical index used for mb imaging, the effect of nonlinear propagation, particularly when bubbles are excited at resonance frequency, is low <cit.>.

§ CONCLUSION

This work introduces bff, a simulation framework for algorithm development and evaluation, with a focus on fast and large dataset creation, which will simplify the development and assesment of deep learning models and traditional algorithms for all the stages in the ulm pipeline.
As a proof of concept, bff was used for the ULTRA-SR competition at ieee ius 2022, proving that it can create random organic networks that resemble physiological structures seen in-vivo, it can yield flow simulations that are accurate representation so of actual microvascular flow, it uses simulation of mb dynamics to generate realistic psf that resemble the ones seen on experiments, and it is suitable for large dataset generation thanks to its speed.

§ ACKNOWLEDGEMENTS

This work was supported by UK epsrc.

The complete simulated dataset for ultra-sr can be found on <https://doi.org/10.5281/zenodo.7271766>

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

abstractCardiac imaging has a pivotal role in the prevention, diagnosis and treatment of ischaemic heart disease. SPECT is most commonly used for clinical myocardial perfusion imaging, whereas PET is the clinical reference standard for the quantification of myocardial perfusion. MRI does not involve exposure to ionizing radiation, similar to echocardiography, which can be performed at the bedside. CT perfusion imaging is not frequently used but CT offers coronary angiography data, and invasive catheter-based methods can measure coronary flow and pressure. Technical improvements to the quantification of pathophysiological parameters of myocardial ischaemia can be achieved. Clinical consensus recommendations on the appropriateness of each technique were derived following a European quantitative cardiac imaging meeting and using a real-time Delphi process. SPECT using new detectors allows the quantification of myocardial blood flow and is now also suited to patients with a high BMI. PET is well suited to patients with multivessel disease to confirm or exclude balanced ischaemia. MRI allows the evaluation of patients with complex disease who would benefit from imaging of function and fibrosis in addition to perfusion. Echocardiography remains the preferred technique for assessing ischaemia in bedside situations, whereas CT has the greatest value for combined quantification of stenosis and characterization of atherosclerosis in relation to myocardial ischaemia. In patients with a high probability of needing invasive treatment, invasive coronary flow and pressure measurement is well suited to guide treatment decisions. In this Consensus Statement, we summarize the strengths and weaknesses as well as the future technological potential of each imaging modality.
journaltitleNature Reviews Cardiology
titleClinical quantitative cardiac imaging for the assessment of myocardial ischaemia
Computed tomography,Echocardiography,Ischaemia,Magnetic resonance imaging,Positron-emission tomography,Radionuclide imaging

abstractAims: Angiogenesis and vessel organisation in laryngeal tumour development and progression were examined to determine characteristics of biological and clinical relevance., Methods: Automated quantitative image analysis was performed on 1451 factor VIII (FVIII) associated blood vessels with regard to occurrence, structure, size, shape, and staining intensity, in addition to vessel direction., Results: Vessel numbers were increased in preneoplastic states and severe dysplasia, in addition to squamous cell carcinomas, being greater in poorly differentiated carcinomas. Small regular vessels predominated in benign conditions and large, irregular vessels in malignant neoplasms. Vessel distribution was related to degree of differentiation in squamous cell carcinomas, with circumferential angiogenesis occurring in well differentiated neoplasms, directional angiogenesis in moderately differentiated tumours, and aberrant angiogenesis in less well differentiated neoplasms. Alterations in vessel shape increased significantly with increasing degree of malignancy. Comparing the characteristics of individual vessels showed vessel shape abnormalities and the intensity of FVIII staining to increase with vessel size., Conclusions: Increased angiogenesis was an early event in laryngeal tumour development, with vessel structure, size, and shape related to the tumour growth pattern and behaviour.
journaltitleJournal of Clinical Pathology
titleSize, shape, structure, and direction of angiogenesis in laryngeal tumour development

titleNeovascularization in Human Atherosclerosis

abstractUltrasound super-resolution imaging through localisation and tracking of microbubbles can achieve sub-wave-diffraction resolution in mapping both micro-vascular structure and flow dynamics in deep tissue in vivo. Currently, it is still challenging to achieve high accuracy in localisation and tracking particularly with limited imaging frame rates and in the presence of high bubble concentrations. This study introduces microbubble image features into a Kalman tracking framework, and makes the framework compatible with sparsity-based deconvolution to address these key challenges. The performance of the method is evaluated on both simulations using individual bubble signals segmented from in vivo data and experiments on a mouse brain and a human lymph node. The simulation results show that the deconvolution not only significantly improves the accuracy of isolating overlapping bubbles, but also preserves some image features of the bubbles. The combination of such features with Kalman motion model can achieve a significant improvement in tracking precision at a low frame rate over that using the distance measure, while the improvement is not significant at the highest frame rate. The in vivo results show that the proposed framework generates SR images that are significantly different from the current methods with visual improvement, and is more robust to high bubble concentrations and low frame rates.
journaltitleIEEE transactions on medical imaging
titleSuper-Resolution Ultrasound Through Sparsity-Based Deconvolution and Multi-Feature Tracking
Animals,Computer Simulation,Humans,Image Processing,Computer-Assisted,Mice,Microbubbles,Microvessels,Ultrasonography

abstractIdentifying and visualizing vasculature within organs and tumors has major implications in managing cardiovascular diseases and cancer. Contrast-enhanced ultrasound scans detect slow-flowing blood, facilitating noninvasive perfusion measurements. However, their limited spatial resolution prevents the depiction of microvascular structures. Recently, super-localization ultrasonography techniques have surpassed this limit. However, they require long acquisition times of several minutes, preventing the detection of hemodynamic changes. We present a fast super-resolution method that exploits sparsity in the underlying vasculature and statistical independence within the measured signals. Similar to super-localization techniques, this approach improves the spatial resolution by up to an order of magnitude compared to standard scans. Unlike super-localization methods, it requires acquisition times of only tens of milliseconds. We demonstrate a temporal resolution of 25 Hz, which may enable functional super-resolution imaging deep within the tissue, surpassing the temporal resolution limitations of current super-resolution methods, e.g., in neural imaging. The subsecond acquisitions make our approach robust to motion artifacts, simplifying in vivo use of super-resolution ultrasound.
journaltitleIEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control
noteConference Name: IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control
titleSUSHI: Sparsity-Based Ultrasound Super-Resolution Hemodynamic Imaging
Blood,Compressed sensing (CS),Hemodynamics,Imaging,Signal resolution,Spatial resolution,Ultrasonic imaging,contrast-enhanced ultrasound (CEUS),sparse representation,super-localization microscopy,super-resolution,super-resolution optical fluctuation imaging

abstractSuper-resolution ultrasound localization microscopy (ULM), based on localization and tracking of individual microbubbles (MBs), offers unprecedented microvascular imaging resolution at clinically relevant penetration depths. However, ULM is currently limited by the requirement of dilute MB concentrations to ensure spatially sparse MB events for accurate localization and tracking. The corresponding long imaging acquisition times (tens of seconds or several minutes) to accumulate sufficient isolated MB events for full reconstruction of microvasculature preclude the clinical translation of the technique. To break this fundamental tradeoff between acquisition time and MB concentration, in this paper we propose to separate spatially overlapping MB events into sub-populations, each with sparser MB concentration, based on spatiotemporal differences in the flow dynamics (flow speeds and directions). MB localization and tracking are performed for each sub-population separately, permitting more robust ULM imaging of high-concentration MB injections. The superiority of the proposed MB separation technique over conventional ULM processing is demonstrated in flow channel phantom data, and in the chorioallantoic membrane of chicken embryos with optical imaging as an in vivo reference standard. Substantial improvement of ULM is further demonstrated on a chicken embryo tumor xenograft model and a chicken brain, showing both morphological and functional microvasculature details at super-resolution within a short acquisition time (several seconds). The proposed technique allows more robust MB localization and tracking at relatively high MB concentrations, alleviating the need for dilute MB injections, and thereby shortening the acquisition time of ULM imaging and showing great potential for clinical translation.
journaltitleScientific Reports
titleShort Acquisition Time Super-Resolution Ultrasound Microvessel Imaging via Microbubble Separation
Cancer imaging,Electrical and electronic engineering,Preclinical research

abstractPerfusion by the microcirculation is key to the development, maintenance and pathology of tissue. Its measurement with high spatiotemporal resolution is consequently valuable but remains a challenge in deep tissue. Ultrasound Localization Microscopy (ULM) provides very high spatiotemporal resolution but the use of microbubbles requires low contrast agent concentrations, a long acquisition time, and gives little control over the spatial and temporal distribution of the bubbles. The present study is the first to demonstrate Acoustic Wave Sparsely-Activated Localization Microscopy (AWSALM) and fast-AWSALM for in vivo super-resolution ultrasound imaging, offering contrast on demand and vascular selectivity. Three different formulations of sono-switchable contrast agents were tested. We demonstrate their use with ultrasound mechanical indices well within recommended safety limits to enable fast on-demand sparse switching at very high agent concentrations. We produce super-localization maps of the rabbit renal vasculature with acquisition times between 5.5 s and 0.25 s, and an 4-fold improvement in spatial resolution. We present the unique selectivity of AWSALM in visualizing specific vascular branches and downstream microvasculature, and we show super-localized kidney structures in systole and diastole with fast-AWSALM. In conclusion we demonstrate the feasibility of fast and selective measurement of microvascular dynamics in vivo with subwavelength resolution using ultrasound and sono-switchable nanodroplets.
notearXiv:2203.04263 [eess]
titleFast and selective super-resolution ultrasound in vivo with sono-switchable nanodroplets
Electrical Engineering and Systems Science - Signal Processing

abstractUltrasound localization microscopy has enabled super-resolution vascular imaging through precise localization of individual ultrasound contrast agents (microbubbles) across numerous imaging frames. However, analysis of high-density regions with significant overlaps among the microbubble point spread responses yields high localization errors, constraining the technique to low-concentration conditions. As such, long acquisition times are required to sufficiently cover the vascular bed. In this work, we present a fast and precise method for obtaining super-resolution vascular images from high-density contrast-enhanced ultrasound imaging data. This method, which we term Deep Ultrasound Localization Microscopy (Deep-ULM), exploits modern deep learning strategies and employs a convolutional neural network to perform localization microscopy in dense scenarios, learning the nonlinear image-domain implications of overlapping RF signals originating from such sets of closely spaced microbubbles. Deep-ULM is trained effectively using realistic on-line synthesized data, enabling robust inference in-vivo under a wide variety of imaging conditions. We show that deep learning attains super-resolution with challenging contrast-agent densities, both in-silico as well as in-vivo. Deep-ULM is suitable for real-time applications, resolving about 70 high-resolution patches ( 128×128 pixels) per second on a standard PC. Exploiting GPU computation, this number increases to 1250 patches per second.
journaltitleIEEE transactions on medical imaging
titleSuper-Resolution Ultrasound Localization Microscopy Through Deep Learning
Contrast Media,Deep Learning,Microbubbles,Microscopy,Ultrasonography

abstractA program for the simulation of ultrasound systems is presented. It is based on the Tupholme-Stepanishen method, and is fast because of the use of a far-field approximation. Any kind of transducer geometry and excitation can be simulated, and both pulse-echo and continuous wave fields can-be calculated for both transmit and pulse-echo. Dynamic apodization and focusing are handled through time lines, and different focusing schemes can be simulated. The versatility of the program is ensured by interfacing it to Matlab. All routines are called directly from Matlab, and all Matlab features can be used. This makes it possible to simulate all types of ultrasound imaging systems.
journaltitleMedical and Biological Engineering and Computing
titleFIELD: A program for simulating ultrasound systems

abstractA new, freely available third party MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields is described. The toolbox, named k-Wave, is designed to make realistic photoacoustic modeling simple and fast. The forward simulations are based on a k-space pseudo-spectral time domain solution to coupled first-order acoustic equations for homogeneous or heterogeneous media in one, two, and three dimensions. The simulation functions can additionally be used as a flexible time reversal image reconstruction algorithm for an arbitrarily shaped measurement surface. A one-step image reconstruction algorithm for a planar detector geometry based on the fast Fourier transform (FFT) is also included. The architecture and use of the toolbox are described, and several novel modeling examples are given. First, the use of data interpolation is shown to considerably improve time reversal reconstructions when the measurement surface has only a sparse array of detector points. Second, by comparison with one-step, FFT-based reconstruction, time reversal is shown to be sufficiently general that it can also be used for finite-sized planar measurement surfaces. Last, the optimization of computational speed is demonstrated through parallel execution using a graphics processing unit.
journaltitleJournal of Biomedical Optics
titlek-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields
Algorithms,Elasticity Imaging Techniques,Electromagnetic Fields,Image Enhancement,Image Interpretation,Computer-Assisted,Imaging,Three-Dimensional,Programming Languages,Reproducibility of Results,Sensitivity and Specificity,Signal Processing,Computer-Assisted,Software

titleSimcenter STAR-CCM+

abstractBy localizing microbubbles (MBs) in the vasculature, ultrasound localization microscopy (ULM) has recently been proposed, which greatly improves the spatial resolution of ultrasound (US) imaging and will be helpful for clinical diagnosis. Nevertheless, several challenges remain in fast ULM imaging. The main problems are that current localization methods used to implement fast ULM imaging, e.g., a previously reported localization method based on sparse recovery (CS-ULM), suffer from long data-processing time and exhaustive parameter tuning (optimization). To address these problems, in this paper, we propose a ULM method based on deep learning, which is achieved by using a modified sub-pixel convolutional neural network (CNN), termed as mSPCN-ULM. Simulations and in vivo experiments are performed to evaluate the performance of mSPCN-ULM. Simulation results show that even if under high-density condition (6.4 MBs/mm2), a high localization precision ( 28 μm in the lateral direction and 24 μm in the axial direction) and a high localization reliability (Jaccard index of 0.66) can be obtained by mSPCN-ULM, compared to CS-ULM. The in vivo experimental results indicate that with plane wave scan at a transmit center frequency of 15.625 MHz, microvessels with diameters of 17 μm can be detected and adjacent microvessels with a distance of 42 μm can be separated. Furthermore, when using GPU acceleration, the data-processing time of mSPCN-ULM can be shortened to 6 sec/frame in the simulations and 23 sec/frame in the in vivo experiments, which is 3-4 orders of magnitude faster than CS-ULM. Finally, once the network is trained, mSPCN-ULM does not need parameter tuning to implement ULM. As a result, mSPCN-ULM opens the door to implement ULM with fast data-processing speed, high imaging accuracy, short data-acquisition time, and high flexibility (robustness to parameters) characteristics.
journaltitleIEEE Transactions on Medical Imaging
noteConference Name: IEEE Transactions on Medical Imaging
titleDeep Learning for Ultrasound Localization Microscopy
Convolution,Convolutional neural networks,In vivo,Microscopy,Spatial resolution,Ultrasonic imaging,deep learning,sub-pixel convolution layer,super-resolution ultrasound imaging,ultrasound localization microscopy

abstractBlood flow velocity and wall shear stress (WSS) influence and are influenced by vascular disease. Their measurement is consequently useful in the laboratory and clinic. Contrast-enhanced ultrasound image velocimetry (UIV) can estimate them accurately but the need to inject contrast agents limits utility. Singular value decomposition and high-frame-rate imaging may render contrast agents dispensable. Here we determined whether contrast agent-free UIV can measure flow and WSS. In simulation, accurate measurements were achieved with a signal-to-noise ratio of 13.5 dB or higher. Signal intensity in the rabbit aorta increased monotonically with mechanical index; it was lowest during stagnant flow and uneven across the vessel. In vivo measurements with contrast-free and contrast-enhanced UIV differed by 4.4% and 1.9% for velocity magnitude and angle and by 9.47% for WSS. Bland–Altman analysis of waveforms revealed good agreement between contrast-free and contrast-enhanced UIV. In five rabbits, the root-mean-square errors were as low as 0.022 m/s (0.81%) and 0.11 Pa (1.7%). This study indicates that with an optimised protocol, UIV can assess flow and WSS without contrast agents. Unlike contrast-enhanced UIV, contrast-free UIV could be routinely employed.
journaltitleUltrasound in Medicine & Biology
titleContrast Agent-Free Assessment of Blood Flow and Wall Shear Stress in the Rabbit Aorta using Ultrasound Image Velocimetry
Atherosclerosis,Blood flow,Blood speckle velocimetry,Echo-particle image velocimetry,Native blood speckle,Ultrasound image velocimetry,Vector flow imaging,Wall shear stress

abstractUltrasound is extensively used in medical imaging, being safe and inexpensive and operating in real time. Its scope of applications has been widely broadened by the use of ultrasound contrast agents (UCAs) in the form of microscopic bubbles coated by a biocompatible shell. Their increased use has motivated a large amount of research to understand and characterize their physical properties as well as their interaction with the ultrasound field and their surrounding environment. Here we review the theoretical models that have been proposed to study and predict the behavior of UCAs. We begin with a brief introduction on the development of UCAs. We then present the basics of free-gas-bubble dynamics upon which UCA modeling is based. We review extensively the linear and non-linear models for shell elasticity and viscosity and present models for non-spherical and asymmetric bubble oscillations, especially in the presence of surrounding walls or tissue. Then, higher-order effects such as microstreaming, shedding and acoustic radiation forces are considered. We conclude this review with promising directions for the modeling and development of novel agents.
journaltitleUltrasound in Medicine & Biology
shorttitleUltrasound Contrast Agent Modeling
titleUltrasound Contrast Agent Modeling: A Review
Microbubble dynamics,Modeling,Ultrasound contrast agents

abstractWe present a model applicable to ultrasound contrast agent bubbles that takes into account the physical properties of a lipid monolayer coating on a gas microbubble. Three parameters describe the properties of the shell: a buckling radius, the compressibility of the shell, and a break-up shell tension. The model presents an original non-linear behavior at large amplitude oscillations, termed compression-only, induced by the buckling of the lipid monolayer. This prediction is validated by experimental recordings with the high-speed camera Brandaris 128, operated at several millions of frames per second. The effect of aging, or the resultant of repeated acoustic pressure pulses on bubbles, is predicted by the model. It corrects a flaw in the shell elasticity term previously used in the dynamical equation for coated bubbles. The break-up is modeled by a critical shell tension above which gas is directly exposed to water.
journaltitleJournal of The Acoustical Society of America - J ACOUST SOC AMER
titleA model for large amplitude oscillations of coated bubbles accounting for buckling and rupture

abstractUltrasound super-resolution techniques use the response of microbubble (MB) contrast agents to visualize the microvasculature. Techniques that localize isolated bubble signals first require detection algorithms to separate the MB and tissue responses. This work explores the three main MB detection techniques for super-resolution of microvasculature. Pulse inversion (PI), differential imaging (DI), and singular value decomposition (SVD) filtering were compared in terms of the localization accuracy, precision, and contrast-to-tissue ratio. MB responses were simulated based on the properties of Sonovue and using the Marmottant model. Nonlinear propagation through tissue was modeled using the k-Wave software package. For the parameters studied, the results show that PI is most appropriate for low frequency applications, but also most dependent on transducer bandwidth. SVD is preferable for high frequency acquisition where localization precision on the order of a few microns is possible. PI is largely independent of flow direction and speed compared to SVD and DI, so is appropriate for visualizing the slowest flows and tortuous vasculature. SVD is unsuitable for stationary MBs and can introduce a localization error on the order of hundreds of microns over the speed range 0-2 mm/s and flow directions from lateral (parallel to probe) to axial (perpendicular to probe). DI is only suitable for flow rates >0.5 mm/s or as flow becomes more axial. Overall, this study develops an MB and tissue nonlinear simulation platform to improve understanding of how different MB detection techniques can impact the super-resolution process and explores some of the factors influencing the suitability of each.
journaltitleIEEE transactions on ultrasonics, ferroelectrics, and frequency control
titleInvestigation of Microbubble Detection Methods for Super-Resolution Imaging of Microvasculature
Algorithms,Computer Simulation,Image Processing,Computer-Assisted,Microbubbles,Microvessels,Models,Biological,Phantoms,Imaging,Transducers,Ultrasonography

abstractFree nonlinear oscillations of a gas bubble in a noncompressible liquid are analysed. Excitation to the growth is considered in the form of an instantaneous delivery of heat into the bubble interior. The characteristic features of the pressure and velocity fields in the liquid are discussed.
journaltitleCzechoslovak Journal of Physics B
titleOn Rayleigh's model of a freely oscillating bubble. I. Basic relations
Basic Relation,Bubble Interior,Characteristic Feature,Nonlinear Oscillation,Velocity Field

abstractA method for simulation of pulsed pressure fields from arbitrarily shaped, apodized and excited ultrasound transducers is suggested. It relies on the Tupholme-Stepanishen method for calculating pulsed pressure fields, and can also handle the continuous wave and pulse-echo case. The field is calculated by dividing the surface into small rectangles and then Summing their response. A fast calculation is obtained by using the far-field approximation. Examples of the accuracy of the approach and actual calculation times are given.<>
journaltitleIEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control
noteConference Name: IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control
titleCalculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers
Biomedical imaging,Biomedical transducers,Closed-form solution,Geometry,Partial differential equations,Pistons,Pulse shaping methods,Ultrasonic imaging,Ultrasonic transducers

abstractRow–column arrays have been shown to be able to generate 3-D ultrafast ultrasound images with an order of magnitude less independent electronic channels than traditional 2-D matrix arrays. Unfortunately, row–column array images suffer from major imaging artifacts due to high sidelobes, particularly when operating at high frame rates. This article proposes a row–column-specific beamforming technique, for orthogonal plane-wave transmissions, row–column-specific frame multiply and sum (RC-FMAS), that exploits the incoherent nature of certain row–column array artifacts. A series of volumetric images is produced using row or column transmissions of 3-D plane waves. The voxelwise geometric mean of the beamformed volumetric images from each row and column pair is taken prior to compounding, which drastically reduces the incoherent imaging artifacts in the resulting image compared to traditional coherent compounding. The effectiveness of this technique was demonstrated in silico and in vitro, and the results show a significant reduction in sidelobe level with over 16-dB improvement in sidelobe to main-lobe energy ratio. Significantly improved contrast was demonstrated with contrast ratio increased by 10 dB and generalized contrast-to-noise ratio increased by 158% when using the proposed new method compared to the existing delay and sum during in vitro studies. The new technique allowed for higher quality 3-D imaging while maintaining high frame rate potential.
journaltitleIEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control
noteConference Name: IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control
titleUltrafast 3-D Ultrasound Imaging Using Row–Column Array-Specific Frame-Multiply-and-Sum Beamforming
3-D ultrasound,Array signal processing,Coherence,Delays,Imaging,In vitro,Probes,Ultrasonic imaging,beamforming,coherence,frame multiply and sum (FMAS),high frame rate/ultrafast imaging,real time,row and column array (RCA)

abstractThis paper investigates two types of nonlinear propagation and their effects on image intensity and contrast-to-tissue ratio (CTR) in contrast ultrasound images. Previous studies have shown that nonlinear propagation can occur when ultrasound travels through tissue and microbubble clouds, making tissue farther down the acoustic path appear brighter in pulse inversion (PI) images, thus reducing CTR. In this study, the effect of nonlinear propagation through tissue or microbubbles on PI image intensity and CTR are compared at low mechanical index. A combination of simulation and experiment with SonoVue microbubbles were performed using a microbubble dynamics model, a laboratory ultrasound system and a clinical prototype scanner. The results show that, close to the bubble resonance frequency, nonlinear propagation through a bubble cloud of a few centimeter thickness with a modest concentration (1:10000 dilution of SonoVue microbubbles) is much more significant than through tissue-mimicking material. Consequently, CTR in regions distal to the imaging probe is greatly reduced for nonlinear propagation through the bubble cloud, with as much as a 12-dB reduction compared with nonlinear propagation through tissue-mimicking material. Both types of nonlinear propagation cause only a small change in bubble PI signals at the bubble resonance frequency. When the driving frequency increases beyond bubble resonance, nonlinear propagation through bubbles is greatly reduced in absolute values. However because of a greater reduction in nonlinear scattering from bubbles at higher frequencies, the corresponding CTR is much lower than that at bubble resonance frequency. (E-mail<EMAIL_ADDRESS>journaltitleUltrasound in Medicine & Biology
titleEffects of Nonlinear Propagation in Ultrasound Contrast Agent Imaging
Contrast agents,Imaging artefacts,Nonlinear propagation,Perfusion quantification