Multi-Plane Ultrafast Compound 3 D Strain Imaging : Experimental Validation in a Carotid Bifurcation Phantom

Strain imaging of the carotid artery (CA) has demonstrated to be a technique capable of identifying plaque composition. This study assesses the performance of volumetric strain imaging derived from multi-plane acquisitions with a single transducer, with and without displacement compounding. These methods were compared to a reference method using two orthogonally placed transducers. A polyvinyl alcohol phantom was created resembling a stenotic CA bifurcation. A realistic pulsatile flow was imposed on the phantom, resulting in fluid pressures inducing 10% strains. Two orthogonally aligned linear array transducers were connected to two Verasonics systems and fixed in a translation stage. For 120 equally spaced elevational positions, ultrasound series were acquired for a complete cardiac cycle and synchronized using a trigger. Each series consisted of ultrafast plane-wave acquisitions at 3 alternating angles. Inter-frame displacements were estimated using a 3D cross-correlation-based tracking algorithm. Horizontal displacements were acquired using the single probe lateral displacement estimate, the single probe compounded by axial displacement estimates obtained at angles of 19.47 and −19.47 degrees, and the dual probe registered axial displacement estimate. After 3D tracking, least squares strain estimations were performed to compare compressive and tensile principal strains in 3D for all methods. The compounding technique clearly outperformed the zero-degree method for the complete cardiac cycle and resulted in more accurate 3D strain estimates.


Introduction
Ultrasound (US) strain imaging [1] is used for a wide variety of biomedical applications.One of them is functional assessment of the carotid vessel wall.Strain imaging holds the ability to quantify the mechanical deformation of the arterial wall induced by the pulsating blood and provides insight into the composition and morphology of the vessel wall, which is unavailable from ultrasonic imaging alone [2][3][4][5].Therefore, strain imaging could reveal ongoing, frequently asymptomatic, pathological processes like plaque development.This is very important, because eventually plaque could develop into either stable (predominantly thick fibrous tissue) or vulnerable plaque (large thrombogenic lipd core sealed with a thin cap).The latter type is prone to rupture [6,7].Since rupture of these plaques has been shown to be one of the main initiators of acute transient ischemic attack or stroke [8], it implies that strain imaging might enable up-front risk stratification of these potentially lethal events.With that it will also provide opportunities for early intervention and hopefully aid in the prevention of these events.
Many arterial strain imaging techniques have already been developed and reported, like intravascular techniques [9], which use an US transducer on the tip of a catheter to image from within the lumen, and non-invasive techniques, which image from outside the body [10][11][12][13][14][15].With these techniques, it has been shown that arterial strain distribution correlates with histological composition of plaques in pigs [16] and in humans [2,17].For intravascular imaging and non-invasive imaging in longitudinal planes, the US beam is aligned with the principal direction of deformation of the artery (radial for a concentric homogeneous artery), which implies a direct and accurate estimation of the strains in the principal direction using algorithms that only estimate the (axial) displacement component along the beam.For non-invasive strain estimation in transverse planes, accurate estimation of the full 2D displacement vector requires more advanced methods [18][19][20][21].
Until now, all in vivo applied techniques have been primarily based on cross-sectional images of the artery and provide a strain estimation in 2D only.However, as illustrated before, 2D estimation can be problematic in the case of complex motion patterns as a result of longitudinal artery wall displacement [22,23].It has been demonstrated that strains are generally smaller in the longitudinal direction than in the radial direction, although the displacement magnitude in the longitudinal direction of the arterial wall equals the magnitude of the radial displacements [24].Reduced longitudinal artery wall motion was even correlated with cardiovascular outcome [25], which makes an estimate of the full 3D tissue motion in horizontal, vertical, and longitudinal direction even more relevant.
Full 3D tissue motion estimation overcomes the following limitations of 2D.First, due to longitudinal wall motion, 2D transverse US acquisitions will suffer from out of plane motion reducing the strain estimation accuracy.A second limitation of 2D scanning is that it is more difficult to find the optimal cross-section needed for appropriate assessment of carotid disease, i.e., the most vulnerable spot might be missed.Also, inter-and intra-operator variability in skill and experience could hamper a reliable judgement about the best response to the therapy of the plaques [26].
Three-dimensional B-mode-based US of the carotid artery [27,28] has been described before and has been shown to improve the visualization of the intima-media thickness.It complements Intima Media Thhickness (IMT), which is a typical 2D US-based biomarker [29] with a measurement of the Total Plaque Volume (TPV) [30] and the Vessel Wall Volume (VWV) [31].However, the implementation of 3D ultrasound in these studies did not allow measurement of 3D tissue motion or strain.The main aim of the present study is to develop a 3D ultrasound technique that enables 3D vascular strain estimation of the carotid artery.As aforementioned, adding this functionality could possibly lead to a significant improvement of the prediction of plaque vulnerability.
A few studies have investigated 3D strain imaging for vascular applications.In vitro studies have shown a 3D representation of the 2D strain tensor derived from multi-slice 2D displacement estimations [32,33].In this approach, the dominant tissue motion in radial direction was aligned with the US beam direction due to the rotational acquisition protocol using the longitudinal axis for rotation.However, the rotation applied in this approach cannot be performed noninvasively, and therefore translation towards clinical practice is not realistic.Furthermore, in reality, the vessel wall tends to deform in three orthogonal directions due to the arterial pressure pulse, and the vessel wall dynamics are affected by the heterogeneous composition and complex geometry in the presence of a stenosis or a plaque [24].Volumetric 2D strain tensor estimation based on 3D motion estimation showed the capability of handling out of plane motion, i.e., the direction perpendicular to the transversal image plane [23].However, measurement of the complete 3D strain tensor for the entire vessel wall of the carotid artery in 3D is desirable in the pursuit of better understanding the underlying true deformation encompassing the horizontal, vertical, and longitudinal directions.This requires 3D tracking of the 3D tissue deformation and the derivation of the full 3D strain tensor.
To our knowledge, only one study reporting 3D strain tensors for arteries using intravascular B-mode US was presented in cylindrical coordinates [34], resulting in radial, circumferential, and longitudinal strain distributions.However, predefining the coordinate space to be cylindrical assumes an equal cylindrical geometry of the vessel wall to relate to the dominant strain direction.In the case of a geometry such as the bifurcation, this relation does not hold.Principal strains based on the principles of continuum mechanics have the intrinsic characteristic of coordinate system independence and are therefore more suitable to report the full 3D strain distribution at the bifurcation.2D principal strains were reported for vascular applications based on 2D strain tensor in previous studies [35,36].In this paper we implemented the acquisition method of our previous simulation study in an experimental setup and extended the strain estimation algorithm to full 3D.To evaluate the performance of this implementation experimentally, we constructed a phantom of the same patient-specific geometry as used in the simulation study of PolyVinyl Alcohol (PVA) solutions, subjected it to a pulsating physiological flow, and evaluated the accuracy of our 3D strain estimation.The main benefit of the simulation study was the availability of a ground truth to compare the performance of different strategies.However, this study did not incorporate all ultrasound artifacts that are present when using real US acquisitions, like reverberations, shadowing, or inhomogeneity of the propagation speed of sound.To overcome these limitations and approximate in vivo situations more closely, the Finite Element Model (FEM) geometry of the simulation study was used to construct a PVA-based phantom using 3D printed molds.This yielded a realistic phantom of the bifurcation, delivering realistic mechanical behavior when subjected to a fluid pressure pulse.However, due to the re-production variability and stability of PVA-cryogel-based phantoms in terms of elastic modulus [37], no real ground truth can be derived from it.Therefore, a multi-perspective dual probe setup was used to derive optimal strain estimates that operate as the ground truth of this work.
Although it has been previously reported that 2D strain estimation based on displacement compounding outperforms zero-degree based 2D strain estimation [38], to our knowledge, the comparison has never been made for a complex geometry using a realistic phantom.

Experimental Setup
A patient-specific vessel wall geometry of a Carotid Artery (CA), including the common carotid, the internal and external carotid, and the bifurcation [23], was used to build an ultrasound-compatible (i.e., matching the acoustical properties of soft tissue) phantom made of PVA-cryogel (Mw 85.0-124.099+% hydrolyzed, Sigma-Aldrich, Zwijndrecht, The Netherlands).Two molds were created using CAD (Autodesk Inventor 2015, Student version, San Rafael, CA, USA) to convert the CA-geometry (Figure 1a) and the surrounding body into 3D printable casts.Subsequently, these molds were filled with a pre-heated ~85 • homogeneous mixture of 54% by weight (wt.%) distilled water, 36 wt.% anti-freezing agent (ethylene glycol, Sigma-Aldrich, Zwijndrecht, The Netherlands), and 10 wt.% PVA including either 2.0 wt.% or 0.5 wt.% Silica gel particles (Merck Kieselgel 60, 0.063-0.100mm, Boom B.V., Meppel, The Netherlands) for the CA-geometry and the surrounding tissue, respectively.After injection and degassing for 1 h, the solutions in the molds were subjected to 4 cycles of freezing (−25 • C) and thawing (21 • C) of 16 and 8 h, respectively.Prior to the last cycle, the surrounding tissue was casted around the pre-fixated carotid artery bifurcation, resulting in the final composed phantom (Figure 1b).This yielded a higher elastic modulus for the CA-geometry compared to the surrounding tissue.The Young's modulus of the CA embedded in surrounding tissue was experimentally estimated to be 225 kPa, which is within range of values corresponding to intima tissue from atherosclerotic plaques [39].The phantom was connected to an inlet flow tube at the common carotid artery side and to two tubes (internal and external carotid artery) at the outlet side.The connections were supported by a mounting bracket also allowing pressure measurements (TruWave pressure transducer, Edwards Lifesciences Corp., Irvine, CA, USA), see Figure 2b.Fluid was circulated in a closed loop circuit by a gear pump driven by in-house software written in LabView (National Instruments, Austin, TX, USA).The CA was pressurized using a realistic pulsatile flow waveform at 60 bpm with a mean flow of 0.62 L min −1 and a peak flow of 1.39 L/min.At both outlet tubes, adjustable flow resistors were set such that the dynamic pressure, measured at the inlet of the CA, was 150 over 100 mmHg, mimicking the pressures encountered in a stenotic carotid artery.With these pressure profiles, strains of up to 10% over the pressure cycle were induced.As can be observed in Figure 2a, the surrounding body contained two orthogonal surfaces for ultrasound imaging, from which the averaged depth to the bifurcation was ~20 mm.Two L12-5 ATL linear array probes were positioned in a probe holder, which allowed fine-tuning to ensure orthogonal fixation and coinciding transversal fields of view.Both probes were closely hovering above the phantom imaging surfaces, while ultrasound gel established the acoustic coupling and simultaneously prevented physical deformation of the phantom due to direct phantom-probe contact.The probe holder was connected to a translation stage enabling elevational movement for 120 equally spaced (0.1 mm) positions while maintaining a constant in plane orientation of the probes.Each transducer was connected to a Verasonics ultrasound system (a Vantage 256 and a V1, Verasonics, Kirkland, WA, USA).

Ultrasound Data Acquisition
Both Verasonics systems and the translation stage were triggered by a common trigger generated by the Labview program at the start of the flow waveform (i.e., the start of the systolic phase).This enabled the time synchronization of the ultrasound data acquired for one cardiac cycle.
At each elevational position ultrasound element data series were acquired at a repetition frequency of 50 Hz for 1 s, i.e., one pressure cycle.Figure 3   The steered transmissions were acquired by applying a linear time delay in transmission.The Pulse Repetition Frequency (PRF) of these series was set to 10 kHz, and the tissue was assumed to be in a quasi-static state of deformation, i.e., intra-acquisition series displacements were assumed to be zero, which seems reasonable given the high PRF and the relatively slow motion of the arterial wall.After completing all acquisitions series for the full cardiac cycle, the probes were translated to the next elevational position and the whole sequence was repeated.In this way, a complete volume over the full cardiac cycle was obtained.Finally, frame-skipping was applied at the diastolic phase, yielding an effective PRF of 25 Hz to increase strain estimation signal to noise ratio [40] and to reduce computation time for this phase of the pressure cycle.

Post Processing RF Element Data
Plane wave imaging allows for the conversion of unfocused element data into beamformed RF data at any given location within the region insonified.In this way, not only the axial and lateral grid ratio can be tailored to optimize the displacement estimation accuracy [41], but also the spacing of the beamforming grid can be adjusted to facilitate a fast combination of the steered acquisitions used in the process of displacement compounding as described before in [23].This beam-forming approach was also used in this study for a steering angle α equal to 19.47 • .An additional orthogonally steered grid (90 • ) was added to enable a match with the 0 • grid of the other probe.For all four beamforming grid orientations, a single displacement grid was defined.Figure 4 shows the grid layout for both probes, indicating the displacement grid defined by the coinciding points of all grid definitions.The spacing of the beamforming grids and the displacement grid (coinciding points) is listed in Table 1.Element data were beamformed using Delay And Sum (DAS) beamforming assuming a speed of sound of 1540 m/s, a dynamic focus F = 0.875, and Hamming apodization window in receive.Please note that the resolution is specified in terms of dx, dy, and dz in Table 1, in which x is oriented in the lateral direction, y is oriented in the axial direction, and z is oriented in the elevational direction.B-mode representations of the beamformed RF-data of the external and internal carotid artery are shown in Figure 4 for both probes and for steered (α = ±19.47• ) and non-steered (α = 0.0 • ) transmissions.The orthogonal positions of both probes relative to each other allows one to utilize the axial zero-degree estimated displacements of probe A as lateral displacement estimates for probe B and vice versa as depicted in the last column of Figure 4.In this way, highly accurate estimates of the lateral displacments can be obtained, because native RF-phase information can be used in the estimation process, which is normally absent when using a single probe.Finally, the 2D beamformed RF data were stacked according to their elevational position to create a volume matrix of RF data over time.Please regard video S1 (Multi-perspective_3D_Bmode.mp4) for a dynamic B-mode presentation of the carotid artery phantom for both probes.

3D Interframe Displacement Estimation
Displacement estimations were performed for the RF data beamformed on each of the four grids for each elevational position and time point for both probes.The resulting displacement estimates were labeled using an A or a B and a subscript.The A and B refer to the probe being used.The subscript refers to the acquisition angle and beamforming grid being used.Thus, A −α , refers to the displacements estimated using probe A and a beamforming grid and acquisition at steering angle −α.
Displacement estimation was performed using a 2-step, normalized Cross-Correlation (CC) based method using 3D kernels defined in Table 2.In the first step, envelope (demodulated RF) data were used to find a global displacement estimate at sample accuracy (i.e., full samples).Subsequently, subsample displacements were determined by 3D spline interpolation of the auto and cross-correlation function based on RF [42].Median filtering was applied for all three directions using 3D kernels (see Table 2).For the elevational direction, additional regularization was applied by allowing only sub-plane displacement estimations.Since probe A and B are positioned orthogonally, and the axial and lateral displacement directions of probe A are equal to the lateral and axial displacement directions of probe B, respectively, and vice versa.The orientation of Probe A is chosen to be the primary orientation throughout the rest of this manuscript.The axial and lateral direction of probe A equals the vertical and horizontal direction of the accumulated displacement estimates as depicted in Figure 5.The displacement estimates from both probes for different directions are combined in three different ways used to derive strain estimates.These methods are described below and named: single-probe zero-degree, single-probe compounding, and Dual-probe reference.

Single-Probe Zero-Degree
In the single-probe zero-degree method, the axial (d 0 ax ), lateral (d 0 lat ), and elevational (d 0 el ) inter-frame displacements estimates (d) based on either A 0 or B 0 were used as basis for strain estimation.Please note that the superscript denotes the displacement estimation direction relative to the corresponding probe.

Single-Probe Compounding
For the single-probe compounding method, the axial displacement (d 0 ax ), together with the lateral inter-frame displacement estimated by displacement compounding, served as basis for strain estimation.In displacement compounding, both angular axial displacement components (d α ax , d −α ax ) obtained using the angled transmissions for each probe (A −α , A α for probe A and B −α , B α for probe B) were used to derive the lateral displacement d 0 comp_lat : Combining the two axial displacements into a single lateral displacement has shown to be beneficial, since a lateral displacement estimate based on the phase information of two angled acquisitions yields a more accurate estimate than the lateral displacement estimation originating from zero-degree acquisitions [18].To determine the elevational displacement component for the single-probe compounding approach, the displacements d 0 ele , d α ele , and d −α ele were combined to a single estimate by means of the median displacement.Please note that these displacements are perpendicular to the elevational direction.

Dual-Probe Reference
The axial (d 0 ax ) displacement estimate used for this method equals the displacement estimate used for the single-probe zero-degree and single-probe compounding method.For the lateral displacement component (d 0 lat ), the axials A 0 & B 0 were registered such that they could be used as RF data originating from A 90 & B 90 , (see Figure 4 last column) preserving phase information in the lateral direction.This method combines 0 respectively.Since this method uses solely zero-degree axial displacement results, it is expected to provide the most accurate displacement estimates in the transverse plane, and therefore it will serve as reference in this work.For the elevational displacement, d 0 el was used.

Regularization Tracking and Segmentation
To accumulate the filtered inter-frame 3D displacements over the complete pressure cycle, displacement tracking was performed, with the end-diastolic frame as reference frame.Hereto axial, lateral, and elevational displacement estimates were accumulated spatially over time using 3D linear interpolation.Additional 3D median filtering was performed on the accumulated displacements, with kernel sizes listed in Table 2.
Finally, a Region of Interest (ROI) at the end diastolic frame was defined by manually contouring the lumen and wall using a composition of all 6 registered B-mode images for all 120 slices in elevational direction.Tracking and strain estimation was only performed for this ROI, representing the vessel wall and its direct surrounding to minimize the computational load.

Strain Estimation
For all three methods, strain estimation was performed using 3D least-squares strain estimation [43] in the axial, lateral, and elevational direction.The applied kernel sizes are listed in Table 2 and were empirically optimized for the reference method, after which they were kept constant for the other two methods.For the calculation of the 3D principal components, an eigenvalue decomposition was performed of the 3 × 3 strain tensor consisting of the normal strains and the shear strain components.At last, the 3 principal components were arranged from maximum to minimum values with corresponding eigenvectors.The maximum and minimum principal strains represent the maximum tensile strain component and the minimum compressive strain component, respectively.

Evaluation
In summary, Figure 6 shows a flowchart indicating the processing steps described in the previous sections from raw ultrasound acquisition, up to principal strain estimation for the zero-degree method, the compounding method, and the reference method.Figure 6.Flowchart indicating the in-line processes performed from the acquisitions to strain estimation to compare the single-probe compounding versus the single-probe zero-degree performance in terms of minimal (compressive) and maximum (tensile) 3D principal strains.Please note that this flowchart describes the situation with respect to probe A. A similar flowchart with respect to probe B was also applied and can be obtained by interchanging A and B.
For performance evaluation, the vertical, horizontal, and longitudinal strain results for both probes were quantitatively and qualitatively reported.The final performance evaluation was done by Root Mean Squared Error (RMSE) analysis of the minimum (compressive) and maximum (tensile) principal strains over time.

Results
Figure 7 shows the vertical and horizontal strain results at peak systolic phase (t = 0.3 s) for six distinct elevational slices for the three different methods.The reference results based on the axial strain as obtained with probe A and the axial strain obtained with probe B are presented in the left column.The compounding results and zero-degree single probe results are shown in the center and right column, respectively.
As indicated by the axial probe direction arrows and described in the methods section, the reference results originate from displacement estimations using solely information along the ultrasound beam direction for both probes.Clear resemblance is shown between the reference strains and the compounded strains, although the strains obtained by the compounding method appear noisier in some regions.This in contrast to the zero-degree lateral component, which depicts a highly non-uniform strain distribution with many strain clipping regions delivering unreliable lateral strain distributions.
Figure 8 depicts the elevational strain distributions at peak systolic pressure, for the reference, compounding, and zero-degree method for both probes.All strain distributions show compressive and tensile strains along the longitudinal direction of the vessel wall at the same longitudinal positions.Strains of up to 10% are observed at the external carotid artery for both probe orientations.Intra-probe strain distributions indicate a high similarity between the different methods.Inter-probe comparison, however, reveals a rotation in the strain distribution orientation at the inlet side (communis) and the outlet side (external CA).  Figure 9 shows the strain distribution over time for one cardiac cycle for the vertical, horizontal, and longitudinal direction for the reference, compounding and zero degree methods for probe A and B. For the vertical strains of probe A and the horizontal strain of probe B, all lines almost perfectly coincide.The 5th and 95th percentile horizontal strain of the zero-degree method show extreme strain values of <−15% and >20% at peak systolic pressure compared to the reference (~−5%, ~8%) strain, respectively.The compounding method outperforms the zero-degree method, which shows a high resemblance to the reference strain distribution for the complete cardiac cycle.For the longitudinal strain distribution, all three methods perform almost equally.Similar trends can be observed for probe B with respect to the reference of probe A. Table 3 lists the absolute median strain and the 5th and 95th percentile strain values for both probes to indicate the strain range at two different time points in the cardiac cycle.For the vertical component of probe A and the horizontal component of probe B, the percentile strains shows almost no difference at maximum systolic pressure for all methods.Please note that these strains originate from axial (Figure 5) displacement estimations for both probes.Also, small residual median strains were measured for these corresponding orientations.For the vertical strain of probe B and horizontal strain of probe A (i.e., lateral direction for each probe, percentile strain values are formatted in bold in Table 3), an increased strain range was found for the compounding technique.However, this range was much more increased for the zero-degree technique.Both strain graphs (Figure 9) show a divergence at t > 0.8 s of the strain range.This divergence is observed also for the longitudinal strain range for the complete pressure cycle leading to almost equal residual strain ranges in between ~−6% and ~11% for all methods and both probes.
Eigenvector decomposition of a 3 × 3 strain tensor provides 3 independent eigenvalues and eigenvectors that are perpendicular to each other.Figure 10 shows the minimum (bottom row) and maximum (top row) principal components, representing the most dominant compressive and tensile deformations as estimated by probe A at peak systolic pressure.For a dynamic development of these components for the complete cardiac cycle, please see video S2 (3D_Principal_MinMax_Strains.mp4).For visualization purposes, principal strains (color) and the corresponding eigenvectors (white arrows) are plotted for three distinct longitudinal positions.The top slice is situated in the common carotid artery.The middle slice is at the bifurcation, and the bottom slice shows the internal and external carotid arteries.In general, the tensile strain distribution is clearly oriented in the circumferential direction of the vessel wall with exception of a region in the bifurcation slice and one in the internal carotid artery slice, in which the orientation of the eigenvectors indicates a strong longitudinal preference.
The compressive strain distribution is mostly oriented in the radial direction and is less dominant than the tensile strain distribution.High tensile strains are observed consistently for all methods at the flow divider of the bifurcation and at ten and two o'clock of the common carotid artery.
Smooth values and eigenvectors can be appreciated for the reference method, indicating high tensile strains at the flow divider.Many outliers of the strain values a chaotic distribution of eigenvectors can be observed for the zero-degree method, clearly indicating the disability in providing consistent and trustworthy results.Especially in the regions in which the reference principal components point in the horizontal direction, erroneous values manifest for the zero-degree method.For the compounding method, the errors in these regions are far less destructive.Although the compounding technique is not perfect, it does yield strain results with less noise compared with the zero-degree method for the complete cardiac cycle, which shows erroneous results primarily in the regions in which the direction of the eigenvectors coincides with the lateral direction for each probe, which for probe A corresponds to the horizontal direction depicted in Figure 10.For the compounding method, the errors in these regions are far less destructive.
Figure 11 shows the RMSE of the minimal (left) and maximum (right) principal strain over time for a single cardiac cycle.The 5th-95th percentile strain ranges for both probes were depicted in the background to relate the RMSE to the overall strain development.In general, the RMSE for the compounding methods is smaller for the complete cycle in comparison with the zero-degree method.However, an up rise of the RMSE and the reference strain range occurs at t > ~0.8 s, which indicates a progressive strain error in the process of strain accumulation for all methods.The RMSE of the compressive and tensile strain components for the complete cardiac cycle for the zero-degree and compounding method for both probes.Please note the double axis for the RMSE minimal principal strain, in which the RMSE is always positive; the 5th-95th percentiles strain range is negative.
The compressive strain RMSE values for probe A and B using only the zero-degree angle at peak systole are 6.6% and 5.9% versus 2.7% and 3.1% for the compounding method.Also, the tensile strain RMSE shows elevated values for the zero-degree method of 14.3% and 12.6%, in comparison to 9.3% and 6.2% for the compounding method.
Looking at the residual compressive strain RMSE, values of 5.4% and 4.1% versus 1.5% and 1.5% were found for the zero-degree method and the compounding method, respectively.The same trend holds for the residual tensile strain RMSE, which was 5.8% and 6.4% versus 1.8% and 2.0%.
The mean RMSE values of the compressive principal strain over the complete cardiac cycle for the zero-degree method are 4.4% and 3.8%, respectively, versus 1.7% and 1.9% for the compounding method.For the tensile principal strain, 6.9% and 6.8% versus 3.3% and 4.0% were found.Therefore, the compounding method clearly outperforms the zero-degree method in terms of overall tensile and compressive principal strain over time.

Discussion
In this study, compressive and tensile 3D principal strain estimation in a pulsating carotid artery bifurcation phantom was performed using a compounding and zero-degree displacement estimation method in a dual probe setup.Previous experimental studies have been performed using straight-forward phantoms (i.e., single tube shaped phantoms) and simulations based on 2D and 3D compound strain imaging [19,23] showing a significant improvement compared to a single angle based strain estimation.Multiple 3D phantoms of the carotid artery have been developed lately [44,45].However, to our knowledge, no other paper has reported 3D principal strain estimation in a complex geometry like the CA bifurcation.The principal finding of this study is that 3D compound strain estimation outperforms single angle zero degree 3D strain estimation, even in a complex geometry, when embedded in surrounding tissue.Compared to other methods, this method is the first showing volumetric 3D strain tensor information for vascular applications retrieved in a non-invasive manner.
The developed dual probe setup approach, which provided a reference for the estimated deformations, is, to our knowledge, also novel and circumvents requiring exact knowledge of mechanical properties and geometry of the tissue and of boundary conditions, which is normally needed when deriving a finite element, model-based ground truth.Therefore, this setup might also be used in future studies to provide reference deformations for in vitro studies using excised tissue in which the mechanical properties are also not known.
To derive the reference strains from the RF data obtained with the dual-probe setup, the RF data acquired with both probes needs to be spatially registered.Since the angle of the probes is set to 90 degrees, only the translation should be resolved, which was performed using a feature-based registration using coherently compounded B-mode images.Although this might seem to be a simple registration, the quality of the registration depends on the choice of registration metric and also on the speed of sound and the assumption of a homogenous speed of sound throughout the phantom.In our case, setting the sound speed to 1540 m/s provided the most optimal registration, which is in line with reported speed of sound values for PVA [37].
In this work, a patient-specific model of the carotid artery at the bifurcation was used to construct a PVA-based phantom.Since the surrounding tissue affects the mechanical behavior of the carotid artery wall [46], the CA was embedded in surrounding PVA instead of water.The addition of surrounding tissue is not only important for obtaining realistic motion and deformation of the carotid artery segment but also resembles realistic in vivo acquisitions, since deterioration of the strain estimates due to side and grating lobes will be present.What was not accounted for in the current phantoms were heterogeneities within the vessel wall and plaque.Especially, the presence of calcifications and induced shadowing is expected to affect strain estimation accuracy.
In the experiment, the cardiac pressure cycle was enforced with 60 bpm.This resulted in a cyclic reproducible deformation of the vessel phantom.In the process of accumulation of inter-frame displacements, inaccuracies and interpolation errors led to an overall progressive displacement error over time, which explains the observed residual strain for all methods at the end of the pressure cycle (Figure 9).The longitudinal progression of the strain distribution (5th-95th percentile strain) over time shows a divergent trend for the full cardiac cycle.Since a ground truth strain distribution in the longitudinal direction is absent, this indicates the lack of a cyclic strain progression for each cardiac cycle.The elevational ultrasound beam width is one order of magnitude larger than the step-size in elevational direction.This results in redundant ultrasound RF data in the elevational direction.In a previous study, we showed the ability to resolve sub-plane elevational displacement estimation [23], but still it remains the direction with the relatively lowest displacement estimation resolution.Furthermore, the geometrical shape of the lumen in combination with this elevational beam overlap between slices might also lead to displacement estimation errors.Inter-slice geometrical lumen changes might dominate the normalized cross-correlation-process, which is supposed to detect inter-frame displacement changes.
Three methods for each probe should yield similar results for the strain estimation in the elevational direction.Figure 8 confirms that, in general, this is the case, since equal strain distributions between the zero-degree and compounding based methods are observed.Nevertheless, local strain differences were found between probe A and B. These were found especially at the inlet side and outlet side, and might be due to different displacement estimation kernel orientations and filter effects.
To calculate the 3D strain components, a 3D Least Square (LSQ) strain estimator was applied.The kernel size determines the tradeoff between outlier suppression and true value detection.Large kernels sizes result in smooth global strain estimates but also in underestimation of the local strain differences.Especially at the edges, underestimation might occur, since the displacement gradient of the true deformation in the vessel wall is the highest at the vessel wall-lumen transition.To compare the performance across the methods, the LSQ kernel size for the single probe compounding and zero degree method was kept equal to that of the reference method.Application of larger kernel sizes for these methods might be favorable to smoothen results without losing essential strain map features.This will be investigated in future studies.
Further research needs to be carried out to fully translate the single probe compounding technique into the clinic.The practical setup needs to be supportive to actually acquire multi-slice US data from the CA in vivo.Due to the curved neck and jaw skin surface, maintaining acoustic contact will be challenging during the complete acquisition range in the elevational direction.A possible solution to this might be to create a thick, pre-shaped layer of ultrasound transparent material, ensuring constant acoustic contact between the transducer and skin surface.
Another challenge for in vivo acquisition is related to the data acquisition.The trigger originating from the pump driver at the start of each systolic phase initiated data acquisition for the complete cardiac cycle.Since the acquisition took 1 second at each position, thus several minutes for the entire acquisition, the phantom motion had to be very reproducible.In vivo, this multi-slice technique might become more difficult to achieve due to heart rate and blood pressure variability and motion artifacts due to, for example, breathing.However, it is already known that heart rate variability in general only affects the diastolic phase and not the systolic phase of the pressure cycle [47], and only one phase of the pressure cycle is required to obtain maximum strain contrast.Thus, for in vivo application we recommend to apply this technique on the systolic phase while the subject is in resting position to minimize blood pressure variations.ECG-triggering will then need to be incorporated to enable synchronization in vivo, and probably the head and neck region should be fixated to minimize motion artefact.
Finally, for in vivo application a reduction in overall acquisition time is favored.The acquisition time can already be halved compared to the acquisition sequence used in this paper, because, as aforementioned, it is not necessary to acquire data for the diastolic phase.With the development of matrix probes with a sufficient footprint, pitch, and number of elements, it might even become possible to acquire volumetric ultrasound and strain data within one cardiac cycle.These probes will also enable displacement compounding in the elevational direction to increase the elevational displacement estimation accuracy.Inherent to the multi-slice acquisition protocol, no beam steering could be applied in the current implementation.

Figure 1 .
Figure 1.(a) 3D CAD-based mold of a realistic patient-specific carotid artery and (b) Polyvinyl alcohol phantom casted in a surrounding body.

Figure 2 .
Figure 2. (a) A schematic drawing of the image acquisition setup in which two probes were fixated orthogonally.The position of the probes was set such that the carotid bifurcation was at a realistic depth and all steered and non-steered insonified regions by the plane wave transmissions encompassed the carotid artery.The phantom was connected to a pressurized realistic dynamic flow circuit; (b) The transducer fixation clamp was connected to a translation stage enabling elevational automated movement in a step size of 0.1 mm.
depicts the temporal acquisition scheme showing subsequent acquisition series.Each series consists of plane wave acquisitions at alternating angles[23] of −19.47 • , 19.47 • , and 0 • for probe A and 0 • and −19.47 • and 19.47 • for probe B. The order in which these acquisitions took place was chosen such that the time between combined acquisitions in post processing was minimized.

Figure 3 .
Figure 3. Schematic of the temporal acquisition protocol applied using an interleaved dual probe setup.

Figure 4 .
Figure 4.The first column depicts the grid configurations for both probes in which the overlay of all steered and non-steered grids results in a displacement grid consisting of coinciding sample-locations.Columns 2, 3, and 4 show the B-mode representations of the actual transmissions performed subsequently overlaid with a displacement grid and a single line indicating the direction of the reconstructed ultrasound beam in receive.The schematic drawing of the aperture indicates the +45 • and −45 • angles of insonification relative to the phantom for probe A and B, respectively.The last column shows the RF data registration in which the axial direction of the probes serves as a native lateral ultrasound transmission of the opposite probe.

Figure 5 .
Figure 5.The axis definition and nomenclature for the probe related directions: axial, lateral, and elevational.Vertical, horizontal, and longitudinal directions are related to the phantom orientation.

Figure 7 .
Figure 7. Vertical and horizontal strain estimates of the bifurcation in 3D at maximum systolic pressure.Please note that for visualization purposes, only 6 distinctive transverse slices are depicted at a 3 times larger elevational scale.

Figure 8 .
Figure 8. Longitudinal strain estimates originating from both probes at peak systolic pressure.

Figure 9 .
Figure 9.The vertical, horizontal, and longitudinal strain range (median and 5th-95th percentile) for the complete cardiac cycle for the reference, compounding, and zero-degree method for probe A and B.

Figure 10 .
Figure 10.Maximum (tensile) and minimum (compressive) principal strain distributions for probe A for the reference, compounding, and zero-degree method in the bifurcation carotid artery region.Principal strains are visualized at the level of the COmmon carotid artery (CO), bifurcation, and distal from the bifurcation, in which the Internal Carotid Artery (ICA) and the External Carotid Artery (ECA) are present.Please note that the lengths of the eigenvectors were scaled to their corresponding eigenvalue.

Figure 11 .
Figure 11.The RMSE of the compressive and tensile strain components for the complete cardiac cycle for the zero-degree and compounding method for both probes.Please note the double axis for the RMSE minimal principal strain, in which the RMSE is always positive; the 5th-95th percentiles strain range is negative.

Table 1 .
Resolution of interlocked steered and non-steered grids.

Table 3 .
Vertical, horizontal, and longitudinal 5th-95th percentile strain values for both probes and all methods.